151: Impedance calculation

(source code)

Same as Example150, but with new and more generic way of passing the parameter.

Impedance calculation for

C ut - (D ux)_x + Ru = 0 in (0,1) u(0,t)=1 + exp(iωt) u(1,t)=0

Measurement: I(t)= D u_x(1,t)

Steady state:

• (D u0x)x + Ru0 = 0

u0(0,t)=1 u0(1,t)=0

Small signal ansatz for ω

u(x,t)= u0(x)+ ua(x) exp(iωt)

iωC ua - (D uax)x + R u_a =0 ua(0)=1 ua(1)=0

module Example151_Impedance1Dusing Printfusing VoronoiFVMusing ExtendableGrids: geomspace, simplexgridusing GridVisualizeusing OrdinaryDiffEqSDIRKfunction main(;        nref = 0, Plotter = nothing, verbose = false,        unknown_storage = :sparse, assembly = :edgewise,        time_embedding = :none,        L = 1.0, R = 1.0, D = 1.0, C = 1.0,        ω0 = 1.0e-3, ω1 = 5.0e1    )

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[1, 2])        return nothing    end    storage = function (f, u, node, data)        f[1] = data.C * u[1]        return nothing    end    reaction = function (f, u, node, data)        f[1] = data.R * u[1]        return nothing    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; region = excited_bc, value = p[1])        boundary_dirichlet!(f, u, node; region = meas_bc, value = 0.0)        return nothing    end

Create discrete system and enable species

    sys = VoronoiFVM.System(        grid; unknown_storage = unknown_storage,        data = data,        flux = flux,        storage = storage,        reaction = reaction,        bcondition = bc,        nparams = 1,        species = 1, assembly = assembly    )

Create test functions for current measurement

    factory = TestFunctionFactory(sys)    measurement_testfunction = testfunction(factory, [excited_bc], [meas_bc])    tend = 1.0    if time_embedding == :builtin        tsol = solve(sys; inival = 0.0, params = [1.0], times = (0.0, tend), force_first_step = true)        steadystate = tsol.u[end]    elseif time_embedding == :ordinarydiffeq        inival = unknowns(sys, inival = 0)        problem = ODEProblem(sys, inival, (0, tend); params = [1.0])        odesol = solve(problem, ImplicitEuler())        tsol = reshape(odesol, sys)        steadystate = tsol.u[end]    elseif time_embedding == :none        steadystate = solve(sys; inival = 0.0, params = [1.0])    else        error("time_embedding must be one of :builtin, :ordinarydiffeq, :none")    end    function meas_stdy(meas, U)        u = reshape(U, sys)        meas[1] = -VoronoiFVM.integrate_stdy(sys, measurement_testfunction, u, params = [1.0])[excited_spec]        return nothing    end    function meas_tran(meas, U)        u = reshape(U, sys)        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)

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)    while ω < ω1

solve 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)

record exact solution

        iω = 1im * ω        z = sqrt(iω * data.C / data.D + data.R / data.D)        eplus = exp(z * L)        eminus = exp(-z * L)        IxL = 2.0 * data.D * z / (eplus - eminus)        push!(allIxL, 1 / IxL)

increase omega

        ω = ω * 1.1    end    vis = GridVisualizer(; Plotter = Plotter)    scalarplot!(        vis, real(allIxL), imag(allIxL); label = "exact", color = :red,        linestyle = :dot    )    scalarplot!(        vis, real(allIL), imag(allIL); label = "calc", show = true, clear = false,        color = :blue, linestyle = :solid    )    return sum(allIL)endusing Testfunction runtests()    testval = 57.92710286186797 + 23.163945443946027im    for unknown_storage in (:sparse, :dense)        for assembly in (:edgewise, :cellwise)            for time_embedding in (:none, :builtin, :ordinarydiffeq)                @test main(; unknown_storage, assembly, time_embedding) ≈ testval            end        end    end    return nothingendend

This page was generated using Literate.jl.