modified integrate calculation + added displacement current calculation

[dilaraabdel]

Features

This pull request introduces the capability to:

• Compute the displacement current $I = \int_\Gamma \partial_t \mathbf{j} \cdot \boldsymbol{\nu} dS$, where, for example, $\mathbf{j} = - \varepsilon \nabla \psi$ via the new method integrate_displacement.
• Access both integral contributions that enter the definition of the total current. Following the test-function approach by Gajewski (WIAS Report No 6, 1993), the current can be written as

$$I = \int_\Gamma \mathbf{j} \cdot \boldsymbol{\nu} dS = \int_\Omega \nabla T \cdot \mathbf{j} dx + \int_\Omega T \nabla \cdot \mathbf{j} dx =: I_1 + I_2.$$

The term $I_1 = \int_\Omega \nabla T \cdot \mathbf{j} dx$ can be obtained via integrate_edgebatch, while $I_2 = \int_\Omega T \nabla \cdot \mathbf{j} dx$ via integral_nodebatch.

Motivation for displacement calculation

Consider a bipolar drift-diffusion model

$$\begin{aligned} - \nabla \cdot (\varepsilon \nabla \psi) &= (z_1 n_1 + z_2 n_2 + C),\\\ z_1 \partial_t n_1 + \nabla \cdot \mathbf{j}_1 &= z_1 R,\\\ z_2 \partial_t n_2 + \nabla \cdot \mathbf{j}_2 &= z_2 R. \end{aligned}$$

For semiconductors (van Roosbroeck), we have $z_1 = - z_2 = 1$, while in electrochemistry (Nernst–Planck) $R = 0$.

The total current is defined as ($\mathbf{j}_\psi = - \varepsilon \nabla \psi$)

$$I_{tot} = \int_\Gamma \mathbf{j}_1 \cdot \boldsymbol{\nu} + \mathbf{j}_2 \cdot \boldsymbol{\nu} + \partial_t \mathbf{j}_\psi \cdot \boldsymbol{\nu} dS.$$

Currently, the displacement current is not available, although it plays an important role in transient drift–diffusion-type equations involving the electric potential.
The newest example 161 shows that, in the non-stationary regime, the total current with and without displacement contributions differs.
Therefore, we introduce the new method integrate_displacement to calculate $\int_\Gamma \partial_t \mathbf{j}_\psi \cdot \boldsymbol{\nu} dS$.

Motivation to have access to both integral contributions

It is also useful to have access to both integral contributions $I_1$, $I_2$.
Using the test-function reformulation of the surface integral into a volume integral, we obtain for the above model

$$\begin{aligned} I_{tot} = &\int_\Omega \nabla T \cdot (\mathbf{j}_1 + \mathbf{j}_2 + \partial_t \mathbf{j}_\psi) dx \\\ & + \int_\Omega T (\nabla \cdot \mathbf{j}_1 + \nabla \cdot \mathbf{j}_2 + \partial_t \nabla \cdot \mathbf{j}_\psi) dx. \end{aligned}$$

The second term vanishes in most applications once the PDE model is inserted:

$$\begin{aligned} &\nabla \cdot \mathbf{j}_1 + \nabla \cdot \mathbf{j}_2 + \partial_t \nabla \cdot \mathbf{j}_\psi\\\ &= (z_1 + z_2) R - z_1 \partial_t n_1 - z_2 \partial_t n_2 + \partial_t (z_1 n_1 + z_2 n_2 + C) = 0, \end{aligned}$$

when $R = 0$ or $z_1 = - z_2$.

Instead of the definition $I_\alpha := \int_\Gamma \mathbf{j}_\alpha \cdot \boldsymbol{\nu} dS$, in semiconductor applications, it is often more natural to define the carrier-dependent total currents as

$$\begin{aligned} I_{tot} = \int_\Omega \nabla T \cdot \mathbf{j}_1 dx + \int_\Omega \nabla T \cdot \mathbf{j}_2 dx + \int_\Omega \nabla T \cdot \mathbf{j}_\psi dx =: I_1 + I_2 + I_\psi. \end{aligned}$$

It would be interesting to verify, in the electrochemical context, whether there is also a discrepancy in the definition of the individual carrier-dependent currents.
In any case, to be at least consistent with the semiconductor device community, access to the integrals $\int_\Omega \nabla T \cdot \mathbf{j}_\alpha dx$ is necessary.
This is now provided via integrate_edgebatch.

[j-fu]

Thank you. I am ok with the logic of splitting of the flux integrals.

However, after looking at the paper, I don't think that we don't need integrate_displacement_nodebatch, as IMHO just need to add the integral
1/t ∫ ∇T * (ε∇ϕ - ε∇ϕ_old) to the calculation.

So IMHO, example II is the correct way.

Another point: I am not sure if we should call this integrate_displacement, because this is very application specific. Something like
integrate_flux_time_derivative would IMHO generic and precise.

When I thought about implementing this, I was about implementing something like

"""
    edgeintegrate(system, T, F,Ut)

Integrate edge function (same signature as flux function)
 `F` of  solution vector multiplied by gradient of test function.
The result is an `nspec` vector.
"""

Here one would pass Ut=(U-Uold)/Δt and calculate ∫∇T⋅∇U_t .