Efficient Numerical Methods for Strongly Anisotropic Elliptic Equations

Abstract

In this paper, we study an efficient numerical scheme for a strongly anisotropic elliptic problem which arises, for example, in the modeling of magnetized plasma dynamics. A small parameter ε induces the anisotropy of the problem and leads to severe numerical difficulties if the problem is solved with standard methods for the case 0<ε≪1. An Asymptotic-Preserving scheme is therefore introduced in this paper in a 2D framework, with an anisotropy aligned to one coordinate axis and an ε-intensity which can be either constant or variable within the simulation domain. This AP scheme is uniformly precise in ε, permitting thus the choice of coarse discretization grids, independent of the magnitude of the parameter ε.

Appendix: 1D Simulations for Heterogeneous Anisotropy Ratios with Steep Gradients

In this section, we detail different procedures to deal with huge variations in the anisotropy intensities. A slightly modified 1D framework is investigated in order to simplify the presentation of the different approaches used in Sect. 3.3 for the 2D simulations. The 1-dimensional SP-model now reads

$$ \left \{ \begin{array}{l@{\quad}l} -\displaystyle\frac{\mathsf{d}}{\mathsf{d}z} \biggl(\frac{1}{\varepsilon (z)}\frac{\mathsf{d}u}{\mathsf{d}z} \biggr)+u=f,&\text{in }\varOmega_z,\\[12pt] \displaystyle\frac{\mathsf{d}}{\mathsf{d}z}u=0,&\text{on }\partial\varOmega_z. \end{array} \right . $$

(A.1)

where u is the unknown of the problem, ε is a positive function of z and Ω _z =[0,L _z ]. Note that Eq. (A.1) is well-posed for ε(z)>0. The uniqueness of the solution in the 1D case is achieved by adding a term in the equation. This was not needed in the 2D case, due to the presence the derivative in the second direction x. Taking now ε(z):=δχ(z) and passing to the limit δ→0, yields the degenerate problem

$$ \left \{ \begin{array}{l@{\quad}l} -\displaystyle\frac{\mathsf{d}}{\mathsf{d}z} \biggl(\displaystyle\frac{1}{\chi(z)}\frac{\mathsf{d}u}{\mathsf{d}z} \biggr)=0,&\text{in}\ \varOmega_z,\\[12pt] \displaystyle\frac{\mathsf{d}}{\mathsf{d}z}u=0,&\text{on}\ \partial\varOmega_z, \end{array} \right . $$

(A.2)

which is ill-posed, since all constants are solutions.

A standard AP-discretization is compared with a non conservative AP-scheme as well as a Harmonic-Mean and a Scharfetter-Gummel AP-scheme.

1.1 A.1 The Standard AP-Scheme

To construct the AP-scheme corresponding to this SP-model, let us decompose u into its mean part $\bar{u}$ and its fluctuating part u′ and follow the guideline of Sect. 2.1. By integrating equation (A.1) in Ω _z , we get the average equation $\bar{u}=\bar{f}$. Subtracting this equation from (A.1) and introducing the Lagrange multiplier λ, we obtain the following AP-formulation

$$ \left \{ \begin{array}{l} \bar{u}=\bar{f},\\[6pt] \displaystyle\int^{L_z}_0 \biggl[\frac{1}{\varepsilon (z)}\frac{\mathsf{d}u'}{\mathsf{d}z}\frac{\mathsf{d}v'}{\mathsf{d}z}+u'v' \biggr] \mathsf{d}z+\lambda\int^{L_z}_0\frac{1}{\varepsilon (z)}v'\mathsf{d}z =\displaystyle\int^{L_z}_0f'v'\mathsf{d}z,\\[12pt] \quad \quad \forall v'\in H^1(\varOmega_z),\\[12pt] \displaystyle\int^{L_z}_0 u'\mathsf{d}z=0. \end{array} \right . $$

(A.3)

This AP scheme will be discretized by a ℙ₁ finite element method. However, note that the difficulty here is the approximation of the high anisotropy-gradients. To face this problem, we shall propose here different methods.

1.2 A.2 Non-Conservative AP-Scheme

In this approach, we try to circumvent the high anisotropy ratio term $\frac{1}{\varepsilon }$ in the fluctuation equation of (A.3). This is done, by developing $\frac{\mathsf{d}}{\mathsf{d}z} (\frac{1}{\varepsilon (z)}\frac{\mathsf{d}u'}{\mathsf{d}z} )$, in order to obtain

$$ \frac{\mathsf{d}}{\mathsf{d}z} \bigl(\ln{ \varepsilon (z)} \bigr)\frac{\mathsf{d}u'}{\mathsf{d}z}- \frac{\mathsf{d}^2u'}{\mathsf{d}z^2}+\varepsilon u'=\varepsilon f'. $$

(A.4)

By injecting (A.4) in the AP-reformulation, we get another reformulation, i.e.

$$ \left \{ \begin{array}{l} \bar{u}=\bar{f},\\[6pt] \displaystyle\int^{L_z}_0 \biggl[\frac{\mathsf{d}}{\mathsf{d}z} (\ln{ \varepsilon (z)} )\frac{\mathsf{d}u'}{\mathsf{d}z}v'+\frac{\mathsf{d}u'}{\mathsf{d}z}\frac {\mathsf{d}v'}{\mathsf{d}z}+\varepsilon u'v' \biggr]\mathsf{d}z+\lambda\int^{L_z}_0v'\mathsf{d}z =\int^{L_z}_0\varepsilon f'v'\mathsf{d}z,\\[12pt] \displaystyle\quad \forall v'\in H^1(\varOmega_z),\\[12pt] \displaystyle\int^{L_z}_0 u'\mathsf{d}z=0. \end{array} \right . $$

(A.5)

Again, we shall discretize (A.5) by the ℙ₁ finite element method.

1.3 A.3 Harmonic Mean AP-Scheme

The harmonic mean AP scheme is just a special discretization of the AP-formulation (A.3) following the ideas detailed in [21, 23]. To present this approach, let us again consider the partition of Ω _z and the basis functions defined in Sect. 2.3.1. Taking now in (A.3) as test functions κ _k (z) and replacing u′(z) by $\sum^{N_{z}+1}_{l=0}\alpha_{l}\kappa_{l}(z)$ give rise for k=1,…,N _z to

where

$$p_{k}=\int^{z_{k}}_{z_{k-1}}\frac{1}{\varepsilon (z)}\mathsf{d}z, \qquad q_k=\int^{z_{k+1}}_{z_{k-1}}\frac{1}{\varepsilon (z)}\kappa_k(z)\mathsf{d}z,\qquad f_k=\int^{z_{k+1}}_{z_{k-1}}f'(z)\kappa_k(z)\mathsf{d}z. $$

If we discretize p _k by a standard numerical quadrature formula, the obtained scheme will just be a ℙ₁ finite element method. However, we will use instead a harmonic mean to approximate p _k , i.e.

$$ p_k\approx\Delta z^2 \biggl(\int^{z_{k}}_{z_{k-1}}{ \varepsilon (z)} \mathsf{d}z \biggr)^{-1}. $$

Finally, the integration $\int^{z_{k}}_{z_{k-1}}{\varepsilon (z)}\mathsf{d}z$ is approximated by a standard numerical quadrature, for example Gauss-Legendre quadrature. By this manner, we obtain the full harmonic mean AP scheme.

1.4 A.4 Scharfetter-Gummel AP-Scheme

The Scharfetter-Gummel AP-scheme is an amelioration of the harmonic mean AP-scheme. In particular a special quadrature formulae is used to approximate p _k [21]. Denoting ε(z _k ) simply by ε _k , we approximate p _k as follows

• if ε _k−1≠ε _k ,

  

• if ε _k−1=ε _k ,

  $$ p_k\approx\Delta z\frac{1}{\varepsilon _k}. $$

  where lnε is approximated by $\sum^{N_{z}+1}_{k=0}\ln \varepsilon _{k}\kappa_{k}$.

1.5 A.5 Numerical Results

To compare these approaches, we take an exact solution of (A.1) defined as

$$ u_e(z):=\varepsilon (z)\cos \biggl(\frac{2\pi}{L_z}z \biggr), $$

(A.6)

with ε(z) defined by expression (3.5). We shall compare the SP-model, the standard AP scheme, the Non-Conservative AP scheme, the Harmonic Mean AP scheme and the Scharfetter-Gummel AP scheme respectively. For this, we choose q=80, ε _max=1 and different values for ε _min. The relative error of the approximations with the exact solution are gathered in Fig. 6. The results reported in the Fig. 6(c) show that the SP-model can not furnish an accurate approximation of the solution. On the contrary the curves related to the approximation computed by the AP-schemes displayed on Figs. 6(a) and 6(b) demonstrate the robustness of these schemes. The standard AP-scheme as well as the HM-AP scheme appear to be rather unprecise in the region near the anisotropy steep gradients, while the NC-AP scheme exhibits oscillations in the same regions. Finally, the SG-AP scheme provides the best approximation above the different methods proposed.

Fig. 6

Comparison between a standard discretization of the SP-model, the AP scheme, the non-conservative AP (NC-AP) scheme, the harmonic mean AP (HM-AP) scheme and the Scharfetter-Gummel AP (SG-AP) scheme on a grid with 25 points. (a) Plots of the approximate solutions computed thanks to the different methods for ε _min=10⁻¹; (b) Same plots for ε _min=10⁻¹⁰; (c) Relative errors between the exact solution and the approximate ones for both cases ε _min=10⁻¹ and ε _min=10⁻¹⁰