same as example 150 but with two species and a different system. This time the exact solution is not known, thus we compare with finite-difference approximation of the impedance. PDE system on x∈(0,L): C∂ₜu₁ - ∂ₓ(D∂ₓ(u₁u₂)) + (Ru₁u₂ - u₂) = 0 C∂ₜu₂ - ∂ₓ(D∂ₓu₂) + Ru₁u₂ = 0 consistent with VoronoiFVM convention ∂ₜs(u)+∇⋅j(u)+r(u)=0 and diffusion flux j=-D∇(⋅). Boundary conditions used here: Dirichlet: u₂(0,t)=1, u₁(0,t)=excitation(t), u₁(L,t)=0 Neumann (natural, zero flux): D∂ₓu₂(L,t)=0
module Example152_Impedance_Multispeciesusing VoronoiFVMusing ExtendableGrids: geomspace, simplexgrid, num_nodesusing GridVisualizeusing OrdinaryDiffEqSDIRKusing Printffunction main(; nref = 0, Plotter = nothing, verbose = false, unknown_storage = :sparse, assembly = :edgewise, L = 1.0, R = 1.0, D = 1.0, C = 1.0, ω0 = 1.0e-3, ω1 = 5.0e1, ω_incfactor = 1.1, N_preliminary_periods::I = 2, Ndt::I = 300, fdtest::Bool = false ) where {I <: Integer} @assert N_preliminary_periods >= 0 "preliminary periods should be non-negative" @assert Ndt > 10 "Ndt should be at least 10 to have a reasonable sampling of the period"Create array which is refined close to 0
h0 = 0.005 / 2.0^nref h1 = 0.1 / 2.0^nref X = geomspace(0, L, h0, h1)Create discretization grid
grid = simplexgrid(X)Create and fill data
data = (R = R, D = D, C = C)Declare constitutive functions
flux = function (f, u, edge, data) f[1] = data.D * (u[1, 1] * u[2, 1] - u[1, 2] * u[2, 2]) return f[2] = data.D * (u[2, 1] - u[2, 2]) end storage = function (f, u, node, data) f[1] = data.C * u[1] return f[2] = data.C * u[2] end reaction = function (f, u, node, data) f[1] = data.R * u[1] * u[2] - u[2] return f[2] = data.R * u[2] * u[1] end excited_bc = 1 excited_bcval = 1.0 excited_spec = 1 meas_bc = 2 bc = function (f, u, node, data) p = parameters(u) boundary_dirichlet!(f, u, node; species = 2, region = 1, value = 1.0) boundary_dirichlet!(f, u, node; species = 1, region = excited_bc, value = p[1]) return boundary_dirichlet!(f, u, node; species = 1, region = meas_bc, value = 0.0) end sys = VoronoiFVM.System( grid; unknown_storage = unknown_storage, data = data, flux = flux, storage = storage, reaction = reaction, bcondition = bc, nparams = 1, assembly = assembly ) enable_species!(sys, 1, [1]) enable_species!(sys, 2, [1]) factory = TestFunctionFactory(sys) measurement_testfunction = testfunction(factory, [excited_bc], [meas_bc]) steadystate = solve(sys; inival = 1.0, params = [1.0]) function meas_stdy(meas, U) if !(typeof(U) <: AbstractMatrix) u = reshape(U, sys) else u = U end meas[1] = -VoronoiFVM.integrate_stdy(sys, measurement_testfunction, u, params = [1.0])[excited_spec] return nothing end function meas_tran(meas, U) if !(typeof(U) <: AbstractMatrix) u = reshape(U, sys) else u = U end meas[1] = -VoronoiFVM.integrate_tran(sys, measurement_testfunction, u, params = [1.0])[excited_spec] return nothing end dmeas_stdy = measurement_derivative(sys, meas_stdy, steadystate) dmeas_tran = measurement_derivative(sys, meas_tran, steadystate) meas_tran_ref = zeros(1) meas_stdy_ref = zeros(1) meas_cos = zeros(1) meas_sin = zeros(1) meas_tran(meas_tran_ref, steadystate) meas_stdy(meas_stdy_ref, steadystate)Create impedance system from steady state
isys = VoronoiFVM.ImpedanceSystem(sys, steadystate)Prepare recording of impedance results
allomega = zeros(0)for calculated data
allI0 = zeros(Complex{Float64}, 0) allIL = zeros(Complex{Float64}, 0)for exact data
allIx0 = zeros(Complex{Float64}, 0) allIxL = zeros(Complex{Float64}, 0) ω = ω0 UZ = unknowns(isys) outflux_ref = zeros(2) outflux_cos = zeros(2) outflux_sin = zeros(2) nnodes = num_nodes(grid) lastedge = (nnodes - 1):nnodes @views flux(outflux_ref, steadystate[:, lastedge], nothing, data) while ω < ω1solve impedance system
solve!(UZ, isys, ω)calculate approximate solution obtain measurement in frequency domain
IL = impedance(isys, ω, steadystate, dmeas_stdy, dmeas_tran)record approximate solution
push!(allomega, ω) push!(allIL, IL) if fdtestcompute reference using finite difference approximation
amplitude = 1.0e-6 data_perturbed = (R = R, D = D, C = C, ω = ω) #change boundary condition to reflect the perturbation bc_cos = function (f, u, node, data) p = parameters(u) boundary_dirichlet!(f, u, node; species = 2, region = 1, value = 1.0) boundary_dirichlet!(f, u, node; species = 1, region = excited_bc, value = 1 + amplitude * cos(data.ω * node.time)) return boundary_dirichlet!(f, u, node; species = 1, region = meas_bc, value = 0.0) end sys_cos = VoronoiFVM.System( grid; unknown_storage = unknown_storage, data = data_perturbed, flux = flux, storage = storage, reaction = reaction, bcondition = bc_cos, nparams = 0, #we no longer track derivative with respect to parameters assembly = assembly ) enable_species!(sys_cos, 1, [1]) enable_species!(sys_cos, 2, [1]) #same for the sine perturbation bc_sin = function (f, u, node, data) p = parameters(u) boundary_dirichlet!(f, u, node; species = 2, region = 1, value = 1.0) boundary_dirichlet!(f, u, node; species = 1, region = excited_bc, value = 1 + amplitude * sin(data.ω * node.time)) return boundary_dirichlet!(f, u, node; species = 1, region = meas_bc, value = 0.0) end sys_sin = VoronoiFVM.System( grid; unknown_storage = unknown_storage, data = data_perturbed, flux = flux, storage = storage, reaction = reaction, bcondition = bc_sin, nparams = 0, assembly = assembly ) enable_species!(sys_sin, 1, [1]) enable_species!(sys_sin, 2, [1]) dt = (2 * π / ω) / (Ndt - 1.0e-8) # without the perturbation we end up sometimes with one extra time step at the end. tend = (N_preliminary_periods + 1) * 2 * π / ωCompute a sufficiently long transient and evaluate the impedance on the last period.
tsol_cos = solve( sys_cos; inival = steadystate, times = (0.0, tend), force_first_step = true, control = VoronoiFVM.SolverControl(Δt_max = dt, Δt_min = dt, Δt = dt, Δu_opt = 1.0e10) ) tsol_sin = solve( sys_sin; inival = steadystate, times = (0.0, tend), force_first_step = true, control = VoronoiFVM.SolverControl(Δt_max = dt, Δt_min = dt, Δt = dt, Δu_opt = 1.0e10) ) @assert length(tsol_cos.t) >= Ndt "Need at least Ndt points to sample last period" @assert length(tsol_sin.t) == length(tsol_cos.t) "Cos and sin solutions should have the same time points" #and use the results to compute the impedance using finite difference approximation time_impedance = zeros(ComplexF64, Ndt) j_last_period = length(tsol_cos.t) - Ndt for i in 1:Ndt j = j_last_period + i time = tsol_cos.t[j] @assert isapprox(time, tsol_sin.t[j], rtol = 1.0e-5) u_cos = tsol_cos.u[j] u_sin = tsol_sin.u[j] #compute flux at the boundary for both solutions and subtract the reference flux to get the flux perturbation endcos = view(u_cos, :, lastedge) endsin = view(u_sin, :, lastedge) flux(outflux_cos, endcos, nothing, data) flux(outflux_sin, endsin, nothing, data) outflux_cos .-= outflux_ref outflux_sin .-= outflux_ref tau = 1 / (X[end] - X[end - 1]) time_impedance[i] = (outflux_cos[1] * tau + 1im * outflux_sin[1] * tau) / (amplitude * exp(1im * ω * time)) end IxL = length(time_impedance) / sum(time_impedance) if verbose ratio = IL / IxL @printf( "Finite difference approximation of impedance at ω = %10.5g: %10.5g%+10.5gi, calculated impedance: %10.5g%+10.5gi, ratio distance to one: %10.5g\n", ω, real(IxL), imag(IxL), real(IL), imag(IL), abs(ratio - 1.0) ) end push!(allIxL, IxL) endincrease omega
ω = ω * ω_incfactor end vis = GridVisualizer(; Plotter = Plotter, legend = :rt) if fdtest scalarplot!( vis, real(allIxL), imag(allIxL); label = "finite difference", color = :red, linestyle = :dot ) end scalarplot!( vis, real(allIL), imag(allIL); label = "calc", show = true, clear = false, color = :blue, linestyle = :solid ) if fdtest println("Ratio of calculated impedance to finite difference impedance, should be close to one: ") avg_ratio = sum(allIL ./ allIxL) / length(allomega) @printf("Minimum distance to one: %.3e, Average value: %.3g%+.3gi, Maximum distance to one: %.3e\n", minimum(abs.(allIL ./ allIxL .- 1)), real(avg_ratio), imag(avg_ratio), maximum(abs.(allIL ./ allIxL .- 1))) #@show minimum(abs.(allIL ./ allIxL)), sum(allIL ./ allIxL) / length(allomega), maximum(abs.(allIL ./ allIxL)) end return sum(allIL)endusing Testfunction runtests() testval = 50.960361928838 + 4.510584656768053im for unknown_storage in (:sparse, :dense) for assembly in (:edgewise, :cellwise) @test main(; unknown_storage, assembly) ≈ testval end end returnendend #end of moduleThis page was generated using Literate.jl.