Stabilized Mixed-Hybrid Finite Elements of Lowest Order in Rⁿ

1. Introduction

Partial differential models of diffusion-advection-reaction (DAR) type are a relevant member of the family of second order elliptic equations because of their ubiquity in computatonal sciences. For example, systems of nonlinearly coupled DAR equations are frequently used to describe the spatio-temporal dynamics of cellular organisms [1] or charge density of electrons and holes in a semiconductor device [2]. Limiting our focus on the stationary case, DAR models are characterized by two important features:

(F.1) 

conservation;

(F.2) 

monotonicity.

To address (F.1), we introduce an arbitrary bounded volume $\mathcal{V}\subset {\mathbb{R}}^n$ in ${\mathbb{R}}^n$ with boundary $\partial V$ on which an outward unit normal vector ${\mathbf{n}}_{\partial V}$ is defined. We use the symbol u to represent the solution of the DAR model, for example, a concentration, and introduce the vector field $\mathbf{Z}$ and the scalar function $\mathcal{R}$, both depending on u, which represent a mathematical model of the motion of u across $\partial \mathcal{V}$ and the net production of u inside $\mathcal{V}$, respectively.

Definition 1 (Conservation of u) 

The following relation expresses theconservationof u (see [3, Chapter 7])

$$\begin{array}{ccc}& {\int }_{\partial \mathcal{V}}\mathbf{Z}·{\mathbf{n}}_{\partial V}d(\partial \mathcal{V})={\int }_{\mathcal{V}}\mathcal{R}d\mathcal{V}.\hfill & \hfill \end{array}$$

(1a)

If $\mathcal{R}=0$, we say that U is conserved, i.e., the amount of u that enters the volume exactly equals the quantity of u that leaves the volume. If $\mathcal{R}>0$, we say that the amount of u that leaves the volume is larger than the amount of u that enters the volume because of production of u inside the volume. If $\mathcal{R}<0$, we say that the amount of u that leaves the volume is smaller than the amount of u that enters the volume because of consumption of u inside the volume.

Definition 2 (Conservation of $u_h$) 

Let $u_h$ denote the approximate solution of the DAR model computed by a numerical scheme and depending by the discretization parameter h. The scheme is said to beconservativeif it satisfies (1a) at the discrete level, namely

$$\begin{array}{ccc}& {\int }_{\partial \mathcal{V}}{\mathbf{Z}}_h·{\mathbf{n}}_{\partial V}d(\partial \mathcal{V})={\int }_{\mathcal{V}}{\mathcal{R}}_{h}d\mathcal{V},\hfill & \hfill \end{array}$$

(1b)

where ${\mathbf{Z}}_h$ and ${\mathcal{R}}_h$ denote the approximations of $\mathbf{Z}$ and $\mathcal{R}$, respectively, both depending on $u_h$ and h.

To address (F.2), we give the following definition.

Definition 3 (Monotonicity)  

Assume that under suitable conditions on the coefficients of the model the solution u satisfies

$$\begin{array}{ccc}& u\left(\mathbf{x}\right)\ge 0\mathbf{x}\in \overline{\Omega }.\hfill & \hfill \end{array}$$

(2a)

Then, we say that the DAR model ismonotone.

Condition (2a) expresses the physical fact that a concentration of a chemical species or electric charge density cannot become negative. A natural request is that also the numerical scheme for the approximation of the solution of the DAR model is monotone. For this, we give the following definition.

Definition 4 (Monotonicity of the scheme) 

Let us assume that there exists a value $h_0>0$ such that, for all $0<h\le h_0$, we have

$$\begin{array}{ccc}& u_h\left(\mathbf{x}\right)\ge 0\mathbf{x}\in \overline{\Omega }.\hfill & \hfill \end{array}$$

(2b)

If $h_0$ is finite, we say that the scheme isconditionally monotone, otherwise, if $h_0=+\infty$, we say that the scheme isunconditionally monotone(or simply,monotone).

Conditional and unconditional monotonicity are the counterpart of conditional and unconditional absolute stability of a numerical scheme for the Cauchy problem (see [4, Chapter 11]). Conditional monotonicity is due to the fact that the advective and/or reactive term of the DAR model dominates in magnitude over the diffusive term. In such a circumstance, the maximum permitted value $h_0$ of the mesh size may become extremely small with, at least, two negative consequences on the process of numerically solving the discrete DAR equation. The first consequence is that the number $N_h$ of degrees of freedom of $u_h$ is inversely proportional to $h_0$, so that the number of flops (floating-point operations) to solve the linear algebraic system scales unfavorably with $N_h$ (as $\mathcal{O}\left(N_h^3\right)$ if a direct solver is adopted). The second consequence is that the condition number of the matrix of the system grows as $h_0^{-2}$, resulting into a significant deterioration of the solution quality and, if an iterative method is used, of the convergence rate. For these two reasons, ensuring the property of unconditional monotonicity of the numerical scheme is of paramount importance for an affordable, stable and accurate approximation of the solution of a DAR equation.

In this article we illustrate a conservative and monotone scheme for the numerical approximation of the stationary DAR model in a polyhedral domain $\Omega \subset {\mathbb{R}}^n$, $n=1,2,3$, subject to Robin boundary conditions. The scheme is a stabilized mixed-hybrid finite element (MHFE) method of lowest order (see [5]), whose output is a piecewise constant vector-valued function ${\mathbf{j}}_h$ and a piecewise linear continuous scalar-valued function $u_h$, that approximate the flux density $\mathbf{j}$ and concentration u on a nonuniform partition ${\mathcal{T}}_h$ of $\overline{\Omega }$ into simplices K of maximum diameter h.

The first step of the construction of the method follows [6] and consists of a change of dependent variable that transforms, for each $K\in {\mathcal{T}}_h$, the non-symmetric DAR differential operator into a symmetric DR differential operator with an exponentially varying diffusion coefficient. The second step of the construction consists of using a quadrature formula to approximate the local bilinear form into a sum of contributions over the edges of the simplex K. This formula can be applied in all dimensions $n=1$, $n=2$ and $n=3$, and gives rise to a decomposition of ${\mathbf{j}}_h$ into its $n(n+1)/2$ tangential components on each edge of K that coincide with the Scharfetter-Gummel representation of the flux density in terms of the nodal values of $u_h$ on the edge (see [7]). This decomposition of ${\mathbf{j}}_h$ on the edges of K resembles the "edge conditions" used to define the divergence-free upwinding method proposed in [8]. The third step of the construction consists of assembling the local stiffness matrix of the stabilized MHFE method into a global stiffness matrix which can be proved to be an M-matrix under suitable geometrical constraints on the partition ${\mathcal{T}}_h$ (see [9, Lemma 2.1]). This algebraic property implies that the hybrid variable $u_h$ computed by the stabilized MHFE method satisfies (4) for all $h>0$ (unconditional monotonicity).

The remainder of the article is organized as follows. Section 2 describes the differential model. Section 3 illustrates the change of dependent variable that locally transforms the DAR equation into a DR equation with an exponentially varying diffusion coefficient. Section 4 describes the construction of the MHFE method and shows that it is conservative in the sense of Definition 2. Section 5 illustrates the matrix form of the scheme whereas Section 6 describes how to stabilize the MHFE method against advective and/or reactive dominating conditions. Stability and error analysis of the MHFE method are conducted in Section 7 whereas Section 8 describes how to recover the computed flux density by post-processing at the level of the single mesh element. Section 9 illustrates a series of simulation tests in all spatial dimensions $n=1,2,3$ to the purpose of validating the performance of the proposed method versus the nonstabilized formulation. Section 10 contains some summarizing conclusions and suggests possible research developments. Appendix A provides the numerical values of the model parameters in the simulation of the case study described in Section 9.2.

2. The Differential Model

Let $\Omega \subset {\mathbb{R}}^n$ denote an open polyhedral domain whose boundary is the surface $\Gamma$ on which an outward normal unit vector $\mathbf{n}$ is defined. A point P belonging to $\overline{\Omega }$ is identified by its coordinates $\mathbf{x}={[x,y,z]}^T$ with respect to a fixed system of reference. The following partial differential equation (PDE) is a diffusion-advection-reaction mathematical model (DAR):

$$\begin{array}{ccc}& \nabla ·\mathbf{j}(u;\mathbf{v})+cu=\mathcal{P}\mathrm{in}\Omega ,\hfill & \hfill \end{array}$$

(3a)

$$\begin{array}{ccc}& \mathbf{j}(u;\mathbf{v})=\mathbf{v}u-D\nabla u\mathrm{in}\Omega ,\hfill & \hfill \end{array}$$

(3b)

$$\begin{array}{ccc}& {\gamma }_R\mathbf{j}(u;\mathbf{v})·\mathbf{n}={\alpha }_{R}u-{\beta }_R\mathrm{on}\Gamma .\hfill & \hfill \end{array}$$

(3c)

The flux density $\mathbf{j}$ (units: $Um{s}^{-1}$) and concentration u (units: U) are the dependent variables whereas c (units: $s^{-1}$), $\mathcal{P}$ (units: $U{s}^{-1}$), $\mathbf{v}$ (units: $m{s}^{-1}$) and D (units: $m^2{s}^{-1}$) are the model coefficients representing the consumption term, the production term, the advective field and the diffusion coefficient, respectively. The function D is strictly positive in $\Omega$ while c and $\mathcal{P}$ are nonnegative in $\Omega$. The coefficients ${\gamma }_R$, ${\alpha }_R$ (units: $m{s}^{-1}$) and ${\beta }_R$ (units: $Um{s}^{-1}$) are given functions on $\Gamma$. If ${\gamma }_R=0$ we need to have ${\alpha }_R\ne 0$ and the boundary condition (3c) becomes the Dirichlet condition $u={\beta }_R/{\alpha }_R$. If ${\gamma }_R\ne 0$ and ${\alpha }_R\ne 0$, the boundary condition (3c) becomes a Robin boundary condition for u, while if ${\gamma }_R\ne 0$ and ${\alpha }_R=0$, the boundary condition (3c) becomes a Neumann boundary condition for u. The vector-valued function $\mathbf{v}={\left[v_x,v_y,v_z\right]}^T$ is a given velocity field which gives rise to a passive drifting force on u. Molecular diffusion according to Fick’s law is instead associated with the diffusion coefficient D.

If $c=0$ and $\mathbf{v}={\mathbf{0}}^T$, system (3) is a diffusion model for u; if $c\ne 0$ and $\mathbf{v}={\mathbf{0}}^T$, system (3) is a diffusion-reaction model for u; conversely, if $c=0$ and $\mathbf{v}\ne {\mathbf{0}}^T$, system (3) is a diffusion-advection model for u. In what follows we assume that the DAR boundary value problem (BVP) admits a unique solution depending with continuity on the data.

3. Local Symmetrization of the DAR Model

Let ${\mathcal{T}}_h$ denote a regular triangulation of the domain $\overline{\Omega }$ into simplices K characterized by a diameter $h_K$ and a baricenter ${\overline{\mathbf{x}}}_K$, with ${\displaystyle h:=\underset{K\in {\mathcal{T}}_h}{max}{h}_K}$ (see [10, Chapter 3]).

Assumption 1. 

From now on, we assume that $\mathbf{v}$ and D are replaced by their values ${\mathbf{v}}_K=\mathbf{v}\left({\overline{\mathbf{x}}}_K\right)$ and $D_K=D\left({\overline{\mathbf{x}}}_K\right)$ for each $K\in {\mathcal{T}}_h$.

Using Assumption 1 we can write

$$\begin{array}{ccc}& {\mathbf{j}}_K=D_{K}u\left(D_K^{-1}{\mathbf{v}}_K-{\displaystyle \frac{\nabla u}{u}}\right)=D_{K}u\nabla \left[{\Psi }_K-ln\left({\displaystyle \frac{u}{\overline{u}}}\right)\right]\forall K\in {\mathcal{T}}_h,\hfill & \hfill \end{array}$$

(4)

where $\overline{u}$ is a positive arbitrary reference concentration (units: U) and ${\Psi }_K$ is defined as

$$\begin{array}{ccc}& {\Psi }_K\left(\mathbf{x}\right)={\displaystyle \frac{{\mathbf{v}}_K·\mathbf{x}}{D_K}}\forall K\in {\mathcal{T}}_h.\hfill & \hfill \end{array}$$

(5)

Let us introduce the function

$$\begin{array}{ccc}& {\phi }_u:={\Psi }_K-ln\left({\displaystyle \frac{u}{\overline{u}}}\right)\forall K\in {\mathcal{T}}_h.\hfill & \hfill \end{array}$$

(6)

Inverting (7), we have

$$\begin{array}{ccc}& u=\overline{u}exp\left({\Psi }_K-{\phi }_u\right)\forall K\in {\mathcal{T}}_h,\hfill & \hfill \end{array}$$

(7)

which can be written as

$$\begin{array}{ccc}& u=\rho exp\left({\Psi }_K\right)\forall K\in {\mathcal{T}}_h,\hfill & \hfill \end{array}$$

(8)

upon introducing the new dependent variable $\rho$. Replacing (8) into the definition of ${\mathbf{j}}_K$ we obtain

$$\begin{array}{ccc}& {\mathbf{j}}_K={\mathbf{v}}_K\rho exp\left({\Psi }_K\right)-D_{K}exp\left({\Psi }_K\right)\nabla \rho -D_K\nabla {\Psi }_K\rho exp\left({\Psi }_K\right)=-D_{K}exp\left({\Psi }_K\right)\nabla \rho \forall K\in {\mathcal{T}}_h,\hfill & \hfill \end{array}$$

(9)

having used (5) to compute $\nabla {\Psi }_K$. The new dependent variable $\rho$ is known as "Slotboom variable" because it was originally introduced in [11] by J. Slotboom and subsequently adopted in [6] to construct mixed finite element approximations of the continuity equations for electrons and holes in the drift-diffusion model for semiconductor devices, We see that the "Slotboom change of variable" (8) has the effect of transforming the DAR differential operator into a symmetric DR differential operator characterized by a modified diffusion coefficient whose range depends exponentially on the dimensionless potential ${\Psi }_K$. Inverting (9) to express $\rho$ as a function of u, the flux density can be written as

$$\begin{array}{ccc}& {\mathbf{j}}_K=-D_{K}exp\left({\Psi }_K\right)\nabla \left[uexp(-{\Psi }_K)\right]\forall K\in {\mathcal{T}}_h.\hfill & \hfill \end{array}$$

(10)

4. The MHFE Discretization of Lowest Order

To simplify the presentation of the method, we assume that ${\gamma }_R\left(\mathbf{x}\right)={\beta }_R\left(\mathbf{x}\right)=0$ for all $\mathbf{x}\in \Gamma$, which means that the DAR model is supplied with homogeneous Dirichlet boundary conditions. We define the set ${\mathcal{N}}_h$ of the indices of all vertices of ${\mathcal{T}}_h$ and we set $N_h=\mathrm{card}\left({\mathcal{N}}_h\right)$. The set ${\mathcal{N}}_h$ is divided into two subsets, namely, the set ${\mathcal{N}}_{h,int}$ of the vertices belonging to the interior of $\Omega$, and the set ${\mathcal{N}}_{h,\Gamma }$ of the vertices belonging to the boundary $\Gamma$. We set $N_{h,int}=\mathrm{card}\left({\mathcal{N}}_{h,int}\right)$ and $N_{h,\Gamma }=\mathrm{card}\left({\mathcal{N}}_{h,\Gamma }\right)$, with $N_h=N_{h,int}+N_{h,\Gamma }$.

Let us introduce the functional spaces where to seek the flux density $\mathbf{j}$ and concentration u, respectively:

$$\begin{array}{ccc}& Q={\left[L^2(\Omega )\right]}^n,\hfill & \hfill \end{array}$$

(11a)

$$\begin{array}{ccc}& V=H_0^1(\Omega ).\hfill & \hfill \end{array}$$

(11b)

Correspondingly, the finite dimensional spaces where to seek the approximations ${\mathbf{j}}_h$ and $u_h$ computed by the MHFE method of lowest order are defined as:

$$\begin{array}{ccc}& Q_h=\left\{\mathbf{q}\in Q\mathrm{s}.\mathrm{t}.\mathbf{q}\in {\left[{\mathbb{P}}_0\left(K\right)\right]}^n\forall K\in {\mathcal{T}}_h\right\},\hfill & \hfill \end{array}$$

(11c)

$$\begin{array}{ccc}& V_h=\left\{v\in V\cap C^0\left(\overline{\Omega }\right)\mathrm{s}.\mathrm{t}.v\in {\mathbb{P}}_1\left(K\right)\forall K\in {\mathcal{T}}_h\right\}.\hfill & \hfill \end{array}$$

(11d)

The MHFE approximation of lowest order of the boundary value problem (3) reads:

Find ${\mathbf{j}}_h\in Q_h$ and $u_h\in V_h$ such that:

$$\begin{array}{ccc}& {\int }_{\Omega }\left[{\left(D_K{e}^{{\Psi }_K}\right)}^{-1}{\mathbf{j}}_h+\nabla \left(u_h{e}^{-{\Psi }_K}\right)\right]·{\tau }_{h}d\Omega =0\forall {\tau }_h\in Q_h,\hfill & \hfill \end{array}$$

(12a)

$$\begin{array}{ccc}& -{\int }_{\Omega }{\mathbf{j}}_h·\nabla v_{h}d\Omega +{\int }_{\Omega }c{u}_h{v}_{h}d\Omega ={\int }_{\Omega }\mathcal{P}{v}_{h}d\Omega \forall v_h\in V_h.\hfill & \hfill \end{array}$$

(12b)

Let $V_i\in {\mathcal{N}}_{h,int}$ denote a vertex of ${\mathcal{T}}_h$ belonging to the interior of $\Omega$, $i=1,\dots ,N_{h,int}$. We associate with $V_i$ the set ${\sigma }_i\subset \Omega$ constituted by the union of the simplexes $K\in {\mathcal{T}}_h$ sharing $V_i$ as a common vertex (see Figure 1 for a graphical representation of ${\sigma }_i$ in the case $n=2$). We indicate by ${\phi }_i:\overline{\Omega }\to \mathbb{R}$ the basis function spanning the vector space $V_h$ that is linear on each simplex $K\in {\sigma }_i$, and is such that ${\phi }_i\left(V_i\right)=1$ and ${\phi }_i\left(V_j\right)=0$ for all $j=1,\dots ,N_h$, $j\ne i$. Based on this definition, the set ${\sigma }_i$ is the support of ${\phi }_i$ (also referred to as the "patch" associated with $V_i$). The following result shows that the MHFE method is conservative in the sense of Definition 2.

Theorem 1 

(Conservation of the MHFE method). The MHFE method satisfies (1b) by setting $\mathcal{V}={\sigma }_i$, ${\mathbf{Z}}_h={\phi }_i{\mathbf{j}}_h$ and ${\mathcal{R}}_h={\phi }_i\left(\mathcal{P}-c{u}_h\right)$, respectively.

Proof. 

The result immediately follows by taking $v_h={\phi }_i$ in (12b) and using integration by parts in the first integral on the left hand side.    □

Remark 1 

(Conservation in 2D). In the case $n=2$, the conservation property expressed by Theorem 1 corresponds to stating that the weighted sum of the jump of the normal flux density across the edges belonging to ${\sigma }_i$ balances the net production over the support. This observation generalizes the considerations developed in [8] to support the stability of a novel divergence-free upwinding finite element approximation of the DAR equation in the case where $c=\mathcal{P}=0$.

Using the fact that ${\mathbf{j}}_h$ is constant for each element K, denoted henceforth ${\mathbf{j}}_{K,h}$, and that ${\tau }_h\ne \mathbf{0}$ over each element K and is equal to zero elsewhere, equation (12a) can be solved independently for each mesh element to obtain

$$\begin{array}{ccc}& {\displaystyle {\mathbf{j}}_{K,h}·{\tau }_h{\int }_K{e}^{-{\Psi }_K}dK=-D_K{\tau }_h·{\int }_K\nabla \left(u_h{e}^{-{\Psi }_K}\right)dK{\tau }_h\in {\left[{\mathbb{P}}_0\left(K\right)\right]}^n,K\in {\mathcal{T}}_h,}\hfill & \hfill \end{array}$$

from which we get

$$\begin{array}{ccc}& {\mathbf{j}}_{K,h}=-D_K{\mathcal{H}}_K\left(e^{{\Psi }_K}\right){\displaystyle \frac{{\int }_K\nabla \left(u_h{e}^{-{\Psi }_K}\right)dK}{\left|K\right|}}K\in {\mathcal{T}}_h.\hfill & \hfill \end{array}$$

(13)

The quantity

$$\begin{array}{ccc}& {\displaystyle {\mathcal{H}}_K\left(e^{{\Psi }_K}\right):={\left[\frac{{\int }_K{e}^{-{\Psi }_K}dK}{\left|K\right|}\right]}^{-1}K\in {\mathcal{T}}_h}\hfill & \hfill \end{array}$$

(14)

is the harmonic average of $e^{{\Psi }_K}$ over the n-dimensional simplex $K\in {\mathcal{T}}_h$, having indicated by $\left|K\right|$ the n-dimensional volume of K (see [12] for the analysis of the role of the harmonic average in the numerical approximation of second-order ellliptic problems). Replacing the expression (13) in (12b), we obtain the final form of the MHFE formulation of lowest order in terms of the sole hybrid variable $u_h$:

Find $u_h\in V_h$ such that

$$\begin{array}{ccc}& B(u_h,v_h)=F\left(v_h\right)\forall v_h\in V_h,\hfill & \hfill \end{array}$$

(15a)

where:

$$\begin{array}{ccc}& B({\mathcal{U}}_h,{\mathcal{V}}_h)=\sum _{K\in {\mathcal{T}}_h}{D}_K{\mathcal{H}}_K\left(e^{{\Psi }_K}\right){\int }_K\nabla \left({\mathcal{U}}_h{e}^{-{\Psi }_K}\right)·\nabla {\mathcal{V}}_{h}dK+{\int }_{\Omega }c{\mathcal{U}}_h{\mathcal{V}}_{h}d\Omega {\mathcal{U}}_h,{\mathcal{V}}_h\in V_h,\hfill & \hfill \end{array}$$

(15b)

$$\begin{array}{ccc}& F\left({\mathcal{V}}_h\right)={\int }_{\Omega }\mathcal{P}{\mathcal{V}}_{h}d\Omega {\mathcal{V}}_h\in V_h.\hfill & \hfill \end{array}$$

(15c)

Assertion 1 

(Well-posedness). We assume henceforth that the discrete problem (15a) is well posed, i.e., it admits a unique solution $u_h\in V_h$ which depends with continuity on the model data.

5. Matrix Form of the MHFE Method

As it is customary in the context of the finite element methodology (see [13]), we introduce the local bilinear form composed of a diffusion-advection contribution and a reaction contribution

$$\begin{array}{ccccc}& B^K({\phi }_j,{\phi }_i)=B^{K,da}({\phi }_j,{\phi }_i)+B^{K,r}({\phi }_j,{\phi }_i)\hfill & & i,j=1,\dots ,n+1,\hfill & \hfill \end{array}$$

(16)

and the local linear functional $F^K\left({\phi }_i\right)$, $i=1,\dots ,(n+1)$, defined as:

$$\begin{array}{ccccc}& B^{K,da}({\phi }_j,{\phi }_i)=D_K{\mathcal{H}}_K\left(e^{{\Psi }_K}\right){\int }_K\nabla \left({\phi }_j{e}^{-{\Psi }_K}\right)·\nabla {\phi }_{i}dK\hfill & & i,j=1,\dots ,n+1,\hfill & \hfill \end{array}$$

(17a)

$$\begin{array}{ccccc}& B^{K,r}({\phi }_j,{\phi }_i)={\int }_{K}c{\phi }_j{\phi }_{i}dK\hfill & & i,j=1,\dots ,n+1,\hfill & \hfill \end{array}$$

(17b)

$$\begin{array}{ccccc}& F^K\left({\phi }_i\right)={\int }_K\mathcal{P}{\phi }_{i}dK\hfill & & i=1,\dots ,n+1.\hfill & \hfill \end{array}$$

(17c)

In all the above integrals, we denote by ${\phi }_q$, $q=1,\dots ,n+1$, the linear function such that ${\phi }_q\left(V_q\right)=1$ and ${\phi }_q\left(V_p\right)=0$ for $p\ne q$, with $p,q=1,\dots ,(n+1)$. Associated with the bilinear and linear forms in (17) we introduce the local diffusion-advection and reaction matrices of size $(n+1)\times (n+1)$:

$$\begin{array}{ccccc}& {\left[{\mathbb{B}}^{K,da}\right]}_{ij}=B^{K,da}({\phi }_j,{\phi }_i)\hfill & & i,j=1,\dots ,n+1,\hfill & \hfill \end{array}$$

(18a)

$$\begin{array}{ccccc}& {\left[{\mathbb{B}}^{K,r}\right]}_{ij}=B^{K,r}({\phi }_j,{\phi }_i)\hfill & & i,j=1,\dots ,n+1,\hfill & \hfill \end{array}$$

(18b)

and the local load vectors of size $(n+1)\times 1$ (column vectors):

$$\begin{array}{ccccc}& {\left[{\mathbb{F}}^K\right]}_i=F^K\left({\phi }_i\right)\hfill & & i=1,\dots ,n+1.\hfill & \hfill \end{array}$$

(19a)

After looping over the mesh elements to assemble the local stiffness matrix and load vector into the corresponding global matrix and right-hand side, and enforcing the boundary conditions, we obtain the following linear system

$$\begin{array}{ccc}& \mathbb{B}\mathbf{U}=\mathbb{F},\hfill & \hfill \end{array}$$

(20)

which represents the algebraic counterpart of Equation (15a).

Assertion 2 (Well-posedness of (20) Analogous to Assertion 1, we assume henceforth that the linear system (20) admits a unique solution $\mathbf{U}\in {\mathbb{R}}^{N_h\times 1}$.

6. Stabilization Via Numerical Quadrature

To implement in a computer code the bilinear forms and the linear functionals defined in (18) and (19), we make use of specific approximations and quadrature formulas. Let us start with the bilinear form (17a) that accounts for the diffusion and advection terms in the DAR model, and replace the function $u_h{e}^{-{\Psi }_K}$ in (13) with its ${\mathbb{P}}_1$-interpolant ${\Pi }_1^K\left(u_h{e}^{-{\Psi }_K}\right)$ over the element K. By doing this, the local approximate flux density becomes

$$\begin{array}{ccc}& {\mathbf{j}}_{K,h}=-D_K{\mathcal{H}}_K\left(e^{{\Psi }_K}\right)\nabla {\Pi }_1\left(u_h{e}^{-{\Psi }_K}\right)K\in {\mathcal{T}}_h.\hfill & \hfill \end{array}$$

(21)

Let us examine the expression (21) in the one-dimensional (1D) case ($n=1)$. Set $K=(x_1,x_2)$ with $\left|K\right|=x_2-x_1:=\ell$ and indicate by ${\Psi }_1$, ${\Psi }_2$ the nodal values of the dimensionless potential, with $\Delta \Psi ={\Psi }_2-{\Psi }_1$, so that ${\Psi }_K\left(s\right)={\Psi }_1+\Delta \Psi {\displaystyle \frac{s}{\ell }}$, with $s\in [0,\ell ]$. The harmonic average of the function $e^{{\Psi }_K}$ is

$$\begin{array}{ccc}& {\mathcal{H}}_K\left(e^{{\Psi }_K}\right)={\left[{\displaystyle \frac{{\int }_0^{\ell }exp(-{\Psi }_1-\Delta \Psi \frac{s}{\ell })ds}{\ell }}\right]}^{-1}={\left[e^{-{\Psi }_1}\left({\displaystyle \frac{exp(-\Delta \Psi )-1}{-\Delta \Psi }}\right)\right]}^{-1}=e^{{\Psi }_1}\mathtt{Be}(-\Delta \Psi ),\hfill & \hfill \end{array}$$

(22)

where $\mathtt{Be}\left(z\right):=z/(e^z-1)$ for $z\ne 0$, such that $\mathtt{Be}\left(0\right)=1$, represents the inverse of the Bernoulli function of $z\in \mathbb{R}$. Inserting (22) in (21), in the case $n=1$, we obtain

$$\begin{array}{ccc}& j_{K,h}=-D_K{e}^{{\Psi }_1}\mathtt{Be}(-\Delta \Psi )\left[{\displaystyle \frac{u_2{e}^{-{\Psi }_2}-u_1{e}^{-{\Psi }_1}}{\ell }}\right]\hfill & \hfill \\ & =-{\displaystyle \frac{D_K}{\ell }}\left[e^{-\Delta \Psi }\mathtt{Be}(-\Delta \Psi )u_2-u_1\mathtt{Be}(-\Delta \Psi )\right]={\displaystyle \frac{D_K}{\ell }}\left[u_1\mathtt{Be}(-\Delta \Psi )-\mathtt{Be}\left(\Delta \Psi \right)u_2\right],\hfill & \hfill \end{array}$$

(23)

having used the property $e^{-z}\mathtt{Be}(-z)=\mathtt{Be}\left(z\right)$. The expression (23) coincides with the Scharfetter-Gummel (SG) formula originally derived in [7] and subsequently analyzed in [14,15]. Because of its second-order accuracy and numerical stability, and the automatic "upwind" weighting of the nodal values $u_1$ and $u_2$, the one-dimensional SG formula (23) is universally employed in the simulation packages of the drift-diffusion model of semiconductor devices in multiple spatial dimensons (see [16,17]).

Motivated by the connection between the MHFE formulation and the SG method in the 1D case, we now address how to use the SG formula to compute the entries of the local matrix ${\mathbb{B}}^{K,da}$ in the case $n>1$. With this aim, we recall the following result proved in [18] and [9].

Lemma 1. 

Let $\mathbf{Z}\in {\left[{\mathbb{P}}_0\left(K\right)\right]}^n$ be a given constant vector in $K\in {\mathcal{T}}_h$. Let also $e\in \partial K$ denote the edge with vertices $V_p$ and $V_q$, length ${\ell }_e$, and define on e a tangent vector ${\mathbf{t}}_e=({\mathbf{x}}_q-{\mathbf{x}}_p)/{\ell }_e$ directed from $V_p$ to $V_q$. Then, we have

$$\begin{array}{ccc}& -{\int }_K\mathbf{Z}·\nabla {\phi }_{i}dK=\sum _{e\in \partial K}{\omega }_e^K{\ell }_e{Z}_e{\delta }_e{\phi }_{i}i=1,\dots ,n+1,\hfill & \hfill \end{array}$$

(24)

having set ${\omega }_e^K:=-{\int }_K\nabla {\phi }_q·\nabla {\phi }_{p}dK$ and $Z_e:=\mathbf{Z}·{\mathbf{t}}_e$, and where ${\delta }_{e}f:=f\left(V_p\right)-f\left(V_q\right)$ for any $f\in C^0\left(\overline{e}\right)$.

For any vertex $V_i$ of K, we denote by $\mathtt{A}\left(\mathtt{i}\right)$ the set of the vertices of K that are adjacent to $V_i$, $i=1,\dots ,n+1$. For example, taking $i=1$, in the case $n=1$ we have $\mathtt{A}\left(\mathtt{i}\right)=\left\{V_2\right\}$, in the case $n=2$ we have $\mathtt{A}\left(\mathtt{i}\right)=\{V_2,V_3\}$, while if $n=3$ we have $\mathtt{A}\left(\mathtt{i}\right)=\{V_2,V_3,V_4\}$. The other expressions of the adjaciencies for $i=2,\dots ,n+1$ can be determined cyclically. Using the definition of ${\omega }_e^K$ and $\mathtt{A}\left(\mathtt{i}\right)$, relation (24) becomes

$$\begin{array}{ccc}& -{\int }_K\mathbf{Z}·\nabla {\phi }_{i}dK=\sum _{j\in \mathtt{A}\left(i\right)}{\gamma }_{ij}^K{Z}_{ij}i=1,\dots ,n+1,\hfill & \hfill \end{array}$$

(25)

where

$$\begin{array}{ccc}& {\gamma }_{ij}^K:=\left(-{\int }_K\nabla {\phi }_j·\nabla {\phi }_{i}dK\right)\left|e_{ij}\right|,\hfill & \hfill \end{array}$$

(26)

and $Z_{ij}$ is the tangential component of $\mathbf{Z}$ along the edge $e_{ij}$ directed from vertex $V_i$ to vertex $V_j$, with length $|e_{ij}|$ and tangent unit vector ${\mathbf{t}}_{ij}$.

Assertion 3. 

Let ${\mathbb{B}}^{K,da}$ and ${\mathbf{U}}_K$ denote the local $(n+1)\times (n+1)$ stiffness matrix accounting for the diffusive and advective contributions and the column vector of size $n+1$ containing the nodal values of $u_h$ at the vertices of K, respectively. Let us also set $\mathbf{Z}:={\mathbf{j}}_{K,h}$. We approximate the matrix-vector product representation of (25) by means of the following quadrature formula

$$\begin{array}{ccc}& {\mathbb{B}}_h^{K,da}(i,:){\mathbf{U}}_K=D_K\sum _{j\in \mathtt{A}\left(i\right)}{\displaystyle \frac{{\gamma }_{ij}^K}{|e_{ij}|}}\left[\mathtt{bnij}{u}_i-\mathtt{bpij}{u}_j\right]i=1,\dots ,n+1,j=1,\dots ,n+1,j\ne i,\hfill & \hfill \end{array}$$

(27a)

where:

$$\begin{array}{ccc}& \mathtt{bnij}=\left\{\begin{array}{cc}\mathtt{Be}(-\Delta \Psi )\hfill & \mathrm{stabilized}\mathrm{case}\hfill \\ {\displaystyle 1+\frac{\Delta \Psi }{2}}\hfill & \mathrm{non}-\mathrm{stabilized}\mathrm{case}\hfill \end{array}\right.\hfill & \hfill \end{array}$$

(27b)

$$\begin{array}{ccc}& \mathtt{bpij}=\left\{\begin{array}{cc}\mathtt{Be}(+\Delta \Psi )\hfill & \mathrm{stabilized}\mathrm{case}\hfill \\ {\displaystyle 1-\frac{\Delta \Psi }{2}}\hfill & \mathrm{non}-\mathrm{stabilized}\mathrm{case}\hfill \end{array}\right.\hfill & \hfill \end{array}$$

(27c)

with $\Delta \Psi :=\Psi \left(V_j\right)-\Psi \left(V_i\right)$, $i,j=1,\dots ,n+1$.

Remark 2. 

The quadrature formula (27a) consists of replacing ${\mathbf{j}}_{K,h}·{\mathbf{t}}_{ij}$ with the 1D SG expression (23). This approach is the analog of the "edge conditions" used to define the divergence-free upwinding method proposed and investigated in [8] and [19].

We complete the characterization of the MHFE method of lowest order by applying the n-dimensional trapezoidal quadrature to compute the entries of the local reaction matrix and load vector, obtaining:

$$\begin{array}{ccc}& {\left[{\mathbb{B}}^{K,r}\right]}_{ij,h}=c\left(V_i\right){\displaystyle \frac{\left|K\right|}{n+1}}{\delta }_{ij}i,j=1,\dots ,n+1K\in {\mathcal{T}}_h,\hfill & \hfill \end{array}$$

(28)

$$\begin{array}{ccc}& {\left[{\mathbb{F}}^K\right]}_{i,h}=\mathcal{P}\left(V_i\right){\displaystyle \frac{\left|K\right|}{n+1}}i=1,\dots ,n+1K\in {\mathcal{T}}_h.\hfill & \hfill \end{array}$$

(29)

The stabilized MHFE method of lowest order for the discretization of the DAR problem (3) reads:

Find $u_h\in V_h$ such that

$$\begin{array}{ccc}& B_h(u_h,v_h)=F_h\left(v_h\right)\forall v_h\in V_h,\hfill & \hfill \end{array}$$

(30a)

where:

$$\begin{array}{ccc}& B_h(u_h,v_h)=B_h^{da}(u_h,v_h)+B_h^r(u_h,v_h),\hfill & \hfill \end{array}$$

(30b)

$$\begin{array}{ccc}& F_h\left(v_h\right)=F_h^{int}\left(v_h\right).\hfill & \hfill \end{array}$$

(30c)

The linear algebraic system of size $N_h$ associated with (30a) is

$$\begin{array}{ccc}& {\mathbb{B}}^{stab}\mathbf{U}={\mathbb{F}}^{stab},\hfill & \hfill \end{array}$$

(31)

where ${\mathbb{B}}^{stab}$ is the global stiffness matrix of size $N_h\times N_h$ whereas $\mathbf{U}$ and ${\mathbb{F}}^{stab}$ are the solution vector and global load vector, both of size $N_h\times 1$, respectively.

Remark 3 

(Equivalence with the EAFE method). The entries of the stiffness matrix and load vector in the linear system (31) coincide with those of the edge-averaged finite element method (EAFE) proposed in [9].

7. Stability and Error Analysis

In this section we leverage the theoretical analysis conducted in previous works to characterize the stability and convergence of the MHFEM illustrated in Section 6.

Theorem 2 

(Stability). Let us assume that c and $\mathcal{P}$ are given nonnegative functions in Ω. Let us also assume that the following geometrical constraint is satisfied for each edge $E\in {\mathcal{T}}_h$

$$\begin{array}{ccc}& {\omega }_E=\sum _{K\cap E\ne \varnothing }{\omega }_e^K\ge 0.\hfill & \hfill \end{array}$$

(32)

Then, the global stiffness matrix ${\mathbb{B}}^{stab}$ is an M-matrix, ${\mathbf{U}}_i>0$, $i=1,\dots ,N_h$ and the stabilized MHFE method of lowest order is monotone.

Proof. 

See the proof of Lemma 3.2 in [9]    □

Remark 4. 

In the case $n=2$, the geometrical constraint (32) is equivalent to requiring that ${\mathcal{T}}_h$ is a Delaunay triangulation (see [20] for the definition).

Theorem 3 

(Convergence). Assume that the solution pair $(u,\mathbf{j})$ of the DAR problem (3) is such that $u\in H^2(\Omega )$ and $\mathbf{j}\in {\left(H^1(\Omega )\right)}^n$, respectively. Then, there exists a positive constant C independent of h such that the following optimal error estimate holds

$$\begin{array}{ccc}& \parallel u-u_h{\parallel }_V+{\parallel \mathbf{j}-{\mathbf{j}}_h\parallel }_Q\le Ch\left({\parallel u\parallel }_{H^2(\Omega )}+{\parallel \mathbf{j}\parallel }_{{\left(H^1(\Omega )\right)}^n}\right),\hfill & \hfill \end{array}$$

(33)

with ${\parallel ·\parallel }_V$ and ${\parallel ·\parallel }_Q$ denoting the norm in the spaces V and Q defined in (11), respectively.

Proof. 

The proof combines the error analysis on the bilinear form $B_h$ conducted in [9, Section 6], the abstract theory for mixed hybrid methods with numerical quadrature developed in [21, Chapter III, Section 11] and Theorem 12.2 of [21, Chapter IV, Section 12].    □

8. Post-Processing to Evaluate ${\mathbf{J}}_h$

Once system (31) is solved, we need to post-process the available nodal values of $u_h$ to determine the approximate flux density ${\mathbf{j}}_h$. A possible approach is suggested in [22], where the flux density ${\mathbf{j}}_h{|}_K$ inside the element $K\in {\mathcal{T}}_h$ is expressed as a function of the tangential flux density along each edge $e\in \partial K$ connecting vertices $V_i$ and $V_j$, $i,j=1,\dots ,n+1$, with orientation from $V_i$ to $V_j$, using the following interpolation formula

$$\begin{array}{ccc}& {\displaystyle \mathbf{j}{|}_K=\sum _{e\in \partial K}{\mathcal{C}}_e\left(\mathbf{j}\right){\mathbf{W}}_e\left({\overline{\mathbf{x}}}_K\right)K\in {\mathcal{T}}_h,}\hfill & \hfill \end{array}$$

(34a)

where

$$\begin{array}{ccc}& {\displaystyle {\mathcal{C}}_e\left(\mathbf{j}\right):={\int }_e\mathbf{j}·{\mathbf{t}}_{e}ds\simeq D_K\left[\mathtt{bnij}{u}_i-\mathtt{bpij}{u}_j\right]}\hfill & \hfill \end{array}$$

(34b)

is the circulation of $\mathbf{j}$ along the edge $e\in \partial K$, with bnij, bp12, $u_i$ and $u_j$ being defined as in (27). In the reconstruction formula (34a), ${\mathbf{t}}_e$ denotes the unit tangent vector along the edge $e\in \partial K$ directed from vertex $V_i$ to vertex $V_j$, and the vector-valued functions ${\mathbf{W}}_e$, also known as Whitney elements (or "edge elements", see [23]), are defined as follows

$$\begin{array}{ccc}& {\mathbf{W}}_e\left(\mathbf{x}\right)={\phi }_i\left(\mathbf{x}\right)\nabla {\phi }_j\left(\mathbf{x}\right)-{\phi }_j\left(\mathbf{x}\right)\nabla {\phi }_i\left(\mathbf{x}\right)\mathbf{x}\in \overline{K}.\hfill & \hfill \end{array}$$

(34c)

9. Numerical Experiments in ${\mathbb{R}}^N$

In this section we conduct a series of numerical tests to validate the computational performance of the MHFE method in its nonstabilized and stabilized versions, by solving the DAR model problem in 1D, 2D and 3D geometries. In all considered cases, model coefficients are selected in such a way that the solution u exhibits the presence of steep internal and/or boundary layers whose accurate resolution would require, at least locally, the adoption of an extremely small mesh size. The main scope of the simulations is to demonstrate the ability of the stabilized MHFE to compute a solution pair $(u_h,{\mathbf{j}}_h)$ with remarkable accuracy even on a coarse-sized mesh without being affected by spurious unphysical oscillations as it happens with the nonstabilized version of the scheme. The post-processing formula (34) is implemented in the codings to solve the DAR model in all spatial dimensions $n=1,2,3$. All computations are performed using the Matlab scientific environment [24] on a HP laptop mounting an Intel Core i7-10510U processor at 1.80GHz.

9.1. The Case $n=1$

In this section we study the DAR model problem (3) in the domain $\Omega =(-1,1)$. The values of the parameters are: $D={10}^{-2}$, $\mathbf{V}={[1,0,0]}^T$, $\mathcal{P}=1$ and $c=5$. A rather coarse uniform partition ${\mathcal{T}}_h$ of 20 elements of size $h=2/20=0.1$ is used in the computations. The boundary conditions are of nonhomogeneous Dirichlet type, with values $u(-1)=0$ and $u(+1)=1$. Figure 2 illustrates the solution $u_h$ and the associated flux density $j_h$ computed by the stabilized MHFE method (solid line, blue color) and non-stabilized MHFE method (solid line, black color). For comparison, the exact solution is superposed (solid line, magenta color). Results indicate that the stabilized formulation captures the essential behavior of u and j even on the coarse mesh utilized in the test. In particular, no spurious oscillations arise in the vicinity of $x=1$, where a steep boundary layer appears in the graphical plot of u and, less markedly, of j. The same comment cannot be applied to the concentration $u_h$ computed by the non-stabilized formulation, which is affected by wiggles that violate the monotonicity of u, reflecting also on the reconstruction of the flux density which shows some over/undershoot with respect to the almost constant profile of the exact solution j.

9.2. The Case $n=2$

In this section we study the stationary continuity equation for the hole concentration p (units: $m^{-3}$) in a two-dimensional (2D) cross-section of a p-n diode according to the so-called "Drift-Diffusion" model of charge transport in semiconductor devices (see [2] and [25]). The mathematical model of the problem consists of the following stationary boundary value problem in conservation form:

$$\begin{array}{ccc}& \nabla ·{\mathbf{j}}_p+cp=\mathcal{P}\mathrm{in}\Omega ,\hfill & \hfill \end{array}$$

(35a)

$$\begin{array}{ccc}& {\mathbf{j}}_p=-q{\mu }_{p}p\nabla \psi -q{D}_p\nabla p\mathrm{in}\Omega ,\hfill & \hfill \end{array}$$

(35b)

$$\begin{array}{ccc}& p=p_D\mathrm{on}{\Gamma }_1,\hfill & \hfill \end{array}$$

(35c)

$$\begin{array}{ccc}& {\mathbf{j}}_p·\mathbf{n}=0\mathrm{on}{\Gamma }_2,\hfill & \hfill \end{array}$$

(35d)

where $\Omega$ is the computational domain shown in Figure 3, whose boundary $\Gamma$ is split into two disjoint subsets, ${\Gamma }_1$ and ${\Gamma }_2$, on which the nonhomogeneous Dirichlet boundary condition (35c) and the homogeneous Neumann boundary condition (35d) are enforced, respectively.

The quantity $q=1.602·{10}^{-19}$ (units: C) is the electron charge, p is the concentration of holes (units: $m^{-3}$), $D_p$ is the hole diffusion coefficient (units: $m^2{s}^{-1}$), ${\mu }_p$ is the hole electric mobility (units: $m^2{V}^{-1}{s}^{-1}$), related to $D_p$ by Einstein’s relation

$$\begin{array}{ccc}& D_p={\mu }_p{V}_{th},\hfill & \hfill \end{array}$$

(36)

with $V_{th}$ denoting the thermal voltage (units: V), and $\psi$ is the electric potential (units: V). The quantities $\tau$ and $p_{eq}$ represent the hole lifetime (units: s) and hole equilibrium concentration (units: $m^{-3}$), respectively, whereas the vector ${\mathbf{j}}_p$ is the hole current density (units: $A{m}^{-2}$). The computational domain $\Omega$ consists of a square whose side is equal to $L=100·{10}^{-6}m$, and represents the 2D cross-section of a p-n junction diode which is the union of three subdomains: the neutral region ${\Omega }_1=(R,\theta )\in (0,R_1)\times (0,\pi /2)$, the annular junction region $\mathrm{JCT}=(R,\theta )\in (R_1,R_2)\times (0,\pi /2)$, and the neutral region ${\Omega }_2=\Omega \setminus ({\Omega }_1\cup \mathrm{JCT})$, having indicated by R and $\theta$ the radial coordinates such that $x=Rcos\left(\theta \right)$ and $y=Rcos\left(\theta \right)$. The neutral regions are characterized by the fact that $\mathcal{P}(R,\theta )=c(R,\theta )=0$ in ${\Omega }_1\cup {\Omega }_2$ while $\psi (R,\theta )={\overline{\psi }}_1$ in ${\Omega }_1$ and $\psi (R,\theta )={\overline{\psi }}_2$ in ${\Omega }_2$, where:

$$\begin{array}{ccc}& {\overline{\psi }}_1=V_1+{\varphi }_{bi,1},\hfill & \hfill \end{array}$$

(37a)

$$\begin{array}{ccc}& {\overline{\psi }}_2=V_2+{\varphi }_{bi,2}.\hfill & \hfill \end{array}$$

(37b)

The quantities $V_1$ and $V_2$ are the externally applied voltages (units: V) at the electric contact ${\Gamma }_1$ and ${\Gamma }_2$, respectively, while ${\varphi }_{bi,1}$ and ${\varphi }_{bi,2}$ are the so-called "built-in" potentials (units: V).

The junction region JCT is characterized by the fact that $c(R,\theta )=1/\tau$, $\mathcal{P}(R,\theta )=q{p}_{eq}/\tau$, with $\tau$ (units: s) and $p_{eq}$ (units: $m^{-3}$) being given positive constants, and

$$\begin{array}{ccc}& {\displaystyle \psi (R,\theta )={\overline{\psi }}_1+\left(\frac{R-R_1}{R_2-R_1}\right)({\overline{\psi }}_2-{\overline{\psi }}_1)(R,\theta )\in [R_1,R_2]\times [0,\frac{\pi }{2}].}\hfill & \hfill \end{array}$$

Remark 5. 

After dividing q out of both sides of the balance equation (35a), we see that system (35) can be regarded as a special instance of system (3) upon setting $u=p$, $D=D_p$, $\mathbf{V}={\mu }_p\nabla \psi$ (which corresponds to ${\Psi }_K{=-\psi |}_K$), $c=1/\tau$ and $\mathcal{P}=p_{eq}/\tau$, with τ equal to $+\infty$ in ${\Omega }_1\cup {\Omega }_2$.

Let us assume that $V_2=0V$, which amounts to setting the electric contact ${\Gamma }_2$ to ground voltage. If $V_1$ is equal to $0V$, the device is working in the so-called "thermal equilibrium" condition which is characterized by a current density ${\mathbf{j}}_p\left(\mathbf{x}\right)=\mathbf{0}$ for all $(x,y)\in \Omega$. If $V_1$ is equal to a positive voltage, the device is working in the so-called "forward bias" condition, and the magnitude of ${\mathbf{j}}_p$ is large in all the computational domain. If $V_1$ is equal to a negative voltage, the device is working in the so-called "reverse bias" condition, and the magnitude of ${\mathbf{j}}_p$ is considerably smaller than in forward bias conditions, being smaller in ${\Omega }_1\cup {\Omega }_2$ than in the JCT region. In all biasing conditions, the hole concentration p is monotone in the sense of Definition 3. We refer to [26] for a physical discussion of the electrical function of the diode and to [27] for an analytical mathematical treatment.

Simulations are conducted on a nonuniform triangulation made of 870 elements and 468 vertices, with $h=8.2046·{10}^{-6}m$. The externally applied voltage at ${\Gamma }_1$ is $V_1=-5V$, corresponding to moderately high reverse bias conditions, while the values of all the physical model parameters are illustrated in Appendix A The reconstruction formula (34a) is used to determine the current density ${\mathbf{j}}_{p,h}$ in the stabilized and non-stabilized case.

Figure 4 shows the spatial distribution of $p_h$ and ${\mathbf{j}}_{p,h}$ computed by the stabilized MH method. We notice that $p_h$ is monotone and exhibits the physically expected behavior: holes move by diffusion in the neutral regions ${\Omega }_1$ and ${\Omega }_2$, while they are drifted by the large electric field in the JCT region towards the electric contact ${\Gamma }_1$. The presence of the steep internal layer is the result of the huge difference between the very small hole concentration in the JCT region (due to the consumption term) and the very large hole concentration in the neutral region ${\Omega }_1$. The streamlines of the hole current density confirm the physical picture described by the plot of the hole concentration. In particular, we notice see that the magnitude of ${\mathbf{j}}_{p,h}$ is small in all the device (because of reverse biasing conditions), except locally in the annular region (where the electric field has a very high strength) and around the point ${(a,0)}^T$ where the so-called "current-crowding" effect occurs due to an abrupt change from nonhomogeneous Dirichlet boundary condition to homogeneous Neumann boundary condition (see [14,35c,35d]).

Figure 5 shows the spatial distribution of $p_h$ and of the streamlines of ${\mathbf{j}}_{p,h}$ computed by the non-stabilized MH method. The hole approximate concentration $p_h$ is not monotone and exhibits significant oscillations between the values $-4.1599·{10}^{22}{m}^{-3}$ and $3.8747·{10}^{22}{m}^{-3}$. The unphysical behavior of the solution computed by the non-stabilzed method is confirmed by the quiver plot of ${\mathbf{j}}_{p,h}$ which is affected by vortices throughout all the computational domain.

9.3. The Case $n=3$

In this section we study the DAR model problem (3) in the domain $\Omega ={(0,1)}^3$ subject to homogeneous Dirichlet boundary conditions. The values of the parameters are: $D={10}^{-6}$, $\mathcal{P}=1$ and $c=0$, while advection is described by the counterclockwise rotating field $\mathbf{V}=-{[y,x,0]}^T$. A rather coarse uniform partition ${\mathcal{T}}_h$ of the domain is generated by subdividing each cube side into 5 intervals of length $h_x=h_y=h_z=1/5=0.2$, resulting into a number of 216 vertices and 750 tetrahedra. The mathematical model under investigation can be regarded as the three-dimensional extension of the first numerical example studied in [9, Section 7]. Figure 6 illustrates the approximate solution $u_h$ and the streamlines of ${\mathbf{j}}_h$ computed by the stabilized MHFEM, plotted in the $x-y$ plane at $z=0.5$. Results indicate the rotation effect due to the advaction field, with the formation of sharp boundary layers along the sides $x=0$ and $y=1$, without the onset of any spurious oscillations. In contrast to this, Figure 7 illustrates the approximate solution $u_h$ and the streamlines of ${\mathbf{j}}_h$ computed by the non-stabilized MHFEM, plotted in the $x-y$ plane at $z=0.5$. Results indicate the non-monotonicity of $u_h$ and the formation of local flow patterns which are totally inconsistent with the physically expected spatial distribution of the flux density vector field.

10. Conclusions and Perspectives

In this article we investigated the solution of the diffusion-advection-reaction boundary value problem. The emphasis was posed to the situation where advection and/or reaction dominate over diffusion because in this case the need of stabilization is particularly demanding in order to enforce the property of monotonicity on the discrete level, thereby preventig the occurrence of spurious unphysical oscillations in the solution computed by a numerical scheme. Conservation was also considered as an important property to be maintained on the discrete level by the numerical scheme chosen for approximation. Guided by these two principal objectives we illustrated the construction of a mixed hybrid finite element method of lowest order endowed with an edge-based stabilization which can be applied to any spatial dimension $n=1,2,3$. The method was analyzed and numerically validated in several case studies characterized by different choices of the boundary conditions. Results demonstrated the attitude of the proposed stabilized method to compute stable and accurate solutions even on coarse domain partitions. This feature of the scheme makes it particularly attractive in the discretization of time-dependent multi-dimensional nonlinear systems of DAR equations for application in computational biology, in which the cost of a simulation is directly related to the repeated solution of linearized DAR boundary value problems required to advance over time levels.