A nonlinear correction finite volume scheme preserving maximum principle for diffusion equations with anisotropic and discontinuous coefficient

1.   Introduction

A discrete scheme that satisfies fundamental properties such as conservation, positivity, and the discrete maximum principle (DMP) for the diffusion equations is of great significance. The finite volume framework can ensure local conservation of mass, an essential physical property, and is a commonly used method in constructing the discretization system (see ^[1,2,3,4]). There is a lot of research on the monotonicity of scheme, as it guarantees the positivity of the approximate solution (see ^{[5,6,7,8,9,10]}). However, a monotone scheme can only maintain the lower bound and cannot simultaneously keep the upper bound. A discrete system preserving the DMP can ensure that there are no spurious oscillations in the numerical solutions and can preserve physical bounds. Therefore, researchers have recently turned to constructing schemes that preserve the DMP (see ^{[11,12,13,14,15]}). A scheme having a specific structure ensures the DMP is developed in ^[16], which provides some theoretical analysis under the coercivity assumption. In ^[17], a nonlinear finite volume scheme preserving the DMP for the anisotropic diffusion equations on distorted meshes is presented. They give a proof of the coercivity of the scheme under some constraints on cell deformation and the diffusion coefficient.

It is challenging to design a precise and dependable scheme on distorted meshes for the diffusion problem with anisotropic and discontinuous coefficients to satisfy the DMP. As is well-known, the existing cell-centered finite volume schemes for the diffusion problem with general discontinuous coefficients on general meshes cannot unconditionally satisfy the DMP. The validity of the DMP for continuous piecewise linear finite element approximations for the Poisson equation with the Dirichlet boundary condition is verified in ^[18]. A nonlinear stabilized Galerkin approximation of the Laplace operator, which preserves the DMP on arbitrary meshes and in arbitrary space dimension without the help of the well-known acute condition or generalizations thereof, was derived in ^[19]. A mesh condition is developed for linear finite element approximation of the anisotropic diffusion–convection–reaction problem to satisfy the DMP in ^[20]. A weak Galerkin discrete scheme preserving the DMP for the boundary value problem of a general anisotropic diffusion problem is created in ^[21]. From all of the above methods, we find that there are invariably certain restrictions on the geometry of mesh-cell or the location of discontinuity in order to preserve the DMP. Desirably, our new scheme can satisfy the DMP unconditionally.

For constructing a scheme satisfying monotone or the DMP, many researchers take a measure of making a correction on flux. In ^[22], a nonlinear method to correct a general finite volume scheme for the anisotropic diffusion problem, which preserves the DMP, was presented. A monotonicity correction was applied to second-order element finite volume methods, obtaining a second-order monotone finite volume scheme for the anisotropic diffusion problem in ^[23]. In ^[24], a new nonlinear method is constructed to preserve the DMP for the diffusion problem with the heterogeneous anisotropic coefficient. Following this idea, we propose a nonlinear correction for the finite volume method to obtain a scheme that preserves the DMP, where the interpolation coefficients of the cell-edge unknowns are not required to be nonnegative. This is an advantage of our scheme.

In this paper, we propose a new nonlinear correction finite volume scheme to preserve the DMP for diffusion equations with anisotropic and discontinuous coefficients. First, two linear and nonconservative fluxes are given on two sides of the cell-edge, respectively. Meanwhile, cell-edge unknowns(auxiliary unknowns) are introduced. When eliminating these cell-edge unknowns by using a weighted combination of their neighboring cell-centered unknowns, we no longer require the weighted coefficient to be nonnegative, but just require the approximation to have second-order accuracy. This idea allows our scheme to adapt to solving the diffusion problem with discontinuous coefficient on universal distorted meshes. Second, the conservative flux is constructed by using a nonlinear weighted combination of these two one-sided edge linear fluxes. Due to the appearance of negative weighted coefficients in the first step, the scheme does not preserve the DMP so far. Then, we come up with a new nonlinear method to transform the conservative flux into a new one which possesses the structure of preserving the DMP. This is one of the key ideas of our work.

The remainder of this article is organized as follows. In Section 2, the construction of the nonlinear finite volume scheme is described. In Section 3, a prior estimation and the proof preserving the DMP of our scheme are proposed. In Section 4, some numerical results are presented to test the accuracy and verify the DMP. Finally, some conclusions are given.

2.   Construction of scheme

Consider the following stationary diffusion problem:

where the diffusion tensor $\kappa$ is a positive-definite matrix, $f\in L^2(\Omega)$, and $g\in C(\partial\Omega)$ are given functions. Let $\Omega$ be an open bounded polygonal domain in $R^2$ with boundary $\partial\Omega$.

2.1.   The one-side edge flux

In this subsection, we present the expression of the one-side edge flux with both cell-center unknowns and cell-edge unknowns. Denote the cell by $K$ or $L$, which also stand for the cell-center, and $\sigma = K|L$ is the common cell-edge of cell $K$ and $L$. Then, integrate (2.1) over the cell $K$ to get

where $\mathcal{E}_K$ is the set of all cell-edges of $K$ and $\mathcal{F}_{K, \sigma}$ is the continuous flux on the edge $\sigma$ whose expression is as follows:

The unit outer normal vector on the edge $\sigma$ of cell $K$(resp.$L$) is denoted by $\mathbf{n}_{K\sigma}$(resp.$\mathbf{n}_{L\sigma}$), and the unit tangential vector on the line $KM_i$(resp.$LM_i$)($i = 1, 2$) is denoted by $\mathbf{t}_{KM_i}$(resp.$\mathbf{t}_{LM_i}$). These notations are shown in Figure 1.

Since $KM_1$ and $KM_2$ are two edges of triangle $KM_1M_2$, the two vectors $\mathbf{t}_{KM_1}$ and $\mathbf{t}_{KM_2}$ cannot be collinear (see Figure 1). Then, there is

Similarly, there is

Substituting (2.5) and (2.6) into (2.4), respectively, we obtain

where $h = \bigg(\sup\limits_{K\in\mathcal{T}}m(K)\bigg)^{\frac{1}{2}}$ and $m(K)$ is the area of cell $K$. Let

The above two expressions (2.7) and (2.8) contain the cell-edge unknowns, which are auxiliary unknowns in this scheme. Then, we will eliminate those cell-edge unknowns. According to experience, those cell-edge unknowns can be locally approximated by some neighboring cell-center unknowns. Therefore, the approximation expression of cell-edge unknowns $u_{M_i}$ can be usually written as

where $u_{K_{M_{i}, j}}$ are some cell-center unknowns, $\omega_{M_{i}, j}$ are some corresponding coefficients and $J_{M_{i}, n}$ is the number of cell-center unknowns involved in the approximation of $u_{M_i}$. The method for selecting the interpolation mainly depends on two aspects in this scheme. On the one hand, the approximation order is required to be second order to preserve the accuracy of the scheme, which is relatively easier to achieve. On the other hand, the interpolation coefficients are allowed to be negative to enable the scheme to be applied to diffusion equations with anisotropic and discontinuous diffusion coefficients on any arbitrary distorted mesh. As is well-known, it is a challenge to realize this goal when constructing a scheme preserving the DMP.

Rewrite the expression (2.7) to

where $\bar{F_1}$ dose not include the term $u_L$, and $a_K$ is a constant which is negative or nonnegative. Similarly, the expression (2.8) can be rewritten to

where $\bar{F_2}$ dose not include the term $u_K$, and $a_L$ is a constant which is negative or nonnegative like $a_K$.

Set $a_{\sigma}$ = $\min\{|a_{K}|, |a_{L}|\}$, which is a nonnegative constant, and we can rewrite (2.10) and (2.11) to be

Therefore, the flux $\mathcal{F}_{K, \sigma}$ and $\mathcal{F}_{L, \sigma}$ can be rewritten as follows:

2.2.   The definition of conservative flux

If $|\hat{F_1}|\leq\varepsilon_1$ and $|\hat{F_2}|\leq\varepsilon_1$, where $\varepsilon_1$ is a sufficiently small positive constant, $\varepsilon_1 < Ch^2$, and the flux $\mathcal{F}_{K, \sigma}$ and $\mathcal{F}_{L, \sigma}$ can be rewritten to be

Hence, the conservative flux can be easily obtained as follows:

If $|\hat{F_1}| > \varepsilon_1$ or $|\hat{F_2}| > \varepsilon_1$, the construction of the conservative flux becomes complicated. We will make a discussion in this case. Let the discrete flux on edge $\sigma$ to be

where $\lambda_1$ and $\lambda_2$ are some positive coefficients satisfying $\lambda_1+\lambda_2 = 1$, which will be determined later. It is obvious that there is $F_{K, \sigma}+F_{L, \sigma} = 0$, which means that the discrete flux on the cell edge $\sigma$ is surely conservative.

Furthermore, there is $|\hat{F_1}|+|\hat{F_2}|\neq0$ under the circumstances. Then, we can take

Hence, the positive coefficients $\lambda_1$ and $\lambda_2$ have been obtained.

If $\hat{F_1}\hat{F_2} < 0$, that is, $|\hat{F_1}|\hat{F_2} = -|\hat{F_2}|\hat{F_1}$, then, Eqs (2.16) and (2.17) can be rewritten as follows:

When eliminating those cell-edge unknowns, the appearance of negative coefficients in (2.9) can potentially lead to negative coefficients in $2\lambda_1\hat{F_1}$ of (2.18). However, nonnegative coefficients are crucial for constructing a scheme that satisfes the DMP. To solve this problem, we present a nonlinear correction to the expression of the conservative flux. Let

where $\mathcal{J}_{K}$ is the set of all cells which share any vertex with the cell $K$.

$\mathbf{Case 1}$: There exist two cell-centered unknowns $u_{K^{'}}$ and $u_{L^{'}}$ such that

where $K^{'}\in\mathcal{J}_{K}$ and $L^{'}\in\mathcal{J}_{L}$. In reality, $u_{K^{'}}$ can be taken to be $u_{K_1}$ or $u_{K_2}$, depending on the circumstances. Similarly, $u_{L^{'}}$ can be treated in the same way.

Therefore, we can rewrite the Eqs (2.18) and (2.19) to be as follows:

where

According to (2.20), we can see that there is

Similarly, there is

$\mathbf{Case 2}$: There does not exist two cells $K^{'}$ and $L^{'}$ such that both Eqs (2.20) and (2.21) hold simultaneously. That is,

for any $K^{'}\in\mathcal{J}_{K}$, or

for any $L^{'}\in\mathcal{J}_{L}$.

Therefore, we give the following conservative flux,

If $\hat{F_1}\hat{F_2}\geq0$, that is, $|\hat{F_1}|\hat{F_2} = |\hat{F_2}|\hat{F_1}$, then, the Eqs (2.16) and (2.17) can be rewritten as follows:

When $\sigma\in\partial\Omega$, we can give the following expression of flux:

where $u_{K, i}$ are cell-center unknowns associated with cell $K$ and $u_{M, j}$ are cell-edge unknowns located on the boundary $\partial\Omega$ related to cell $K$.

2.3.   The finite volume scheme

The finite volume scheme is as follows:

where $\mathcal{P}_{in}$ is the set of all cell centers which are in the domain $\Omega$ and $\mathcal{P}_{out}$ is the set of all midpoints of boundary edge.

Obviously, the scheme results in a nonlinear system. We solve it by using the Picard iteration, and we linearize the flux as follows:

(ⅰ) when $|\hat{F_1}|\leq\varepsilon_1$ and $|\hat{F_2}|\leq\varepsilon_1$,

(ⅱ) when $|\hat{F_1}| > \varepsilon_1$ or $|\hat{F_2}| > \varepsilon_1$, $\hat{F_1}\hat{F_2} < 0$ for Case 1,

(ⅲ) when $|\hat{F_1}| > \varepsilon_1$ or $|\hat{F_2}| > \varepsilon_1$, $\hat{F_1}\hat{F_2} < 0$ for Case 2,

(ⅳ) when $|\hat{F_1}| > \varepsilon_1$ or $|\hat{F_2}| > \varepsilon_1$, $\hat{F_1}\hat{F_2}\geq0$,

Here, $s$ is the nonlinear iteration index.

Let $A(U)$ be the matrix of this system, where $U$ is the vector of discrete unknowns. The global discrete nonlinear system reads as:

We choose a small value $\varepsilon_{non} > 0$ and initial vector $U^0$ in the Picard iteration, and repeat for $s = 1, 2, \ldots$,

1) solve $A(U^{s-1})U^s = B$,

2) stop if $\|A(U^s)U^s-B\|\leq\varepsilon_{non}$.

In this paper, we take $\varepsilon_{non} = 10^{-5}$.

3.   Analysis

In this section, we give the prior estimation and analysis of maximum principle. Let $\mathcal{E}$ be the set of all cell-edges. Define the norms

Before giving the prior estimation, we need to give the assumption (H1) that the scheme (2.30) and (2.31) satisfies the coercivity property. That is, the following inequality holds:

where $C_1$ is a constant independent of $h$.

3.1.   The prior estimation

For the nonlinear correction scheme in this paper, we can get the priori estimation.

$\mathbf{Theorem}$ $\mathbf{1.}$ Under the assumption (H1), for any solution $u$ of the scheme (2.30) and (2.31), there holds

where $C_3$ is a positive constant independent of $h$.

Proof. Multiply (2.30) by $u_K$ and sum up these products for all $K$ to obtain

Noticing the conservation of the normal flux, we can get

Using the assumption (H1), we can get the following inequality:

Hence, there is

Using the discrete Poincare inequality,

We can obtain the inequality as follows:

Then, it follows that

□

3.2.   The discrete maximum principle

From the structure of $F_{K, \sigma}$ in Eqs (2.14), (2.22), (2.26) and (2.28), we figure out that the discrete flux has the form:

where $N_{K, \sigma}$ is the number of cell-center unknowns involved in the flux $F_{K, \sigma}$. We can easily know that $A_{K, \sigma, i}\geq0$. Then, the finite volume scheme (2.30) reads as:

where $N_{K}$ is the number of cell-center unknowns related to cell $K$, and $A_{K, i}\geq0$, which is the sum of some $A_{K, \sigma, i}$. So far, we conclude that our scheme satisfies the DMP.

$\mathbf{Theorem}$ $\mathbf{2.}$ The finite volume scheme (2.30) and (2.31) preserves the discrete maximum principle. Let $u_{min} = \min\limits_{K\in\mathcal{P}_{in}\cup\mathcal{P}_{out}}u_K$ and $u_{max} = \max\limits_{K\in\mathcal{P}_{in}\cup\mathcal{P}_{out}}u_K$, then $u_{min}$ and $u_{max}$ are only obtained on the boundary unless $u$ is a constant in the whole domain. That is, $u_{min}$ can be only obtained on the boundary of domain unless $u$ is a constant for $f\geq0$, and $u_{max}$ can be only obtained on the boundary of domain unless $u$ is a constant for $f\leq0$.

The proof of this theorem is similar to some existing papers, so we do not prove it any more.

4.   Numerical experiment

As usual, we use discrete $L_2-$norms to evaluate approximate errors. The $L_2-$norm for the solution $\mathrm{u}$ is taken to be $\varepsilon_2^\mathrm{u} = \bigg[\sum_{K\in\mathcal{T}} (\mathrm{u}_K-\mathrm{u}(K))^2m(K)\bigg]^{\frac{1}{2}}$ and the $L_2-$norm for the flux $\mathrm{F}$ is taken to be $\varepsilon_2^\mathrm{F} = \bigg[\sum_{K\in\mathcal{T}} (\mathrm{F}_{K, \sigma}-\mathcal{F}_{K, \sigma})^2\bigg]^{\frac{1}{2}}$.

4.1.   Test 1: The linear elliptic equation with continuous coefficients

Test the linear elliptic equation with continuous coefficient in unit square $\Omega = (0, 1)^2$. The diffusion coefficient is

where $k_1 = 1+2x^2+y^2$, $k_2 = 1+x^2+2y^2$, and $\theta = \frac{5\pi}{12}$. The exact solution is taken as $u(x, y) = sin(\pi x)sin(\pi y)$.

We apply our scheme and the scheme in ^[13] to solve this problem on quadrilateral meshes shown in Figure 2 and triangular meshes shown in Figure 3. Tables 1 and 2 exhibit the error and convergence of these two schemes for test 1 on these two kinds of meshes, respectively. From these two tables, it's easy to find that the convergence rate of our new scheme for the solution is almost second-order and the flux is more than first-order. Meanwhile, we discover that the error of the solution of our new method is similar to the scheme in ^[13] by this test example, but the accuracy of the flux remarkably reduced for our new scheme compared to the scheme in ^[13].

4.2.   Test 2: The linear elliptic equation with discontinuous coefficients

Consider the problem (2.1) and (2.2) with discontinuous coefficients in the unit square $\Omega = (0, 1)^2$, where

and

The exact solution is

By using our scheme and the scheme in ^[13], we test the discontinuous problem on random quadrilateral and triangular meshes with a discontinuity given in Figures 4 and 5. The errors are exhibited in Tables 3 and 4, respectively. We can easily know that our scheme gains nearly second-order accuracy for the solution and more than first-order accuracy for the flux on both kinds of meshes.

From Table 3, we can see that the accuracy of the solution of our scheme reduces a little compared to the scheme in ^[13], but the accuracy of the flux reduces a lot. From Table 4, the accuracy of the solution and the flux of our new scheme are both remarkably higher than the scheme in ^[13]. Overall, our new nonlinear scheme is relatively better suited for solving diffusion equations with anisotropic and discontinuous coefficients on distorted meshes.

4.3.   Test 3: The linear elliptic equation with linear solution

Consider the problem (2.1) and (2.2) with diffusion coefficients

and source term $f = 0$. Take the following linear solution:

We set the stopping criterion $\varepsilon_{non} = 10^{-14}$ for this example, and test it on random quadrilateral meshes and random triangular meshes, respectively. Computational results are shown in Tables 5 and 6. Both tables demonstrate that our scheme can accurately recover the linear solution.

4.4.   Test 4: The linear elliptic equation on the domain with a hole

Consider the problem (2.1) and (2.2) on the domain $\Omega = (0, 1)^2\backslash[4/9, 5/9]^2$, whose boundary $\partial\Omega$ is composed of two disjoint parts $\Gamma_1$ and $\Gamma_0$ exhibited in Figures 6 and 7, where $\Gamma_1$ is the interior boundary and $\Gamma_0$ is the exterior boundary. Set $f = 0$, $g = 0$ on $\Gamma_0$, and $g = 2$ on $\Gamma_1$. Take the anisotropic diffusion tensor $\kappa$ as follows:

where $k_1 = 1$, $k_2 = \left\{ \begin{array}{lcl} 100, \quad (x, y)\in(0, \frac{2}{3}]\times(0, 1)\\ 10, \quad (x, y)\in(\frac{2}{3}, 1)\times(0, 1) \end{array} \right.$, and $\theta = \frac{\pi}{6}$.

We test this example on random quadrilateral meshes with a hole shown in Figure 6 and random triangular meshes with a hole shown in Figure 7. The contour of the numerical solution of the scheme in ^[25] on these two kinds of meshes are shown in Figures 8 and 9, respectively. The analysis of both figures reveals that the numerical solution exhibits negative minima and surpasses the value of 2 at its maxima. Consequently, the scheme in ^[25] violates the DMP.

The numerical solution of our new scheme on these two kinds of meshes are shown in Figures 10 and 11, respectively. Both figures show that the numerical solution maintains a minimum value of 0 and reaches a maximum value of 2. These results testify that our scheme satisfies the DMP.

4.5.   Test 5: The linear elliptic equation with anisotropic solution

Consider the problem (2.1) and (2.2) with discontinuous and anisotropic diffusion tensor coefficients as follows:

The computational domain is $\Omega = (0, 16)\times(0, 16)$. We set $f = 0$, $g(x, 0) = 0$, $g(16, y) = 0$,

and

Though the solution of the problem is unknown, the maximum principle states that the solution of this problem should stay between $0$ and $1$.

We first use the positivity-preserving scheme in ^[26] to test this example on the random quadrilateral meshes with $48\times48$ cells and the random triangular meshes with $48\times48\times2$ cells. The contour of the numerical solution on random quadrilateral meshes is shown in Figure 12, and that on random triangular meshes is shown in Figure 13. Both figures reveal that the peak value of the numerical solution surpasses 1. The numerical findings indicate that the positivity-preserving scheme fails to maintain the upper bound constraints of the solution.

We use our nonlinear scheme to test this example on the same two kinds of meshes. Results are presented in Figures 14 and 15, respectively. From these two figures, we see that our new scheme avoids undershoot and overshoot, which testify that our scheme satisfies the DMP.

5.   Conclusions

An improved nonlinear finite volume scheme which preserves the DMP for the diffusion problem with anisotropic and discontinuous coefficients is developed. This new scheme can solve diffusion equations with the discontinuous coefficient on any arbitrary distorted meshes, overcoming the drawback of existing schemes that require the interpolation coefficients of cell-edge unknowns to be nonnegative. We provide a prior estimation based on an assumption of the coercivity property. Numerical results demonstrate that the accuracy of this new scheme surpasses that of the existing scheme. Furthermore, numerical results testify that this new scheme satisfies the DMP. When solving the nonlinear algebraic system, we find that the Picard iteration is unsatisfactory in terms of computational efficiency. We hope to detect a useful iteration method that is adapted for this system in the future.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is partially supported by Shandong Provincial Natural Science Fundation (ZR2021QA109), National Natural Science Foundation of China (No. 12101536) and the LCP Fund for Young Scholar (No. 6142A05QN22004).

Conflict of interest

The authors declare there is no conflicts of interest.