Extracting integral data from Finite Volume solutions

After calculating solutions based on the finite volume method, it may be interesting to obtain information about the solution besides of the graphical representation.

Here, we focus on the following data:

• integrals of the solution

• flux through parts of the boundary

Sample problem

Here, we define a sample problem for discussing these issues, which could be formulated in a more general way as well.

Let us define the following reaction - diffusion system in $\mathrm{\Omega }$:

$$\begin{array}{rl}{\partial }_t{u}_1-\mathrm{\nabla }\cdot \mathrm{\nabla }{u}_1+r(u_1,u_2)& =f=1.0\\ {\partial }_t{u}_2-\mathrm{\nabla }\cdot \mathrm{\nabla }{u}_1-r(u_1,u_2)& =0\end{array}$$

with boundary conditions $u_2=0$ on ${\mathrm{\Gamma }}_2\subset \partial \mathrm{\Omega }$ and $r(u_1,u_2)=u_1+0.1u_2$.

The source $f$ creates species $u_1$ which reacts to $u_2$, $u_2$ then leaves the domain at boundary ${\mathrm{\Gamma }}_2$.

... be careful with the sign: reaction is on the left hand side, source on the right hand side.

Create the system, enable species, set boundary condition, solve, create initial value:

Stationary case

For this problem, we have the following flux balances derived from the equations and from Gauss' theorem:

$$\begin{array}{rl}{\int }_{\mathrm{\Omega }}r(u_1,u_2)d\omega & ={\int }_{\mathrm{\Omega }}fd\omega \\ {\int }_{\mathrm{\Omega }}-r(u_1,u_2)d\omega & ={\int }_{{\mathrm{\Gamma }}_2}\mathrm{\nabla }{u}_2\cdot \overrightarrow{n}ds\end{array}$$

The volume integrals can be approximated based on the finite volume subdivision $\mathrm{\Omega }={\cup }_{i\in \mathcal{N}}{\omega }_i$:

$$\begin{array}{rl}{\int }_{\mathrm{\Omega }}r(u_1,u_2)d\omega & \approx \sum _{i\in \mathcal{N}}|{\omega }_i|r(u_{1,i},u_{2,i})\\ {\int }_{\mathrm{\Omega }}fd\omega & \approx \sum _{i\in \mathcal{N}}|{\omega }_i|f_i\end{array}$$

The integrate method of VoronoiFVM provides a possibility to calculate the volume integral of a function of a solution as described above. It returns a num_species$\times$num_regions matrix of the integrals of the function of the unknowns over the different subdomains (here, we have only one):

• Amount of $u_1$ and $u_2$ in the domain aka integral over identity storage function:

• Amount of species created by source term per unit time:

• Amount of reaction per unit time:

Reaction == creation ?

Let us check our first identity: creation rate = reaction rate:

Outflow == reaction ?

The test function trick

This trick goes back to work of H. Gajewski in the field of semiconductor device simulation.

But what about the boundary integral ? Here, we use a trick to cast the surface integral to a volume integral with the help of a test function:

Let $T(x)$ be the solution of the Laplace problem $-{\mathrm{\nabla }}^{2}T=0$ in $\mathrm{\Omega }$ and the boundary conditions

$$\begin{array}{rl}T& =0\text{at}{\mathrm{\Gamma }}_4\\ T& =1\text{at}{\mathrm{\Gamma }}_2\\ {\partial }_{n}T& =0\text{at}{\mathrm{\Gamma }}_1,{\mathrm{\Gamma }}_3\end{array}$$

VoronoiFVM.jl provides a special API for obtaining such a test function:

Write $\overrightarrow{j}=-\mathrm{\nabla }u$. and assume $\mathrm{\nabla }\cdot \overrightarrow{j}+r=f$

$$\begin{array}{rl}{\int }_{{\mathrm{\Gamma }}_2}\overrightarrow{j}\cdot \overrightarrow{n}ds& ={\int }_{{\mathrm{\Gamma }}_2}T\overrightarrow{j}\cdot \overrightarrow{n}ds\text{due to}T=1\text{on}{\mathrm{\Gamma }}_2\\ & ={\int }_{\partial \mathrm{\Omega }}T\overrightarrow{j}\cdot \overrightarrow{n}ds\text{due to}T=0\text{on}{\mathrm{\Gamma }}_4,\overrightarrow{j}\cdot \overrightarrow{n}=0\text{on}{\mathrm{\Gamma }}_1,{\mathrm{\Gamma }}_3\\ & ={\int }_{\mathrm{\Omega }}\mathrm{\nabla }\cdot (T\overrightarrow{j})d\omega \text{(Gauss)}\\ & ={\int }_{\mathrm{\Omega }}\mathrm{\nabla }T\cdot \overrightarrow{j}d\omega +{\int }_{\mathrm{\Omega }}T\mathrm{\nabla }\cdot jd\omega \\ & ={\int }_{\mathrm{\Omega }}\mathrm{\nabla }T\cdot \overrightarrow{j}d\omega +{\int }_{\mathrm{\Omega }}T(f-r)d\omega \end{array}$$

and we approximate

$$\begin{array}{r}{\int }_{\mathrm{\Omega }}\mathrm{\nabla }T\cdot \overrightarrow{j}d\omega \approx \sum _{k,l}\frac{|{\omega }_k\cap {\omega }_l|}{h_{k,l}}g(u_k,u_l)(T_k-T_l)\end{array}$$

where the sum runs over pairs of neigboring control volumes.

The integrate method with a test function parameter returns a value for each species, the sign convention assumes that species leaving the domain lead to negative values.

Check that none of $u_1$ leaves the domain through the boundary:

Check that creation of $u_2$ in the reaction is balanced by $u_2$ leaving the domain through ${\Gamma }_2$:

So we indeed can confirm the requirement for the right balance of source, reaction and outflow.

Transient problem

For the transient case, in addition, we need to consider the time derivative part along with reaction and source. In the derivation of the test function procedure, under the assumption of the implicit Euler time discretization method, this can be achieved by handling the finite difference in time along with source and reaction.

The call to solve corresponds to a new API for transient solutions. It returns a solution object which is compatible to that from DifferentialEquations.jl.

In particular, a TransientSolution tsol can be accessed as follows:

• tsol[it] contains the solution for timestep i

• tsol[ispec,:,it] contains the solution for component ispec at timestep i

• tsol(t) returns a (linearly) interpolated solution value for t.

• tsol.t[it] is the corresponding time

• tsol[ispec,ix,it] refers to solution of component ispec at node ix at moment it

Time:

6.7

From the solution we now can calculate the normal flux via our test function "trick", once again through the API provided by VoronoiFVM:

For increasing time, the outflow rate should approach the value we calculated from the stationary solution:

The overall amount of species which left the domain is can be calculated integrating the discrete outflow rate over time

$$I_{all}={\int }_{t_0}^{t_{end}}{\int }_{{\mathrm{\Gamma }}_2}\mathrm{\nabla }{u}_2\cdot \overrightarrow{n}ds$$

The amount of species created via the source term (measured in F) integrated over time should be equal to the sum of the amount of species left in the domain at the very end of the evolution and the amount of species which left the domain:

$${\int }_{t_0}^{t_{end}}{\int }_{\mathrm{\Omega }}fd\omega dt={\int }_{\mathrm{\Omega }}(u_1+u_2)d\omega +I_{all}$$