A positivity-preserving finite volume element method for anisotropic diffusion problems on quadrilateral meshes

TL;DR

This paper introduces a nonlinear finite volume element scheme that preserves positivity for anisotropic diffusion problems on quadrilateral meshes, demonstrating high accuracy and efficiency even on distorted meshes and with discontinuous coefficients.

Contribution

A novel positivity-preserving nonlinear FVE scheme using a two-point flux technique for anisotropic diffusion on quadrilateral meshes, avoiding complex decompositions.

Findings

Effective on distorted quadrilateral meshes

Achieves approximate second-order accuracy

More efficient than standard FVE methods

Abstract

In this paper, we propose a nonlinear positivity-preserving finite volume element(FVE) scheme for anisotropic diffusion problems on quadrilateral meshes. Based on an overlapping dual partition, the one-sided flux is approximated by the iso-parametric bilinear element. A positivity-preserving nonlinear scheme with vertex-centered unknowns is obtained by a new two-point flux technique, which avoids the convex decomposition of co-normals and the introduction of intermediate unknowns. The existence of a solution is proved for this nonlinear system by applying the Brouwer's theorem. Numerical results show that the proposed positivity-preserving scheme is effective on distorted quadrilateral meshes and has approximate second-order accuracy for both isotropic and anisotropic diffusion problems. Moreover, the presented scheme is applied on an equilibrium radiation diffusion problem with discontinuous coefficients. Numerical results show that the new scheme is much more efficient than the standard FVE method.

Tables

Table 1. Table 1: Standard FVE method: Errors and convergence rates for 𝐄𝐱𝐚𝐦𝐩𝐥𝐞𝟏 𝐄𝐱𝐚𝐦𝐩𝐥𝐞𝟏 \mathbf{Example1} .

	Uniform meshes				Random distorted meshes
Mesh level	${\Vert u-u_h\Vert }_{0,h}$	Order	${\mid u-u_h\mid }_{1,h}$	Order	${\Vert u-u_h\Vert }_{0,h}$	Order	${\mid u-u_h\mid }_{1,h}$	Order
1	$4.3080\times {10}^{-3}$	$-$	$1.2970\times {10}^{-2}$	$-$	$4.1090\times {10}^{-3}$	$-$	$1.7767\times {10}^{-2}$	$-$
2	$1.0642\times {10}^{-3}$	2.0172	$4.0432\times {10}^{-3}$	1.6816	$1.3476\times {10}^{-3}$	1.6985	$8.0991\times {10}^{-3}$	1.1968
3	$2.6630\times {10}^{-4}$	1.9986	$1.2056\times {10}^{-3}$	1.7457	$3.1712\times {10}^{-4}$	2.1088	$3.7088\times {10}^{-3}$	1.1384
4	$6.6784\times {10}^{-5}$	1.9954	$3.4669\times {10}^{-4}$	1.7980	$1.0000\times {10}^{-4}$	1.6895	$1.8126\times {10}^{-3}$	1.0480
5	$1.6738\times {10}^{-5}$	1.9963	$9.7073\times {10}^{-5}$	1.8365	$2.7870\times {10}^{-5}$	1.8788	$9.2189\times {10}^{-4}$	0.9942

Table 2. Table 2: Monotone FVE scheme: Errors and convergence rates for 𝐄𝐱𝐚𝐦𝐩𝐥𝐞𝟏 𝐄𝐱𝐚𝐦𝐩𝐥𝐞𝟏 \mathbf{Example1} .

	Uniform meshes				Random distorted meshes
Mesh level	${\Vert u-u_h\Vert }_{0,h}$	Order	${\mid u-u_h\mid }_{1,h}$	Order	${\Vert u-u_h\Vert }_{0,h}$	Order	${\mid u-u_h\mid }_{1,h}$	Order
1	$3.1348\times {10}^{-3}$	$-$	$1.4965\times {10}^{-2}$	$-$	$4.8452\times {10}^{-3}$	$-$	$2.3929\times {10}^{-2}$	$-$
2	$8.0581\times {10}^{-4}$	1.9598	$4.7442\times {10}^{-3}$	1.6573	$1.5094\times {10}^{-3}$	1.7768	$1.0179\times {10}^{-2}$	1.3022
3	$2.0610\times {10}^{-4}$	1.9670	$1.4366\times {10}^{-3}$	1.7235	$3.7475\times {10}^{-4}$	2.0307	$4.9751\times {10}^{-3}$	1.0434
4	$5.1667\times {10}^{-5}$	1.9960	$4.1797\times {10}^{-4}$	1.7811	$1.1363\times {10}^{-4}$	1.7469	$2.4630\times {10}^{-3}$	1.0292
5	$1.2361\times {10}^{-5}$	2.0634	$1.1804\times {10}^{-4}$	1.8241	$3.0842\times {10}^{-5}$	1.9177	$1.2319\times {10}^{-3}$	1.0188

Table 3. Table 3: Standard FVE Method: Errors and convergence rates for 𝐄𝐱𝐚𝐦𝐩𝐥𝐞𝟐 𝐄𝐱𝐚𝐦𝐩𝐥𝐞𝟐 \mathbf{Example2} .

	Uniform meshes				Random distorted meshes
Mesh level	${\Vert u-u_h\Vert }_{0,h}$	Order	${\mid u-u_h\mid }_{1,h}$	Order	${\Vert u-u_h\Vert }_{0,h}$	Order	${\mid u-u_h\mid }_{1,h}$	Order
1	$5.6362\times {10}^{-3}$	$-$	$1.8126\times {10}^{-2}$	$-$	$6.8585\times {10}^{-3}$	$-$	$7.2047\times {10}^{-2}$	$-$
2	$1.3960\times {10}^{-3}$	2.0134	$4.5117\times {10}^{-3}$	2.0063	$1.8195\times {10}^{-3}$	2.0216	$3.3136\times {10}^{-2}$	1.1833
3	$3.4820\times {10}^{-4}$	2.0033	$1.1266\times {10}^{-3}$	2.0016	$4.5899\times {10}^{-4}$	2.0075	$1.5466\times {10}^{-2}$	1.1106
4	$8.7001\times {10}^{-5}$	2.0008	$2.8159\times {10}^{-4}$	2.0003	$1.1934\times {10}^{-4}$	1.9719	$7.7407\times {10}^{-3}$	1.0132
5	$2.1747\times {10}^{-5}$	2.0002	$7.0393\times {10}^{-5}$	2.0009	$3.0536\times {10}^{-5}$	2.0045	$3.8458\times {10}^{-3}$	1.0287

Table 4. Table 4: Monotone FVE scheme: Errors and convergence rates for 𝐄𝐱𝐚𝐦𝐩𝐥𝐞𝟐 𝐄𝐱𝐚𝐦𝐩𝐥𝐞𝟐 \mathbf{Example2} .

	Uniform meshes				Random distorted meshes
Mesh level	${\Vert u-u_h\Vert }_{0,h}$	Order	${\mid u-u_h\mid }_{1,h}$	Order	${\Vert u-u_h\Vert }_{0,h}$	Order	${\mid u-u_h\mid }_{1,h}$	Order
1	$1.8732\times {10}^{-2}$	$-$	$2.1715\times {10}^{-2}$	$-$	$1.0505\times {10}^{-2}$	$-$	$1.3502\times {10}^{-1}$	$-$
2	$4.6227\times {10}^{-3}$	2.0186	$5.4135\times {10}^{-3}$	2.0040	$2.8031\times {10}^{-3}$	2.0127	$6.5878\times {10}^{-2}$	1.0933
3	$1.1518\times {10}^{-3}$	2.0048	$1.3700\times {10}^{-3}$	1.9823	$8.2891\times {10}^{-4}$	1.7759	$3.7008\times {10}^{-2}$	0.8405
4	$2.8767\times {10}^{-4}$	2.0014	$3.6002\times {10}^{-4}$	1.9280	$2.4248\times {10}^{-4}$	1.7994	$1.9862\times {10}^{-2}$	0.9110
5	$7.1895\times {10}^{-5}$	2.0004	$1.0489\times {10}^{-4}$	1.7791	$7.1023\times {10}^{-5}$	1.8057	$1.0479\times {10}^{-2}$	0.9403

Table 5. Table 5: Standard FVE method: Errors and convergence rates for 𝐄𝐱𝐚𝐦𝐩𝐥𝐞𝟑 𝐄𝐱𝐚𝐦𝐩𝐥𝐞𝟑 \mathbf{Example3} .

	Uniform meshes				Random distorted meshes
Mesh level	${\Vert u-u_h\Vert }_{0,h}$	Order	${\mid u-u_h\mid }_{1,h}$	Order	${\Vert u-u_h\Vert }_{0,h}$	Order	${\mid u-u_h\mid }_{1,h}$	Order
1	$3.1045\times {10}^{-3}$	$-$	$4.0033\times {10}^{-3}$	$-$	$3.9770\times {10}^{-3}$	$-$	$2.5651\times {10}^{-2}$	$-$
2	$7.7297\times {10}^{-4}$	2.0058	$9.9901\times {10}^{-4}$	2.0026	$1.2873\times {10}^{-3}$	1.6441	$1.6449\times {10}^{-2}$	0.6293
3	$1.9304\times {10}^{-4}$	2.0015	$2.4971\times {10}^{-4}$	2.0002	$3.3317\times {10}^{-4}$	1.9787	$8.1944\times {10}^{-3}$	1.0200
4	$4.8248\times {10}^{-5}$	2.0003	$6.2500\times {10}^{-5}$	1.9983	$8.8884\times {10}^{-5}$	1.9431	$4.0906\times {10}^{-3}$	1.0217
5	$1.2063\times {10}^{-5}$	1.9998	$1.5841\times {10}^{-5}$	1.9801	$2.2002\times {10}^{-5}$	2.0020	$2.1065\times {10}^{-3}$	0.9759

Table 6. Table 6: Monotone FVE scheme: Errors and convergence rates for 𝐄𝐱𝐚𝐦𝐩𝐥𝐞𝟑 𝐄𝐱𝐚𝐦𝐩𝐥𝐞𝟑 \mathbf{Example3} .

	Uniform meshes				Random distorted meshes
Mesh level	${\Vert u-u_h\Vert }_{0,h}$	Order	${\mid u-u_h\mid }_{1,h}$	Order	${\Vert u-u_h\Vert }_{0,h}$	Order	${\mid u-u_h\mid }_{1,h}$	Order
1	$2.7071\times {10}^{-3}$	$-$	$3.0775\times {10}^{-3}$	$-$	$4.2274\times {10}^{-3}$	$-$	$3.2208\times {10}^{-2}$	$-$
2	$6.7108\times {10}^{-4}$	2.0121	$1.2207\times {10}^{-3}$	1.3340	$1.4420\times {10}^{-3}$	1.6386	$2.1269\times {10}^{-2}$	0.6322
3	$1.6554\times {10}^{-4}$	2.0193	$3.8384\times {10}^{-4}$	1.6691	$3.9136\times {10}^{-4}$	1.9009	$1.0640\times {10}^{-2}$	1.0095
4	$4.2234\times {10}^{-5}$	1.9707	$1.1397\times {10}^{-4}$	1.7518	$9.8116\times {10}^{-5}$	2.0253	$5.1350\times {10}^{-3}$	1.0665
5	$1.0510\times {10}^{-5}$	2.0066	$2.9818\times {10}^{-5}$	1.9343	$2.5751\times {10}^{-5}$	1.9672	$2.6463\times {10}^{-3}$	0.9749

Table 7. Table 7: L 2 subscript 𝐿 2 L_{2} -norm of the solution for the proposed scheme.

	$L_2$-norm on uniform meshes	$L_2$-norm on random quadrilateral meshes
Standard FVE method	0.8854	0.8868
Monotone FVE scheme	0.8854	0.8867

Table 8. Table 8: Average number of iterations and CPU time on random quadrilateral meshes.

	Number of nonlinear iterations	Number of linear iterations	CPU time
Standard FVE method	20.64075	32.654	3.85(h)
Monotone FVE scheme	20.67585	20.000	1.12(h)

Equations

$${\partial }_{t}u-\nabla \cdot (\kappa \nabla u)$$

$$\gamma (\kappa \nabla u)\cdot \mathbf{n}+\delta u$$

$$\forall (x,t)\in \partial \Omega \times (0,T),\delta (x,t)\ge 0,\text{and}\gamma (x,t)\ge {\gamma }_0>0.$$

$$u(x,0)=u_0(x)\in H^1(\Omega ).$$

$$\frac{h_K}{{\rho }_K}\le C,$$

$$\mid \cos {\theta }_K\mid \le C_1<1,\text{where}{\theta }_K\text{is any interior angle of}K.$$

$$F_{K}=\left\{\begin{array}[]{ll}x=x_{1}+a_{1}\xi+a_{2}\eta+a_{3}\xi\eta,\\ y=y_{1}+b_{1}\xi+b_{2}\eta+b_{3}\xi\eta,\end{array}\right.$$

$$a_1=x_2-x_1,a_2=x_4-x_1,a_3=x_1-x_2+x_3-x_4,$$

$$b_1=y_2-y_1,b_2=y_4-y_1,b_3=y_1-y_2+y_3-y_4,$$

$$J_{K}(\xi,\eta)=\left(\begin{array}[]{cc}\frac{\partial x}{\partial\xi}&\frac{\partial x}{\partial\eta}\\ \frac{\partial y}{\partial\xi}&\frac{\partial y}{\partial\eta}\\ \end{array}\right)=\left(\begin{array}[]{cc}a_{1}+a_{3}\eta&a_{2}+a_{3}\xi\\ b_{1}+b_{3}\eta&b_{2}+b_{3}\xi\\ \end{array}\right).$$

$${\mathcal{J}}_Q=S_K.$$

$${\mathcal{J}}_Q\ge C{h}_K^2.$$

$$U_h=\{u_h\in H^1(\overline{\Omega }):u_h{\mid }_K={\hat{u}}_h\circ F_K^{-1},{\hat{u}}_h\in U_h(\hat{K}),K\in {\mathcal{T}}_h\},$$

$$V_h=\{v_h\in L^2(\Omega ):v_h{\mid }_{K_P^{\ast }}=\text{constant},\forall P\in {\Omega }_h\}.$$

$$({\partial }_t{u}_h,v_h)+a_h(u_h,v_h)=(f,v_h)\forall v_h\in V_h$$

$$a_h(u_h,v_h)=-\sum _{K_P^{\ast }\in {\mathcal{T}}_h^{\ast }}{\int }_{\partial K_P^{\ast }}\kappa \nabla u_h\cdot \mathbf{n}{v}_{h}ds=-\sum _{K_P^{\ast }\in {\mathcal{T}}_h^{\ast }}{v}_h(P){\int }_{\partial K_P^{\ast }}\kappa \nabla u_h\cdot \mathbf{n}ds,$$

$$(f,v_h)=\sum _{K_P^{\ast }\in {\mathcal{T}}_h^{\ast }}{\int }_{K_P^{\ast }}f{v}_{h}dx=\sum _{K_P^{\ast }\in {\mathcal{T}}_h^{\ast }}{v}_h(P){\int }_{K_P^{\ast }}fdx,$$

$$({\partial }_{t}u,v_h)+\sum _{K_P^{\ast }\in {\mathcal{T}}_h^{\ast }}{v}_h(P)\sum _{\sigma \in \partial K_P^{\ast }}{\mathcal{F}}_{P,\sigma }=\sum _{K_P^{\ast }\in {\mathcal{T}}_h^{\ast }}{v}_h(P){\int }_{K_P^{\ast }}fdx\forall v_h\in V_h$$

$${\mathcal{F}}_{P,\sigma }=-{\int }_{\sigma }\kappa \nabla u\cdot \mathbf{n}ds=-{\int }_{\sigma }\nabla u\cdot {\kappa }^T\mathbf{n}ds$$

$$u_h(x,y)={\hat{u}}_h(\xi ,\eta )=u_1{N}_1(\xi ,\eta )+u_2{N}_2(\xi ,\eta )+u_3{N}_3(\xi ,\eta )+u_4{N}_4(\xi ,\eta )$$

$$N_1(\xi ,\eta )=(1-\xi )(1-\eta ),N_2(\xi ,\eta )=\xi (1-\eta )$$

$$N_3(\xi ,\eta )=\xi \eta ,N_4(\xi ,\eta )=(1-\xi )\eta .$$

$$\nabla u_h{\mid }_Q=\frac{u_1-u_3}{2{\mathcal{J}}_Q}{\mathbf{q}}_1+\frac{u_2-u_4}{2{\mathcal{J}}_Q}{\mathbf{q}}_2,$$

$$F_{P_1}$$

$$A_{\sigma }=\frac{{\mathbf{q}}_1\cdot {\kappa }_Q^T{\mathbf{q}}_1}{2{\mathcal{J}}_Q}\text{is the main part of the flux,}$$

$$r_{\sigma }(u_2,u_4)=\frac{{\mathbf{q}}_2\cdot {\kappa }_Q^T{\mathbf{q}}_1}{2{\mathcal{J}}_Q}(u_2-u_4)\text{is the residual of the flux,}$$

$$0<A_{\sigma }\le \frac{{\lambda }_{\text{max}}{h}_K^2}{C{h}_K^2}=\frac{{\lambda }_{\text{max}}}{C},$$

$$F_{P_3}$$

$$r_{\sigma }^{+}=\frac{\mid r_{\sigma }(u_2,u_4)\mid +r_{\sigma }(u_2,u_4)}{2},r_{\sigma }^{-}=\frac{\mid r_{\sigma }(u_2,u_4)\mid -r_{\sigma }(u_2,u_4)}{2}.$$

$$F_{P_1}=A_{\sigma }(u_1-u_3)+r_{\sigma }^{+}-r_{\sigma }^{-},$$

$$F_{P_3}=A_{\sigma }(u_3-u_1)+r_{\sigma }^{-}-r_{\sigma }^{+}.$$

$$F_{P_1,\sigma }=(A_{\sigma }+\frac{r_{\sigma }^{+}M}{M{u}_1+C{h}^2})u_1-(A_{\sigma }+\frac{r_{\sigma }^{-}M}{M{u}_3+C{h}^2})u_3.$$

$$F_{P_3,\sigma }=(A_{\sigma }+\frac{r_{\sigma }^{-}M}{M{u}_3+C{h}^2})u_3-(A_{\sigma }+\frac{r_{\sigma }^{+}M}{M{u}_1+C{h}^2})u_1.$$

Full text

††journal: Journal of Computational Physicsmytitlenotemytitlenotefootnotetext: This work was partially supported by the Postdoctoral Science Foundation of China (no.2017M620689), the National Science Foundation of China (nos. 11571048 and 11401034), and the CAEP Developing Fund of Science and Technology (no. 2014A0202009).

1 Introduction

Positivity-preserving or monotonicity is a significant requirement for discrete schemes of diffusion equations. A scheme without such a property can lead to non-physical solutions or oscillation. In some applications, such as Lagrangian computation of hydrodynamic problems and reservoir simulations, meshes are usually distorted and the coefficients are heterogeneous and anisotropic. These two factors may cause diffusion schemes more readily to produce non-physical negative solutions. To avoid this problem, considerable attention has been given to the development of numerical schemes which respect the positivity-preserving property.

Owing to local conservation of finite volume(FV) methods, several positivity-preserving FV schemes have been proposed[2, 11, 15, 16, 19, 29, 35]. It is shown in [24] that no linear and monotone nine-point scheme with second-order accuracy exists on any distorted quadrilateral mesh or for any anisotropic diffusion. Therefore, to get a second-order monotone scheme, one must reference nonlinear monotone schemes. In [26], a nonlinear FV scheme is proposed for a highly anisotropic diffusion operator on unstructured triangular meshes. Since its proposal, this approach has been further developed to obtain schemes that preserve positivity[2, 15, 16, 29, 35] and maximum value principle[11, 19] on general grids. In [35], a nonlinear FV scheme with little severe restriction to the collocation points is presented, which has been further improved in [29, 30] and extended to non-equilibrium radiation diffusion problems[31] and advection diffusion problems[33]. In these works, the convex decomposition of co-normals and the positivity-preserving interpolation of auxiliary unknowns are two key components. However, the convex decomposition of co-normals can result in a non-fixed stencil. Moreover, the positivity-preserving interpolation of auxiliary unknowns is hard to guarantee on general distorted meshes or for anisotropic diffusion problems if we strive to obtain a scheme with accuracy higher than the first order.

In recent years, improved FV schemes have been further developed. In [30], a nonlinear FV scheme with second-order accuracy is proposed, which avoids the assumption that values of auxiliary unknowns are nonnegative. In [17], the authors present an interpolation-free approach based on local repositioning of cell centers, though it is applicable for cells with only one discontinuous face. In [6], a nonlinear scheme with two sets of primary unknowns is constructed, which is interpolation-free. Nevertheless, this scheme can not achieve the second-order accuracy in the case of discontinuous diffusion coefficients. Nonlinear two-point flux approximations that differ from the above approaches are further developed in [1, 10, 38], where neither the positive interpolation of auxiliary unknowns nor the convex decomposition of co-normals is required.

The finite volume element(FVE) method[22], also called the generalized difference method[20], has the properties of simple calculation and local conservation. Compared with FV methods, an advantage of FVE methods is the theoretical basis. Recently, many studies have been devoted to developing higher-order FVE schemes and establishing their error estimations. The readers can refer to the papers[7, 12, 13, 39] and the references therein. However, almost all of their construction and corresponding theoretical analysis are executed for problems with scalar or even constant diffusion coefficients, and they focus on the proofs of coercivity and optimal convergence rates. The development of FVE schemes satisfying positivity preservation is far more difficult to conduct, especially for problems with anisotropic diffusion. This difficulty has impeded the application of FVE methods on some complicated problems. Although there exist some works devoted to the application of FVE methods on complicated numerical simulations[14, 37], these schemes are not monotone and the non-physical solutions are just modified by a simple cutoff method[21] or repair techniques[23]. To our best knowledge, few monotone FVE schemes have been found until now. Therefore, the primary motivation of the present study is to construct a nonlinear monotone FVE scheme with approximate second-order accuracy, which can be applied to real practical problems.

In this paper, we construct a nonlinear monotone FVE scheme in the framework of the iso-parametric bilinear FVE method on quadrilateral meshes. The traditional barycenter dual partition is replaced by an overlapping one, which make it possible to construct a positivity-preserving FVE scheme. Firstly, we employ the iso-parametric bilinear element to approximate the gradients on the barycenters, which is used to get a continuous one-sided flux on each edge of dual elements. Finally, a new two-point flux technique is applied to obtain a monotone scheme with vertex-centered unknowns. In order to make the resulted scheme regular with respect to the unknowns, two parameters are introduced which can still keep the scheme’s second-order accuracy. Accordingly, we prove the existence of a solution for this nonlinear system by the Brouwer’s theorem. It should be pointed out that the nonlinear two-point flux technique employed in this paper differs from those considered in [1, 10, 38] and is also applicable to other second-order schemes. Compared with the standard FVE scheme[14, 28], the proposed monotone FVE scheme has comparable accuracy and can significantly reduce the computational costs. Besides, the construction and implementation of the new scheme are simpler than FV schemes, because the auxiliary unknowns interpolation and the decomposition of co-normals are not needed.

The remainder of this paper is organized as follows. We briefly introduce the standard FVE method on quadrilateral meshes in Section 2. In Section 3, we describe the construction of the monotone FVE scheme and prove the existence of a solution for this new scheme. In Section 4, numerical results are presented to evaluate the performance of the proposed FVE scheme. Finally, conclusions are provided in Section 5.

2 Preliminaries of the standard FVE method

Consider the unsteady diffusion problem:

[TABLE]

where

(a)

$u=u(x,t)$ is the solution of the diffusion equation;

(b)

$\Omega$ is an open-bounded, convex-connected polygonal domain in $\mathbb{R}^{2}$ with boundary $\partial\Omega$, and $T$ is a positive constant;

(c)

The source term $f\in L^{2}(\Omega\times(0,T))$ and the boundary data $g\in H^{1/2}(\partial\Omega\times(0,T))$;

(d)

$\kappa$ is a symmetric tensor such that $\kappa$ is piecewise continuous on $\Omega\times(0,T)$. It may be discontinuous on some interfaces, and the set of eigenvalues of $\kappa$ is included in $[\lambda_{\text{min}},\lambda_{\text{max}}]$ with $\lambda_{\text{min}}>0$.

(e)

$\gamma$ and $\delta$ are smooth functions such that

[TABLE]

The initial data is given by

[TABLE]

Divide $\overline{\Omega}$ into a set of strictly convex quadrilaterals such that different quadrilaterals have no common interior point. Accordingly, the vertexes of any quadrilateral do not lie on the interior of an edge of any other quadrilateral, and any vertex on the boundary is a vertex of some quadrilateral. Thus, we obtain a primary partition $\mathcal{T}_{h}$, where $h$ is the largest diameter of all quadrilaterals. For any element $K\in\mathcal{T}_{h}$ with vertexes $P_{i},i=1,2,3,4,$ as shown in Figure 1, let $M_{i}$ be the midpoint of edge $P_{i}P_{i+1}(P_{5}=P_{1})$ and $Q$ be the barycenter of $K$. Let $\Omega_{h}$ be the set of nodes of the mesh. The number of all nodes is denoted by $\mathcal{N}$.

Let $h_{K}$ denote the diameter of the element $K\in\mathcal{T}_{h}$, $\rho_{K}$ denote the maximum diameter of circles contained in $K$, and $S_{K}$ denote the area of K. Suppose $\mathcal{T}_{h}$ is regular, i.e., there exist two positive constants $C$ and $C_{1}$ such that for any element $K\in\mathcal{T}_{h}$

[TABLE]

To define the trial function space $U_{h}$, we take the unit square $\widehat{K}=[0,1]\times[0,1]$ on the $(\xi,\eta)$ plane as a reference element. For any convex quadrilateral $K$, there exists a unique invertible bilinear mapping $F_{K}$, which maps $\widehat{K}$ onto $K$.

[TABLE]

where

[TABLE]

and $(x_{i},y_{i}),i=1,2,3,4,$ are the coordinates of vertices $P_{i}(i=1,2,3,4)$ on element $K$. Denote the Jacobi matrix of the mapping $F_{K}$ by $J_{K}(\xi,\eta)$, then

[TABLE]

Let $\mathcal{J}_{Q}$ be the determinant of the Jacobi matrix with respect to $Q$. Simple calculation shows that

[TABLE]

Due to the regularity of $\mathcal{T}_{h}$, there exists a positive constant $C$ such that

[TABLE]

Denote by $U_{h}(\widehat{K})$ the iso-parametric bilinear polynomial space on $\widehat{K}$, then the trial function space $U_{h}$ is defined as

[TABLE]

where $\widehat{u}_{h}\circ F_{K}^{-1}$ denotes the compound function of $\widehat{u}_{h}$ and the inverse mapping $F_{K}^{-1}.$

The control volume of vertex $P_{i}$ in $K$, as shown in Figure 1, is the subregion $M_{i-1}P_{i}M_{i}Q,i=1,2,3,4$, where $M_{0}=M_{4}.$ Hence, the control volume associated with vertex $P$ is denoted by $K_{P}^{\ast}$, which is the union of the above subregions sharing the vertex $P$. The set of all control volumes forms a dual partition $\mathcal{T}_{h}^{\ast}$.

Define the test function space as

[TABLE]

The standard FVE scheme for (1)-(2) on quadrilateral meshes is to find $u_{h}(x,\cdot)\in U_{h}$, such that

[TABLE]

with

[TABLE]

[TABLE]

where $P\in\Omega_{h}$ and $\mathbf{n}$ is the outward unit normal vector of $\partial K_{P}^{\ast}$.

In [20], detailed proofs are provided for the coercivity and error estimations of the FVE scheme (10). However, a study on the monotonicity of FVE methods is difficult to conduct and few works can be found. Owing to the nonlinearity of the iso-parametric bilinear element gradient, constraints on computational meshes and material properties prevent FVE methods from practical applications. The numerical results presented in [37] have shown that the standard FVE schemes can lead to negative solutions which are non-physical and make the computation terminate. To overcome this drawback, authors in [37] employ two nonnegative modifications, the cutoff method in [21] and repair techniques in [23], to handle this problem. However, the cutoff method suffers from the conservation error, and the local repair technique increases the computational cost and breaks the governing equations. Therefore, it is urgent to construct a positivity-preserving FVE scheme.

3 Monotone FVE scheme

In this section, we construct a monotone FVE scheme for the spacial discretization. For simplicity and without confusion, we omit the variate $t$ in the following description.

To overcame the difficulties of standard FVE method, we define a new dual partition formed by the diagonals of all quadrilateral elements. For a node $P$, as shown in Figure 2, let the shaded area be the control volume associated with $P$, which is denoted by $K_{P}^{\ast}.$ The set of all such control volumes forms a new dual partition, which is also denoted by $\mathcal{T}_{h}^{\ast}$. Obviously, this new dual partition is overlapping.

We still denote by $V_{h}$ the test function space associated with this new dual partition. The definition of $V_{h}$ is similar to (9), except its support set.

Applying $V_{h}$ to test problem (1)-(2), we obtain the corresponding weak formulation by Green’s theorem:

[TABLE]

with

[TABLE]

where $P\in\Omega_{h}$ and $\mathbf{n}$ is the outward unit normal vector of $\sigma$.

We present a discrete flux $F_{P,\sigma}$ to approximate $\mathcal{F}_{P,\sigma}$. Unlike any FV method, we do not need to decompose the co-normal $\kappa^{T}\mathbf{n}$. We directly discretize $\nabla u$ by the iso-parametric bilinear element.

We take nodes $P_{1},P_{3}$ and an edge $\sigma=P_{2}P_{4}$ (as shown in Figure 3) to present the discretization of continuous fluxes $\mathcal{F}_{P_{1},\sigma},\mathcal{F}_{P_{3},\sigma}$. In Figure 3, $\mathbf{q}_{i},i=1,2,$ denote the normal vectors with respect to the diagonals. Assume that these vectors have the same lengths as their corresponding diagonals. For a given $u_{h}\in U_{h}$, let $u_{i}$ be the value of $u_{h}$ at node $P_{i},i=1,2,3,4.$

On an element as shown in Figure 3, the trial function $u_{h}$ takes the following form

[TABLE]

where

[TABLE]

By simple calculation, the gradient of $u_{h}$ at $Q$ is formulated as follows

[TABLE]

where $\mathcal{J}_{Q}$ is the determinant of the Jacobi matrix corresponding to $Q$.

3.1 Construction of the monotone FVE scheme

In this subsection, we describe the construction of the new monotone FVE scheme, which is a nonlinear two-point flux approximation.

Substitute $\nabla u_{h}\mid_{Q}$ given in (15) into (14), we obtain a one-sided flux $F_{P_{1}}$ with respect to node $P_{1}$

[TABLE]

where

[TABLE]

and $\kappa_{Q}^{T}=\kappa(Q)^{T}$. Firstly, we have the following estimation for $A_{\sigma}$

[TABLE]

where the positive definiteness of $\kappa_{Q}$ and the regularity (7) have been used. Hence, the coefficient $A_{\sigma}$ is bounded. With respect to the node $P_{3}$, continuing to use $\nabla u_{h}\mid_{Q}$ to approximate $\nabla u$, we obtain a one-sided flux $F_{P_{3}}$

[TABLE]

Let

[TABLE]

be the positive and negative parts of $r_{\sigma}(u_{2},u_{4})$, then the one-sided fluxes $F_{P_{1}}$ and $F_{P_{3}}$ can be further rewritten as

[TABLE]

Assume that $u_{h}$ is positive at all vertexes and let $M,C$ be two positive constants independent of $u_{h}$. The value of $M$ will be discussed in next subsection and the value of $C$ is less than or equal to 1 in our numerical tests.

Based on $F_{P_{1}},F_{P_{3}}$ given by (18) and (19), we obtain the two-point flux approximations $F_{P_{1},\sigma},F_{P_{3},\sigma}$ by the following nonlinear technique:

For the node $P_{1}$, we have

[TABLE]

• 2.

For the node $P_{3}$, we have

[TABLE]

Obviously, $F_{P_{1},\sigma},F_{P_{3},\sigma}$ still satisfy the local conservation because $F_{P_{1},\sigma}+F_{P_{3},\sigma}=0$.

Specially, if $P_{1}$ is a boundary node and $\sigma\in\partial K_{P_{1}}^{\ast}$ is a boundary edge, then it holds that from the boundary condition (2)

[TABLE]

We employ the following discretization $F_{P_{1},\sigma}$ to approximate $\mathcal{F}_{P_{1},\sigma}$

[TABLE]

where $u_{P_{1}}$ is a unknown, $\mid\sigma\mid$ is the length of edge $\sigma$, and $\delta_{\sigma},\gamma_{\sigma}$ and $g_{\sigma}$ are the approximations of $\delta,\gamma$ and $g$ on the edge $\sigma$ respectively. There are several ways to calculate $\delta_{\sigma},\gamma_{\sigma}$ and $g_{\sigma}$, such as the midpoint value or area mean on the edge $\sigma$.

Finally, the monotone FVE scheme can be formulated as: Find $u_{h}(x,\cdot)\in U_{h}$, such that

[TABLE]

Obviously, this monotone FVE scheme is nonlinear and there are five nonzero elements in each row. If the diffusion coefficient $\kappa$ is linear, the solution of this scheme is less time-consumming than linear schemes. Nevertheless, numerical tests show that, if the diffusion coefficient $\kappa$ is nonlinear, the monotone FVE scheme (22) needs fewer linear iterations per time step and reduce the costs by almost half compared with standard FVE method.

Remark 3.1**.**

*Although the monotone FVE scheme (22) is constructed for unsteady diffusion equations, it is likewise applicable for steady diffusion equations. For convenience, we take several steady diffusion problems to examine the monotonicity and convergence of this new scheme in our numerical tests. *

Remark 3.2**.**

*The monotone FVE scheme (22) is valid for the case that the diffusion coefficients are discontinuous, provided that the cell sides fit the interfaces. *

3.2 Analysis of the truncation error

In the nonlinear technique (20)-(21), the two parameters $M$ and $C$ are introduced to ensure the scheme’s regularity with respect to the unknowns. However, we will show that they do not break the second-order accuracy in numerical experiments. Next, we give the analysis of the truncation error for this new monotone FVE scheme.

Under the assumption that $u_{h}$ is positive at all vertexes, there exists a positive constant $M$ which is large enough such that

[TABLE]

Similar to the standard FVE method, under the assumptions (3)-(4), the one-sided flux $(\ref{eqF24})$ satisfies

[TABLE]

which is a classical estimation of the finite element space $U_{h}$ [20]. For the nonlinear two-point flux approximation $F_{P_{1},\sigma}$ given by $(\ref{eqF27})$, we have

[TABLE]

Due to (23) and $r_{\sigma}(u_{2},u_{4})=O(h)$, we further get that

[TABLE]

Therefore, the truncation error of the monotone FVE scheme (22) remains $O(h^{2})$, which is a necessary condition for a scheme to guarantee the second-order accuracy.

As analysed in [3], the resulted flux is not consistent on non-rectangular meshes, if employing a difference quotient to approximate the gradient directly. Hence, usually a combination of two one-sided fluxes is employed to obtain a second-order monotone FV scheme. Nevertheless, the iso-parametric bilinear finite element space has a first-order approximate capability to the gradient. Then, the associated discrete flux presented in this paper still keeps a second-order accuracy even if the combination is omitted.

3.3 Monotonicity and nonlinear iteration

Let $\triangle t$ be the time size, $U^{n}=(U^{n}_{P})_{P\in\Omega_{h}}\in\mathrm{R}^{\mathcal{N}}$ be the vector of discrete unknowns at the $n-$th time level. Denote the area of dual element $K_{P}^{\ast}$ by $S_{K_{P}^{\ast}}$. Let $L$ be the diagonal matrix of the areas of dual elements.

Taking the basis functions of $V_{h}$ as test functions in (22) and using an implicit discretization in time, we obtain the full discretization for system (1)-(2):

[TABLE]

Here, $A(U^{n+1})$ is a matrix corresponding to the spacial discretization, and $f^{n+1}=(f_{P}^{n+1})_{P\in\Omega_{h}}$ and $g^{n+1}=(g_{P}^{n+1})_{P\in\Omega_{h}}$ are two vectors corresponding to the source term and boundary condition. Finally, the equation (24) can be further rewritten as

[TABLE]

For any vector $U\in\mathrm{R}^{\mathcal{N}}$ and $U\geq 0$, the matrix $A(U)$ has the following properties:

All diagonal entries of matrix $A(U)$ are positive; 2. 2.

All off-diagonal entries of $A(U)$ are non-positive; 3. 3.

Column sum corresponding to the interior node is 0 and column sum corresponding to the boundary node is positive.

There exist many nonlinear iteration methods to solve system (25). Since the efficiency of nonlinear iteration methods is not the focused issue in this paper, we just employ a simple nonlinear iteration, the relaxed version of the Picard iteration method, to solve the nonlinear system (25). This relaxed version demonstrates more robust behavior [18] and the nonlinear iterations are all convergent in our numerical tests. We choose a small value $\varepsilon_{non}>0$ and initial vector $\widehat{U}^{0}\geq 0$, and repeat for $k=1,2,\cdots$,

Solve $(L+\triangle tA(\widehat{U}^{k-1}))\widetilde{U}^{k}=\triangle t(f^{n+1}+g^{n+1})+LU^{n}$; 2. 2.

$\widehat{U}^{k}=\widehat{U}^{k-1}+\omega(\widetilde{U}^{k}-\widehat{U}^{k-1})$, where $0<\omega\leq 1$ is the damping factor; 3. 3.

Stop if $\parallel\widehat{U}^{k}-\widehat{U}^{k-1}\parallel\leq\varepsilon_{non}\parallel\widehat{U}^{0}\parallel$; 4. 4.

The value of $U^{n+1}$ is then given by the last $\widehat{U}^{k}$.

According to the definition of M-matrix [34], it is noted that $A(U)^{T}$ is an M-matrix for any vector $U\geq 0$ and hence $(L+\triangle tA(U))^{T}$ is also an M-matrix. It holds from the properties of M-matrix [34] that the matrix $(L+\triangle tA(U))^{-T}$ has non-negative elements. Under the assumptions that $(f^{n+1}+g^{n+1})\geq 0,U^{0}\geq 0$ and the linear systems in the above Picard iterations are exactly solved, it is readily evident that all iterative solutions are non-negative vectors.

3.4 Existence of a solution for the nonlinear FVE scheme

It is easy to observe that the monotone FVE scheme (22) is regular with respect to $u_{h}$ for any fixed positive constant $M$. Accordingly, the existence of a solution for nonlinear system (25) could be obtained by the Brouwer’s theorem.

Theorem 3.1**.**

*If $f^{n+1}+g^{n+1}\geq 0$ and $U^{n}\geq 0$, then system (25) has a solution $U^{n+1}$ for any fixed positive constant $M.$ *

Proof. 1**.**

Let $C_{0}$ be a non-negative constant. Define a compact set

[TABLE]

and a mapping $\phi:\mathcal{C}\mapsto\mathrm{R}^{\mathcal{N}}$ such that

[TABLE]

We need to prove that $\phi$ has a fixed point in order to prove that equation (25) has a solution.

In view of the discussion in the above subsection, we prove that the matrix $(L+\triangle tA(w))^{-1}$ has non-negative elements. Since $\triangle t(f^{n+1}+g^{n+1})+LU^{n}\geq 0$, we have

[TABLE]

Rewrite (26) as

[TABLE]

Multiplying (27) by the constant vector $(1,\ldots,1)$ on the left, we deduce that

[TABLE]

where we have used the facts that $\phi(w)\geq 0$ and all column sums in $A(w)$ are no-negative.

The righthand side of (28) is a non-negative constant. We take $C_{0}=\sum_{P\in\Omega_{h}}(\triangle t(f_{P}^{n+1}+g_{P}^{n+1})+S_{K_{P}^{\ast}}U^{n}_{P})$. Then, $\phi$ maps $\mathcal{C}$ into itself.

*The set $\mathcal{C}$ is a convex compact subset of $\mathrm{R}^{\mathcal{N}}$. Since each coefficient of $A(w)$ is a continuous function of $w$ and $A(w)\mapsto(L+\triangle tA(w))^{-1}$ is continuous from the set M-matrices to the set of matrices, $\phi$ is continuous. Hence, we may apply Brouwer’s theorem, which implies that $\phi$ has a fixed point in $\mathcal{C}$, i.e., the system (25) has a solution. The theorem 3.1 is proved. *

4 Numerical experiments

In this section, we examine numerical performance of the monotone FVE method proposed in this paper. The monotonicity and convergence study of the discrete solution are presented for steady diffusion problems on both uniform and random distorted quadrilateral meshes. Moreover, to illustrate the efficiency of our scheme, we apply it to solve an equilibrium radiation diffusion equation. In all examples, we compare this new monotone FVE scheme with the standard FVE method[14, 20, 28].

In the nonlinear iterations, the prescribed tolerance $\varepsilon_{non}$ is $10^{-7}$, and we use the GMRES method[4] to solve the linear systems. In our numerical experiments, we use the following discrete $L_{2}$-norm and discrete $H^{1}$-semi-norm to evaluate the approximation errors:

[TABLE]

where $S_{K^{\ast}_{P}}$ and $S_{K}$ are the areas of $K^{\ast}_{P}$ and $K$, respectively. The rate of convergence is obtained by the following formula:

[TABLE]

where $h_{1},h_{2}$ denote the mesh sizes of two successive meshes, and $E(h_{1}),E(h_{2})$ are the corresponding errors.

The numerical tests are performed on uniform meshes and random distorted quadrilateral meshes, as shown in Figure 4. The distorted meshes are constructed from the uniform meshes by adding a random perturbation to the interior nodes. Unlike the convergence study for the standard FVE methods in [22], the distortion in this paper is performed on each refinement level.

4.1 Monotonicity

The computational domain is a bi-unit square. The forcing function is taken as

[TABLE]

The diffusion tensor is given by

[TABLE]

Consider a steady problem (1)-(2) with boundary condition $10^{-9}(\kappa\nabla u)\cdot\mathbf{n}+u=0.$ In this paper, we take $\alpha=0.01$.

We test the proposed monotone FVE scheme on uniform meshes and random quadrilateral meshes (see Figure 4). The exact solution $u(x,y)$ is unknown; nevertheless, the maximum principle indicates that it is non-negative. The numerical solutions obtained by the standard FVE scheme are shown in Figure 5, and the counterparts for the monotone FVE scheme are shown in Figures 6. In these figures, we observe that the standard FVE scheme produces negative values; however, the monotone FVE scheme preserves the positivity of the continuous solution.

4.2 Convergence

$\mathbf{Example~{}1}$: Convergence study for the isotropic diffusion problem. Let $\Omega=(0,1)^{2}$, the exact solution and diffusion coefficient be as follows:

[TABLE]

The source term $f$ and the boundary data $(\kappa\nabla u)\cdot\mathbf{n}+u=g$ are set accordingly to the exact solution. Since the positivity of the discrete unknowns is a requirement for the proposed monotone FVE scheme, the exact solution, and $f+g$ are non-negative in this example. Numerical errors and convergence rates for standard FVE method and monotone FVE scheme are presented in Table 1 and Table 2, respectively.

$\mathbf{Example~{}2}$: Convergence study for the anisotropic diffusion problem. Let $\Omega=(0,1)^{2}$, the exact solution and diffusion coefficient be as follows:

[TABLE]

[TABLE]

where $k_{1}=0.1,k_{2}=4,\theta=\pi/6$. The source term $f$ and the boundary data $g=(\kappa\nabla u)\cdot\mathbf{n}+1000u$ are set accordingly to the exact solution. Likewise, the exact solution, and $f+g$ are non-negative in this example. Numerical errors and convergence rates for standard FVE method and monotone FVE scheme are presented in Table 3 and Table 4, respectively.

$\mathbf{Example~{}3}$: Convergence study for the diffusion problem with a discontinuous coefficient. Let $\Omega=(0,1)^{2}$ and the exact solution and diffusion coefficient be as follows:

[TABLE]

[TABLE]

The source term $f$ and the boundary data $g=(\kappa\nabla u)\cdot\mathbf{n}+10000u$ are set accordingly to the exact solution. We have examined that the exact solution, and $f+g$ are non-negative. Numerical errors and convergence rates for standard FVE method and monotone FVE scheme are presented in Table 5 and Table 6, respectively.

From the numerical results presented in Tables 1-6, it is evident that the proposed monotone FVE scheme can achieve nearly the same convergence rates as the standard FVE method in both the discrete $L_{2}$-norm and discrete $H^{1}$-semi-norm, and the numerical errors are comparable for these two methods. On both uniform meshes and random distorted quadrilateral meshes, the proposed method have an approximate second-order convergence rate in the discrete $L_{2}$-norm and a higher than first-order convergence rate in the discrete $H^{1}$-semi-norm. Moreover, super-convergence in the discrete $H^{1}$-semi-norm can be still observed on uniform meshes.

4.3 Equilibrium radiation diffusion problem

To examine the computational efficiency of the new monotone FVE scheme, we apply it to the equilibrium radiation diffusion problem. We provide a brief description of this equation. See [27, 36] for a more detailed discussion. Consider the following nonlinear parabolic partial differential equation:

[TABLE]

where $\Omega=(0,1)^{2}$, $E$ is the dimensionless gray radiation energy density, and $D_{L}(E)$ is the Larsen s form of the flux-limited diffusion coefficient defined as

[TABLE]

Here, $D(E)=Z^{-3}E^{\frac{3}{4}}$, and $Z$ is the atomic number of the medium, which can be a discontinuous function. The value for $Z$ is 1 everywhere, except in the two obstacles defined by

[TABLE]

and

[TABLE]

where the value for $Z$ is 10. The boundary condition is

[TABLE]

The initial radiation energy is given by

[TABLE]

where $r=\sqrt{x^{2}+y^{2}}.$

For the time discretization, we use the backward Euler method. We take the time step size as $5.0e-4$ and let the final state be $T=1.0$.

Solution profiles for standard FVE method and monotone FVE scheme are presented in Figures 7 and 8, respectively. In comparing Figures 7 and 8, we observe that the solution obtained by these two methods are nearly the same. Moreover, the variation tendencies of the solutions on random quadrilateral meshes are in accordance with those on rectangular meshes for monotone FVE scheme. In Table 7, we give the $L_{2}$-norm of the solution on rectangular and random quadrilateral meshes. From this table, it is apparent that the $L_{2}$-norm of the solution for each scheme on random quadrilateral meshes is close to that on rectangular meshes.

Next, we give the comparison of computational costs between our scheme and the standard FVE scheme. Table 8 provides the average number of nonlinear iterations per time step, the average number of linear iterations per nonlinear iteration and total CPU time on random quadrilateral meshes. There are nine nonzero elements in each row for the linear system of the standard FVE method; however, there are only five nonzero elements in each row for the linear system of our scheme. Hence, our scheme reduces the costs by almost half compared with the standard FVE method when the total number of linear iterations is the same. Table 8 shows that the proposed monotone FVE scheme requires fewer linear iterations in each time step than the standard FVE method and the total CPU time is reduced significantly. Hence, the proposed monotone FVE scheme is more efficient than the standard FVE method.

5 Conclusion

In this paper, we construct a nonlinear positivity-preserving FVE scheme on quadrilateral meshes, where the one-sided flux is approximated by the standard iso-parametric bilinear element. The new scheme is positivity-preserving by applying a nonlinear two-point flux technique and is applicable for both continuous and discontinuous anisotropic problems. Compared with the FV schemes, the main feature of this monotone FVE scheme is that the decomposition of co-normals and the introduction of auxiliary unknowns are both avoided. Numerical results show that the proposed monotone FVE scheme do not produce non-physical solutions and can reach nearly optimal convergence rates in both the $L_{2}$-norm and $H^{1}$-semi-norm. From the performance in solving the equilibrium radiation diffusion equation, it can be observed that the presented scheme can significantly reduce computational costs compared with the standard FVE method.