GitHub

Postprocessing

Plotting

Plotting can be performed using the package GridVisualize.jl. This package extends the API with a couple of methods for systems:

GridVisualize.scalarplot — Function
scalarplot(
    sys::VoronoiFVM.AbstractSystem,
    sol::AbstractMatrix;
    species,
    scale,
    kwargs...
) -> Any

Plot one species from solution

source
scalarplot(
    sys::VoronoiFVM.AbstractSystem,
    sol::VoronoiFVM.AbstractTransientSolution;
    species,
    scale,
    tscale,
    kwargs...
) -> Any

Plot one species from transient solution

source
GridVisualize.scalarplot! — Function
scalarplot!(
    vis,
    sys::VoronoiFVM.AbstractSystem,
    sol::AbstractMatrix;
    species,
    scale,
    kwargs...
) -> Any

Plot one species from solution

source
scalarplot!(
    vis,
    sys::VoronoiFVM.AbstractSystem,
    sol::VoronoiFVM.AbstractTransientSolution;
    species,
    scale,
    tscale,
    tlabel,
    kwargs...
) -> Any

Plot one species from transient solution

source
VoronoiFVM.plothistory — Function
plothistory(tsol; 
            plots=[:timestepsizes,
                   :timestepupdates,
                   :newtonsteps,
                   :newtonupdates], 
            size=(700,600), 
            logmin=1.0e-20,
            Plotter=GridVisualize.default_plotter(),
            kwargs...)

Plot solution history stored in tsol. The plots argument allows to choose which plots are shown.

source

Grid verification

VoronoiFVM.nondelaunay — Function
nondelaunay(grid;tol=1.0e-14, verbose=true)

Return non-Delaunay edges. Returns a vector of tuples: Each tuple consists of (node1, node2, edge factor, region)

If the vector has length 0, the grid is boundary conforming Delaunay with respect to each cell region. This means that up to tol, all edge form factors are nonnegative.

source

Norms & volumes

LinearAlgebra.norm — Function
norm(system, u)
norm(system, u, p)

Calculate Euklidean norm of the degree of freedom vector.

source
norm(system, u)
norm(system, u, p)

Calculate Euklidean norm of the degree of freedom vector.

source
VoronoiFVM.lpnorm — Function
lpnorm(sys, u, p)
lpnorm(sys, u, p, species_weights)

Calculate weighted discrete $L^p$ norm of a solution vector.

source
VoronoiFVM.l2norm — Function
l2norm(sys, u)
l2norm(sys, u, species_weights)

Calculate weighted discrete $L^2(\Omega)$ norm of a solution vector.

source
VoronoiFVM.w1pseminorm — Function
w1pseminorm(sys, u, p)
w1pseminorm(sys, u, p, species_weights)

Calculate weighted discrete $W^{1,p}(\Omega)$ seminorm of a solution vector.

source
VoronoiFVM.h1seminorm — Function
h1seminorm(sys, u)
h1seminorm(sys, u, species_weights)

Calculate weighted discrete $H^1(\Omega)$ seminorm of a solution vector.

source
VoronoiFVM.w1pnorm — Function
w1pnorm(sys, u, p)
w1pnorm(sys, u, p, species_weights)

Calculate weighted discrete $W^{1,p}(\Omega)$ norm of a solution vector.

source
VoronoiFVM.h1norm — Function
h1norm(sys, u)
h1norm(sys, u, species_weights)

Calculate weighted discrete $H^1(\Omega)$ norm of a solution vector.

source
VoronoiFVM.lpw1pseminorm — Function
lpw1pseminorm(sys, u, p)
lpw1pseminorm(sys, u, p, species_weights)

Calculate weighted discrete $L^p([0,T];W^{1,p}(\Omega))$ seminorm of a transient solution.

source
VoronoiFVM.l2h1seminorm — Function
l2h1seminorm(sys, u)
l2h1seminorm(sys, u, species_weights)

Calculate weighted discrete $L^2([0,T];H^1(\Omega))$ seminorm of a transient solution.

source
VoronoiFVM.lpw1pnorm — Function
lpw1pnorm(sys, u, p)
lpw1pnorm(sys, u, p, species_weights)

Calculate weighted discrete $L^p([0,T];W^{1,p}(\Omega))$ norm of a transient solution.

source
VoronoiFVM.l2h1norm — Function
l2h1norm(sys, u)
l2h1norm(sys, u, species_weights)

Calculate weighted discrete $L^2([0,T];H^1(\Omega))$ norm of a transient solution.

source

Solution integrals

VoronoiFVM.integrate — Method
integrate(system,F,U; boundary=false, data=system.physics.data)

Integrate node function (same signature as reaction or storage) F of solution vector U region-wise over domain or boundary. The result is an nspec x nregion matrix $\int_{\Omega_j} F_i(U)\, dx$ or $\int_{\Gamma_j} F_i(U)\, ds$.

source
VoronoiFVM.integrate — Method
integrate(system,F, Ut; boundary=false, data=system.physics.data)

Integrate transient solution vector region-wise over domain or boundary. The result is an nspec x nregion x nsteps DiffEqArray $\int_{\Omega_j} F_i(U_k)\, dx$ or $\int_{\Gamma_j} F_i(U_k)\, ds$.

source
VoronoiFVM.integrate — Method
integrate(system,U; boundary=false)

Integrate solution vector region-wise over domain or boundary. The result is an nspec x nregion matrix $\int_{\Omega_j} U_i\, dx$ or $\int_{\Gamma_j} U_i\, ds$.

source
VoronoiFVM.edgeintegrate — Function
edgeintegrate(system,F,U; boundary=false)

Integrate edge function (same signature as flux function) F of solution vector region-wise over domain or boundary. The result is an nspec x nregion matrix.

source

Boundary fluxes

Continuum level motivation

For the calculation of boundary fluxes, we use a method inspired by H. Gajewski, WIAS Report No 6, 1993: "The regularity of the (weak) solutions of the drift-diffusion equations, which are guaranteed by existence statements, is not sufficient to justify a naive calculation of the currents as boundary integrals. Rather, it has also proven to be practical and necessary to reduce contact currents to volume integrals in accordance with the definition of weak solutions using suitable test functions."

We start with the continuous case for a problem

\[ \partial_t s(u) + \nabla\cdot \vec j(u) + r(u) - f =0 \]

defined in a domain $\Omega$ with a boundary $Γ=⋃_{i\in B}Γ_i$ subdivided into non-overlapping parts where the set of boundary regions $B$ is subdivided into $B=B_N\cup B_0 \cup B_1$. Assume that $\vec j\cdot \vec n=0$ on $\Gamma_i$ for $i\in B_N$.

Let $T(x)$ be a smooth test function such that

  • $\nabla T \cdot \vec n|_{Γ_i}= 0$ for $i\in B_N$,
  • $T|_{Γ_i} = 0$ for $i\in B_0$,
  • $T|_{Γ_i} = 1$ for $i\in B_1$.

Here, we obtain it by solving the Laplace equation

\[ -\Delta T =0 \quad \text{in}\; \Omega\]

To obtain as a quantity of interest the flux $Q$ through the boundaries $\Gamma_i$ for $i\in B_1$ of, calculate:

\[ \begin{aligned} Q=Q(t)&=\sum_{i \in B_1}\int\limits_{\Gamma_i} T\vec j(u)\cdot\vec n \,d\gamma\\ &=\sum_{i \in B_1}\int\limits_{\Gamma_i} T\vec j(u)\cdot\vec n \,d\gamma + \sum_{i \in B_N}\int\limits_{\Gamma_i} T\vec j(u)\cdot\vec n \,d\gamma + \sum_{i \in B_0}\int\limits_{\Gamma_l} T\vec j(u)\cdot\vec n \,d\gamma\\ &=\int\limits_{\Gamma} T\vec j(u)\cdot\vec n \,d\gamma\\ &=\int\limits_{\Omega}\nabla\cdot(T\vec j(u))\,d\omega\\ &=\int\limits_{\Omega}\nabla T \cdot \vec{j(u)} \,d\omega + \int\limits_{\Omega}T \nabla\cdot \vec{j(u)} \,d\omega\\ &=\int\limits_{\Omega}\nabla T \cdot \vec{j}(u) \,d\omega - \int\limits_{\Omega}T r(u) \,d\omega + \int\limits_{\Omega}T f \,d\omega - \int\limits_{\Omega}T \partial_ts(u) \,d\omega \end{aligned}\]

Discrete approximations of the integrals

The discete versions of these integrals for evaluation at $t=t^n$ are as follows:

\[ \begin{aligned} \int\limits_{\Omega}T r(u) \,d\omega &\approx \sum_{k\in N} T_kr(u^n_k)|\omega_k| =: I_{func}(r, T,u^n)\\ \int\limits_{\Omega}T f \,d\omega &\approx \sum_{k\in N} T_k f_k|\omega_k| =: I_{src}(f, T)\\ \int\limits_{\Omega}T \partial_ts(u) \,d\omega &\approx \frac{ I_{func}(s,T,u^n) -I_{func}(s,T,u^{n-1})}{t^n-t^{n-1}}\\ \int\limits_{\Omega}\nabla T \cdot \vec{j}(u^n) d\omega &\approx \sum\limits_{\stackrel{k,l}{\partial \omega_k \cup \partial \omega_l\neq\emptyset}} (T_{k}-T_{l})\frac{\sigma_{kl}}{h_{kl}}g(u^n_k, u^n_l) =: I_{flux}(j,T,u) \end{aligned}\]

General flux calculation API

VoronoiFVM.TestFunctionFactory — Type
mutable struct TestFunctionFactory{Tu, Tv}

Data structure containing DenseSystem used to calculate test functions for boundary flux calculations.

Type parameters:

  • Tu: value type of test functions
  • Tv: Default value type of system

Fields:

  • system::VoronoiFVM.AbstractSystem{Tv} where Tv: Original system
  • state::VoronoiFVM.SystemState{Tu, Tp, TMatrix, TGenericMatrix, TSolArray} where {Tu, Tp, TMatrix<:AbstractMatrix{Tu}, TGenericMatrix, TSolArray<:AbstractMatrix{Tu}}: Test function system state
  • control::SolverControl: Solver control
source
VoronoiFVM.testfunction — Function
testfunction(factory::TestFunctionFactory, B_0, B_1)

Create a test function $T$ defined in the domain $Ω$ given by the system behind the factory object. Assume that the boundary $Γ=∂Ω=⋃_{i∈ B}Γ_i$ is subdivided into non-overlapping parts where the set of boundary regions $B$ is subdivided into $B=B_N∪ B_0 ∪ B_1$. Create $T$ by solving

\[ -Δ T =0 \quad \text{in}\; Ω\]

such that

  • $∇ T ⋅ \vec n|_{Γ_i}= 0$ for $i∈ B_N$
  • $T|_{Γ_i} = 0$ for $i∈ B_0$,
  • $T|_{Γ_i} = 1$ for $i∈ B_1$.

Returns a vector representing T.

source
VoronoiFVM.integrate — Method
 integrate(system, T, U::AbstractMatrix)

Calculate test function integral for a steady state solution $u$ $∫_Γ T \vec j(u) ⋅ \vec n dω ≈ I_{flux}(j, T,u)-I_{func}(r,T,u) + I_{src}(f,T)$, see the definition of test function integral contributions.

The result is a nspec vector giving one value of the integral for each species.

source
VoronoiFVM.integrate — Method
integrate(system,T, U::AbstractTransientSolution; rate=true, params, data)

Calculate test function integrals for the transient solution $∫_Γ T \vec j(u^n) ⋅ \vec n dω ≈ I_{flux}(j,T,u^n)-I_{func}(r,T,u^n)-\frac{I_{func}(s,T,u^n)-I_{func}(s,T,u^{n-1})}{t^n-t^{n-1}} + I_{src}(f,T)$ for each time interval $(t^{n-1}, t^n)$, see the definition of test function integral contributions.

Keyword arguments:

  • params: vector of parameters used to calculate U. Default: [].
  • data: user data used to calculate U. Default: data(system)
  • rate: If rate=true (default), calculate the flow rate (per second) through the corresponding boundary. Otherwise, calculate the absolute amount per time interval.

The result is a nspec x (nsteps-1) DiffEqArray.

source

Special purpose cases

VoronoiFVM.integrate — Method
 integrate(system, T, U::AbstractMatrix, Uold::AbstractMatrix, Δt; kwargs...)

Calculate test function integrals for the implicit Euler time step results $u^n≡$ U and $u^{n-1}≡$ Uold as $∫_Γ T \vec j(u^n) ⋅ \vec n dω ≈ I_{flux}(j,T,u^n)-I_{func}(r,T,u^n)-\frac{I_{func}(s,T,u^n)-I_{func}(s,T,u^{n-1})}{Δt} + I_{src}(f,T)$, see the definition of test function integral contributions.

Keyword arguments:

  • params: vector of parameters used to calculate U. Default: [].
  • data: user data used to calculate U. Default: data(system)

The result is a nspec vector giving one value of the integral for each species.

source
VoronoiFVM.integrate_TxSrc — Function
integrate_TxSrc(system, T, f!; kwargs...)

Calculate $I_{src}(f!,T)= ∫_Ω T⋅ f! dω$ for a test function T. The functionf!` shall have the same signature as a storage or reaction function.

source
VoronoiFVM.integrate_∇TxFlux — Function
  integrate_∇TxFlux(system, T, f!, U; kwargs...)

Calculate $I_{flux}(f!,T, U)=∫_Ω ∇T⋅ f!(U) dω$ for a test function T and unknown vector U. The function f! shall have the same signature as a flux function.

source
  integrate_∇TxFlux(system, T, U; kwargs...)

Calculate $I_{flux}(f!,T, U)=∫_Ω ∇T⋅ f!(U) dω$ for a test function T and unknown vector U, where flux! is the system flux.

source
VoronoiFVM.integrate_TxEdgefunc — Function
  integrate_TxEdgeFunc(system, T, f!, U; kwargs...)

Calculate $∫_Ω T⋅ f!(U) dω$ for a test function T and unknown vector U. The function f! shall have the same signature as a flux function, but is assumed to describe a reaction term given on edges.

source
VoronoiFVM.integrate_tran — Function
integrate_tran(system, T, U; kwargs...)

Calculate storage term contribution to test function integral. Used for impedance calculations.

source

Nodal flux reconstruction

VoronoiFVM.nodeflux — Function
nodeflux(system, F,U; data)

Reconstruction of an edge function as vector function on the nodes of the triangulation. The result can be seen as a piecewiesw linear vector function in the FEM space spanned by the discretization nodes exhibiting the flux density.

The reconstruction is based on the "magic formula" R. Eymard, T. Gallouët, R. Herbin, IMA Journal of Numerical Analysis (2006) 26, 326−353, Lemma 2.4 (also: hal.science/hal-00004840 ).

The return value is a dim x nspec x nnodes vector.

CAVEAT: there is a possible unsolved problem with the values at domain corners in the code. Please see any potential boundary artifacts as a manifestation of this issue and report them.

source
nodeflux(system, U; data)

Reconstruction of the edge flux as vector function on the nodes of the triangulation. See nodeflux(system,F,U;data).

The flux of species i can e.g. plotted via GridVisualize.vectorplot.

Example:

    ispec=3
    vis=GridVisualizer(Plotter=Plotter)
    scalarplot!(vis,grid,solution[ispec,:],clear=true,colormap=:summer)
    vectorplot!(vis,grid,nf[:,ispec,:],clear=false)
    reveal(vis)
source

Impedance calculatiom

Impedance calculation can be seen as a postprocessing step after the solution of the unexcited stationary system.

VoronoiFVM.ImpedanceSystem — Type
mutable struct ImpedanceSystem{Tv} <: VoronoiFVM.AbstractImpedanceSystem{Tv}

Concrete type for impedance system.

  • sysnzval::AbstractArray{Complex{Tv}, 1} where Tv: Nonzero pattern of time domain system matrix
  • storderiv::AbstractMatrix: Derivative of storage term
  • matrix::AbstractArray{Complex{Tv}, 2} where Tv: Complex matrix of impedance system
  • F::AbstractArray{Complex{Tv}, 2} where Tv: Right hand side of impedance system
  • U0::AbstractMatrix: Stationary state
source
VoronoiFVM.ImpedanceSystem — Method
ImpedanceSystem(system, U0, excited_spec, excited_bc)

Construct impedance system from time domain system sys and steady state solution U0 under the assumption of a periodic perturbation of species excited_spec at boundary excited_bc.

source
VoronoiFVM.freqdomain_impedance — Method
freqdomain_impedance(
    impedance_system,
    ω,
    U0,
    excited_spec,
    excited_bc,
    excited_bcval,
    dmeas_stdy,
    dmeas_tran
)

Calculate reciprocal value of impedance.

  • excitedspec,excitedbc,excited_bcval are ignored.
Warning

This is deprecated: use impedance.

source
VoronoiFVM.impedance — Method
impedance(impedance_system,ω, U0 ,
          excited_spec, excited_bc, excited_bcval,
           dmeas_stdy,
           dmeas_tran 
           )

Calculate impedance.

  • ω: frequency
  • U0: steady state slution
  • dmeas_stdy: Derivative of steady state part of measurement functional
  • dmeas_tran Derivative of transient part of the measurement functional
source
VoronoiFVM.measurement_derivative — Method
measurement_derivative(system, measurement_functional, U0)

Calculate the derivative of the scalar measurement functional at steady state U0

Usually, this functional is a test function integral. Initially, we assume that its value depends on all unknowns of the system.

source