# Nonlinear Finite Volume Scheme Preserving Positivity for 2D Convection-Diffusion Equations on Polygonal Meshes

Authors: Bin Lan; Jianqiang Dong
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, a nonlinear finite volume scheme preserving positivity for solving 2D steady convection-diffusion equation on arbitrary convex polygonal meshes is proposed. First, the nonlinear positivity-preserving finite volume scheme is developed. Then, in order to avoid the computed solution beyond the upper bound, the cell-centered unknowns and auxiliary unknowns on the cell-edge are corrected. We prove that the present scheme can avoid the numerical solution beyond the upper bound. Our scheme is locally conservative and has only cell-centered unknowns. Numerical results show that our scheme preserves the above conclusion and has second-order accuracy for solution.


---

## Body

## 1. Introduction

Convection-diffusion equations are widely used in the fields of solid mechanics, material science, image processing, and so on. So, it is both theoretically and practically important to investigate numerical methods for such equations. An accurate numerical method must maintain the fundamental properties of practical problems. The extremum principle is an important property of solutions for the convection-diffusion equation. It includes minimum principle and maximum principle. The authors of [1, 2] pointed out that the discrete maximum principle (DMP) plays an important role in proving the existence and uniqueness of discrete solution, enforcing numerical stability, and deriving convergence for a sequence of approximate solutions [3]. Pert [4] pointed out that a scheme violating extremum principle can lead to two problems: fully implicit discretization with large time-steps has relatively poor accuracy, and spurious negative values are generated. Moreover, it is proved that a linear operator, resulting from the discretization of diffusion equations, satisfies extremum principle if and only if it is both differential and nonnegativity maintaining.In general cases, the discrete extremum principle (DEP) is more restrictive than monotonicity (positivity-preserving). However, it is difficult to construct a reliable discretization method that satisfies the DEP on arbitrary convex polygonal meshes. Hence, positivity-preserving is one of the key requirements to discrete schemes for the convection-diffusion equation, which says that it can only guarantee nonnegative bound of the numerical solution. Sheng and Yuan [2] pointed out that the scheme without positivity-preserving can lead to the violation of the entropy constraints of the second law of thermodynamics, causing heat to flow from regions of lower temperature to higher temperature. In regions of large temperature variations, this can cause the temperature to become negative.The finite volume methods (FVM) guarantee the local conservation. But many classical schemes fail to maintain positivity for strong anisotropic diffusion tensors or on distorted meshes [5–8]. Some nonlinear methods have been developed [9–21] for general diffusion or convection-diffusion equations, which guarantee the positivity on general or distorted meshes for general tensor coefficients.Bertolazzi and Manzini [22] proposed a MUSCL-like cell-centered finite volume method, where the discretization of advective fluxes is based on a least-square reconstruction of the vertex values from cell averages. Lipnikov et al. [23] proposed a new slope limiting technique based on a specially minimal nonlinear correction, which follows the ideas of the monotonic upstream-centered scheme for conservation laws (MUSCL). Then, in the studies by Wang et al. and Zhang et al. [16, 24], the limiting technique is used to avoid nonphysical oscillation. In the study by Lan et al. [21], a new upwind scheme is used to discretize the convective flux, and the method did not introduce any slope limiting technique.In this paper, we develop a nonlinear FV scheme, which satisfies DEP for convection-diffusion problems on arbitrary convex polygonal meshes. Following the idea of the discretization for diffusive flux [18] and convective flux [21], the adaptive approach of choosing stencil is applied. Positivity-preserving scheme can only guarantee nonnegative bound of the numerical solution. Considering that the computation of value on the cell edge and the value of cell-centered unknowns may be out of bound, a correcting technique is introduced. Our scheme is constructed by a nonlinear combination technique and has second-order accuracy for the solution and first-order for the flux.The article is organized as follows. The model problem is described, and some notations are introduced in Section2. The main process of construction for the 2D steady convection-diffusion equation is given in Section 3. In Section 4, several numerical tests are exhibited to illustrate the features of our scheme. At last, some conclusions are given in Section 5.


## 2. The Problem and Notation

Consider the following stationary convection-diffusion problem for unknown functionu=ux:(1)−∇⋅κ∇u−v⟶u=f,inΩ,(2)u=g,on∂Ω,where Ω is a bounded polygonal domain in ℝ×ℝ with boundary ∂Ω, κ=κx is a known diffusive coefficient, and v⟶=v⟶x is a velocity vector field.Assume that the functionsv⟶x,fx, and gx satisfy the constraints listed as follows:(3)∇⋅v⟶≥0,v⟶∈C1Ω¯2,f∈L2Ω,g∈H1/2∂Ω∩C∂Ω,and there are two positive constants λ1 and λ2 such that(4)λ1ξ2≤κxξ⋅ξ≤λ2ξ2,∀ξ∈ℝ×ℝ.The solvability of the problem (1)-(2) has been given, and the maximum and minimum principle can be found in the study by Gilbarg and Trudinger [25].We use a mesh onΩ made up of arbitrary convex polygon cells. The set of all cells, edges, and nodes are denoted by T,ε, and N, respectively.We denote the cell byK,L,…, and the cell center is also denoted by K,L,…,. In addition, the common edge of two cells K and L is denoted by σ, i.e., σ=KL∈Ε. The cell-edge σ is also denoted by AB, and the midpoint of σ is denoted by M. Moreover, we denote P1 and P2 are two adjacent midpoints of cell K (Figure 1).Figure 1
Local stencil 1.Letn⟶K,σ (or n⟶L,σ) be the unit outer normal vector on the cell-edge σ of cell K (or L), κT be the transpose of matrix κ, and ℰK be the set of all edges of cell K. Denote εint=ε∩Ω and εext=ε∩∂Ω. Denote h=supK∈TmK1/2, where mK is the area of cell K.Integrating (1) over the cell K, we obtain(5)∑σ∈ΕKℱK,σ+GK,σ=∫Kfxdx,where the diffusive and convective flux are defined as(6)ℱK,σ=−∫σ∇ux⋅κTxn⟶K,σdl,(7)GK,σ=∫σv⟶ux⋅n⟶K,σdl.


## 3. Construction of the Scheme

### 3.1. The Diffusive Flux

Following the idea in the study by Sheng and Yuan [18], the adaptive approach of choosing stencil is applied for the approximation of the diffusive flux (equation (6)) together with a nonlinear combination technique. In the method, the two nonnegative parameters are introduced to define a nonlinear two-point flux. Then, the continuity of normal flux on the cell edge is used to give the final discretization of the diffusive flux. At last, the continuity of normal flux is used to obtain the value of uM. First, we give a brief review of the construction.Figure1 shows that a ray originating at the point K along the direction κTn⟶K,σ must intersect one segment connecting two neighboring midpoints of edge of cell K, where the two midpoints are denoted by P1 and P2, and the cross point is denoted by O1. Similarly, a ray originating at the midpoint M along the direction −κTn⟶K,σ must intersect one certain KP4, where P4 must be one vertex of σ, and the cross point is denoted by O2.Lett⟶KP1,t⟶KP2,t⟶MK, and t⟶MP4 be some unit tangential vectors along their corresponding directions, respectively. θii=1,…,4 are some corresponding angles. Hence, we established the following relations:(8)κTn⟶K,σκTn⟶K,σ=sinθ2sinθ1+θ2t⟶KP1+sinθ1sinθ1+θ2t⟶KP2,(9)−κTn⟶K,σκTn⟶K,σ=sinθ4sinθ3+θ4t⟶MK+sinθ3sinθ4+θ4t⟶MP4.Substituting equation (8) into equation (6) and neglecting the high-order terms, we have(10)F¯1=−κKTn⟶K,σσsinθ2sinθ1+θ2uP1−uKKP1+sinθ1sinθ1+θ2uP2−uKKP2=a1uK−uP1+a2uK−uP2,where(11)a1=κKTn⟶K,σσKP1sinθ2sinθ1+θ2,a2=κKTn⟶K,σσKP1sinθ1sinθ1+θ2.Similarly, substituting equation (9) into equation (6), we have(12)F¯2=κKTn⟶K,σσsinθ4sinθ3+θ4uK−uMMK+sinθ3sinθ3+θ4uP4−uMMP4=a3uK−uM+a4uP4−uM,where(13)a3=κKTn⟶K,σσMKsinθ4sinθ3+θ4,a4=κKTn⟶K,σσMP4sinθ3sinθ3+θ4.Combined with equations (10)–(12), the discrete normal flux on σ can be defined as follows:(14)FK,σ=μ1F¯1+μ2F¯2=μ1a1uK−uP1+μ1a2uK−uP2+μ2a3uK−uM+μ2a4uP4−uM,where μ1 and μ2 are some nonlinear coefficients with μ1+μ2=1, which will be given later.In order to assureμ1 and μ2 are positive, two additional parameters ω1 and ω2 are introduced later. According to the different positions of P1 and P2, three cases exist.We assume thatuP1≠uM and uP2≠uM, the normal flux (14) can be rewritten as(15)FK,σ=μ1a1+a2+μ2a3uK−μ2a3+a4uM−μ1a1uP1+a2uP2+μ2a4uP4=μ1a1+a21+ω1+μ2a3uK−μ2a3+a41+ω2uM−μ1a1uP1+ω1uK+a2uP2+ω1uK+μ2a4uP4+ω2uM.In order to obtain the two-point flux approximation, the last two terms of the above expression should vanish; hence,μ1 and μ2 are given as follows:(16)μ1=a4uP4+ω2uMa1uP1+ω1uK+a2uP2+ω1uK+a4uP4+ω2uM,μ2=a1uP1+ω1uK+a2uP2+ω1uKa1uP1+ω1uK+a2uP2+ω1uK+a4uP4+ω2uM.Hence, (15) can be expressed as follows:(17)FK,σ=AK,σ,1uK−AK,σ,2uM,where(18)AK,σ,1=μ1a1+a21+ω1+μ2a3,AK,σ,2=μ2a3+a41+ω2.In order to assureμ1>0 and μ2>0, two parameters ω1 and ω2 can be chosen such that(19)a1uP1+ω1uK+a2uP2+ω1uK≥0,a4uP4+ω2uM≥0.If(20)a1uP1+ω1uK+a2uP2+ω1uK=a4uP4+ω2uM=0.We letμ1=μ2=1/2.IfuP1=uM or uP2=uM, equation (14) can be expressed in the similar form (17) by using the above method.Similarly, on the edgeσ of the cell L, we have(21)FL,σ=AL,σ,1uL−AL,σ,2uM.Using the continuity of normal fluxFK,σ+FL,σ=0 on edge σ, we obtain(22)uM=AK,σ,1uK+AL,σ,1uLAK,σ,2+AL,σ,2.Substitute equation (22) into equation (17) to obtain the nonlinear two-point diffusive flux on σ=KL:(23)FK,σ=AK,σuK−AL,σuL,where AK,σ=AK,σ,1AL,σ,2/AK,σ,2+AL,σ,2, and AL,σ=AK,σ,2AL,σ,1/AK,σ,2+AL,σ,2.From the computation of vertex unknowns, a method with second-order accuracy has been proposed in the study by Sheng and Yuan [26]. We know that uM>0 as long as uK>0 and uL>0 in equation (22).


### 3.2. The Convective Flux

We focus on the expression of convective flux in equation (7) for ∀σ=KL∈εint and obtain(24)GK,σ=uM∫σv⟶⋅n⟶K,σdl+Oh2=uMvK,σ+−vK,σ−+Oh2,where uM is the value of midpoint M on the cell-edge σ, and(25)vK,σ+=12vK,σ+vK,σ,vK,σ−=12vK,σ−vK,σ,vK,σ=∫σv⋅n⟶K,σdl.In order to ensure that the discretization of equation (24) has the same structure as the scheme (23), we divide the integral term into positive part vK,σ+ and negative part vK,σ−. Moreover, the property of upwind is also considered.Neglecting the high-order terms, we have the following approximate expression of the upwind formula [27]:(26)GK,σ≈uMvK,σ+−vK,σ−.In order to approximate the continuous fluxGK,σ on the cell-edge σ with second-order accuracy, we propose the following method.A local stencil is given in Figure2. For the cell K, M is the midpoint of an arbitrary edge and M1,M2 are the other two midpoints adjacent to it. We denote S△MM1K be the area of triangle MM1K and define(27)u¯K=∑M∈σ,σ∈ℰKωMuM,where ωM=S△MM1K+S△KM2M−1/∑M∈σ,σ∈εKS△MM1K+S△KM2M−1. For a special case (Figure 3), we set the cell K as a triangle and define u¯K=ωMuM+ωM1uM1+ωM2uM2, where(28)ωM=s1s,ωM1=s2s,ωM2=s3s,s=s1+s2+s3,s1=S△MM1K+S△KM2M−1,s2=S△MM1K+S△KM1M2−1,s3=S△KM1M2+S△MKM2−1.Figure 2
Local stencil 2.Figure 3
Local stencil of triangle.It is obvious thatu¯K is a second-order approximation to uK, i.e., uK−∑M∈σ,σ∈ℰKωMuM=Oh2 if the solution u∈C2K.Then, the approximation ofGK,σ on cell-edge σ can be defined as follows.Forσ∈ℰint, we define(29)GK,σ=BK,σuK−BL,σuL,where(30)BK,σ=vK,σ+uMu¯K≥0,BL,σ=vL,σ+uMu¯L≥0.Forσ∈ℰext, we define(31)GK,σ=BK,σuK−bK,σ,where(32)BK,σ=vK,σ+uMu¯K≥0,bK,σ=vK,σ−gM.


### 3.3. The Finite Volume Scheme and Picard Iteration

By using the definition of discretization of diffusive and convective flux, the finite volume scheme can be constructed as follows:(33)∑σ∈εKFK,σ+GK,σ=KFK+∑σ∈εK∩εextAK,σ,2uM+bK,σ,K∈J.(34)uMi=gMi,∀Mi∈∂Ω,where fK=fK and gMi=gxMi.LetU be the discrete unknown vector and AU be the coefficient matrix. A nonlinear algebraic system of the schemes (33) and (34) can be formed: AUU=F. The AU is assembled by the coefficients of diffusive term FK,σ and convective term GK,σ. We use the Picard nonlinear iteration method to solve the system: choose a small value ℰnon>0 and initial vector U0>0 and repeat for nonlinear iteration index s=1,2,…,:(1)
SolveAUs−1Us=F(2)
Stop ifAUsUs−F2≤ΕnonAU0U0−F2The linear algebraic system with coefficient matrixAUs−1 is solved by the biconjugate gradient stabilized (BiCGSTab) method, and the linear iterations are terminated when relative norm of the initial residual becomes smaller than εlin.


### 3.4. The Algorithm

In this subsection, we describe the detailed algorithm.Step 1. InitializeU0>0, εnon, and εlin.Step 2. Whens=0,(a)
computeAK,σ0 and BK,σ0;(b)
compute initial residualAU0U0−F2.Step 3. Whens=1,2,…,(a)
solveAUs−1Us=F;(b)
correctUs, see Remark 1;(c)
computeuPs and correct uMs, ∀M∈σ,σ∈Ε,P∈N, see Remark 2;(d)
computeAK,σs and BK,σs;(e)
compute residualAUsUs−F2;(f)
whetherAUsUs−F2≤εnonAU0U0−F2, if true, then go to (Step 4), otherwise, go to (Step 3).Step 4. Stop.Remark 1.
For∀uKs, if uKs>maxuMs−1,M∈σ,σ∈ΕK, let uKs=u¯Ks−1.Remark 2.
The value of∀uMs can be obtained by equation (22). If ∀uMs>maxuKs,uLs,uAs,uBs, where the common edge σ of two cells K and L is also denoted by AB (Figure 1). Let uMs=uMs−1.
It should be noted that the algorithm in Remark1 is important to avoid the numerical solution beyond the upper bound where the numerical results need to depend on nonnegative initial values of nonlinear iteration.Theorem 1.
LetF≥0,U0≥0, and linear systems in Picard iterations are solved exactly. Then,(35)Us≥0,s=1,2,3,…,.The detailed proof of positivity is given in the study by Yuan and Sheng [11].Now, we state our conclusion, which says that our scheme can avoid the numerical solution beyond the upper bound. Denoteumax=max0,uK,uM,∀K∈T,∀M∈ε.We assumeuK0=umax. Using Remark 1, we know that uK0≤maxuM,M∈σ,σ∈εK0.


## 3.1. The Diffusive Flux

Following the idea in the study by Sheng and Yuan [18], the adaptive approach of choosing stencil is applied for the approximation of the diffusive flux (equation (6)) together with a nonlinear combination technique. In the method, the two nonnegative parameters are introduced to define a nonlinear two-point flux. Then, the continuity of normal flux on the cell edge is used to give the final discretization of the diffusive flux. At last, the continuity of normal flux is used to obtain the value of uM. First, we give a brief review of the construction.Figure1 shows that a ray originating at the point K along the direction κTn⟶K,σ must intersect one segment connecting two neighboring midpoints of edge of cell K, where the two midpoints are denoted by P1 and P2, and the cross point is denoted by O1. Similarly, a ray originating at the midpoint M along the direction −κTn⟶K,σ must intersect one certain KP4, where P4 must be one vertex of σ, and the cross point is denoted by O2.Lett⟶KP1,t⟶KP2,t⟶MK, and t⟶MP4 be some unit tangential vectors along their corresponding directions, respectively. θii=1,…,4 are some corresponding angles. Hence, we established the following relations:(8)κTn⟶K,σκTn⟶K,σ=sinθ2sinθ1+θ2t⟶KP1+sinθ1sinθ1+θ2t⟶KP2,(9)−κTn⟶K,σκTn⟶K,σ=sinθ4sinθ3+θ4t⟶MK+sinθ3sinθ4+θ4t⟶MP4.Substituting equation (8) into equation (6) and neglecting the high-order terms, we have(10)F¯1=−κKTn⟶K,σσsinθ2sinθ1+θ2uP1−uKKP1+sinθ1sinθ1+θ2uP2−uKKP2=a1uK−uP1+a2uK−uP2,where(11)a1=κKTn⟶K,σσKP1sinθ2sinθ1+θ2,a2=κKTn⟶K,σσKP1sinθ1sinθ1+θ2.Similarly, substituting equation (9) into equation (6), we have(12)F¯2=κKTn⟶K,σσsinθ4sinθ3+θ4uK−uMMK+sinθ3sinθ3+θ4uP4−uMMP4=a3uK−uM+a4uP4−uM,where(13)a3=κKTn⟶K,σσMKsinθ4sinθ3+θ4,a4=κKTn⟶K,σσMP4sinθ3sinθ3+θ4.Combined with equations (10)–(12), the discrete normal flux on σ can be defined as follows:(14)FK,σ=μ1F¯1+μ2F¯2=μ1a1uK−uP1+μ1a2uK−uP2+μ2a3uK−uM+μ2a4uP4−uM,where μ1 and μ2 are some nonlinear coefficients with μ1+μ2=1, which will be given later.In order to assureμ1 and μ2 are positive, two additional parameters ω1 and ω2 are introduced later. According to the different positions of P1 and P2, three cases exist.We assume thatuP1≠uM and uP2≠uM, the normal flux (14) can be rewritten as(15)FK,σ=μ1a1+a2+μ2a3uK−μ2a3+a4uM−μ1a1uP1+a2uP2+μ2a4uP4=μ1a1+a21+ω1+μ2a3uK−μ2a3+a41+ω2uM−μ1a1uP1+ω1uK+a2uP2+ω1uK+μ2a4uP4+ω2uM.In order to obtain the two-point flux approximation, the last two terms of the above expression should vanish; hence,μ1 and μ2 are given as follows:(16)μ1=a4uP4+ω2uMa1uP1+ω1uK+a2uP2+ω1uK+a4uP4+ω2uM,μ2=a1uP1+ω1uK+a2uP2+ω1uKa1uP1+ω1uK+a2uP2+ω1uK+a4uP4+ω2uM.Hence, (15) can be expressed as follows:(17)FK,σ=AK,σ,1uK−AK,σ,2uM,where(18)AK,σ,1=μ1a1+a21+ω1+μ2a3,AK,σ,2=μ2a3+a41+ω2.In order to assureμ1>0 and μ2>0, two parameters ω1 and ω2 can be chosen such that(19)a1uP1+ω1uK+a2uP2+ω1uK≥0,a4uP4+ω2uM≥0.If(20)a1uP1+ω1uK+a2uP2+ω1uK=a4uP4+ω2uM=0.We letμ1=μ2=1/2.IfuP1=uM or uP2=uM, equation (14) can be expressed in the similar form (17) by using the above method.Similarly, on the edgeσ of the cell L, we have(21)FL,σ=AL,σ,1uL−AL,σ,2uM.Using the continuity of normal fluxFK,σ+FL,σ=0 on edge σ, we obtain(22)uM=AK,σ,1uK+AL,σ,1uLAK,σ,2+AL,σ,2.Substitute equation (22) into equation (17) to obtain the nonlinear two-point diffusive flux on σ=KL:(23)FK,σ=AK,σuK−AL,σuL,where AK,σ=AK,σ,1AL,σ,2/AK,σ,2+AL,σ,2, and AL,σ=AK,σ,2AL,σ,1/AK,σ,2+AL,σ,2.From the computation of vertex unknowns, a method with second-order accuracy has been proposed in the study by Sheng and Yuan [26]. We know that uM>0 as long as uK>0 and uL>0 in equation (22).


## 3.2. The Convective Flux

We focus on the expression of convective flux in equation (7) for ∀σ=KL∈εint and obtain(24)GK,σ=uM∫σv⟶⋅n⟶K,σdl+Oh2=uMvK,σ+−vK,σ−+Oh2,where uM is the value of midpoint M on the cell-edge σ, and(25)vK,σ+=12vK,σ+vK,σ,vK,σ−=12vK,σ−vK,σ,vK,σ=∫σv⋅n⟶K,σdl.In order to ensure that the discretization of equation (24) has the same structure as the scheme (23), we divide the integral term into positive part vK,σ+ and negative part vK,σ−. Moreover, the property of upwind is also considered.Neglecting the high-order terms, we have the following approximate expression of the upwind formula [27]:(26)GK,σ≈uMvK,σ+−vK,σ−.In order to approximate the continuous fluxGK,σ on the cell-edge σ with second-order accuracy, we propose the following method.A local stencil is given in Figure2. For the cell K, M is the midpoint of an arbitrary edge and M1,M2 are the other two midpoints adjacent to it. We denote S△MM1K be the area of triangle MM1K and define(27)u¯K=∑M∈σ,σ∈ℰKωMuM,where ωM=S△MM1K+S△KM2M−1/∑M∈σ,σ∈εKS△MM1K+S△KM2M−1. For a special case (Figure 3), we set the cell K as a triangle and define u¯K=ωMuM+ωM1uM1+ωM2uM2, where(28)ωM=s1s,ωM1=s2s,ωM2=s3s,s=s1+s2+s3,s1=S△MM1K+S△KM2M−1,s2=S△MM1K+S△KM1M2−1,s3=S△KM1M2+S△MKM2−1.Figure 2
Local stencil 2.Figure 3
Local stencil of triangle.It is obvious thatu¯K is a second-order approximation to uK, i.e., uK−∑M∈σ,σ∈ℰKωMuM=Oh2 if the solution u∈C2K.Then, the approximation ofGK,σ on cell-edge σ can be defined as follows.Forσ∈ℰint, we define(29)GK,σ=BK,σuK−BL,σuL,where(30)BK,σ=vK,σ+uMu¯K≥0,BL,σ=vL,σ+uMu¯L≥0.Forσ∈ℰext, we define(31)GK,σ=BK,σuK−bK,σ,where(32)BK,σ=vK,σ+uMu¯K≥0,bK,σ=vK,σ−gM.


## 3.3. The Finite Volume Scheme and Picard Iteration

By using the definition of discretization of diffusive and convective flux, the finite volume scheme can be constructed as follows:(33)∑σ∈εKFK,σ+GK,σ=KFK+∑σ∈εK∩εextAK,σ,2uM+bK,σ,K∈J.(34)uMi=gMi,∀Mi∈∂Ω,where fK=fK and gMi=gxMi.LetU be the discrete unknown vector and AU be the coefficient matrix. A nonlinear algebraic system of the schemes (33) and (34) can be formed: AUU=F. The AU is assembled by the coefficients of diffusive term FK,σ and convective term GK,σ. We use the Picard nonlinear iteration method to solve the system: choose a small value ℰnon>0 and initial vector U0>0 and repeat for nonlinear iteration index s=1,2,…,:(1)
SolveAUs−1Us=F(2)
Stop ifAUsUs−F2≤ΕnonAU0U0−F2The linear algebraic system with coefficient matrixAUs−1 is solved by the biconjugate gradient stabilized (BiCGSTab) method, and the linear iterations are terminated when relative norm of the initial residual becomes smaller than εlin.


## 3.4. The Algorithm

In this subsection, we describe the detailed algorithm.Step 1. InitializeU0>0, εnon, and εlin.Step 2. Whens=0,(a)
computeAK,σ0 and BK,σ0;(b)
compute initial residualAU0U0−F2.Step 3. Whens=1,2,…,(a)
solveAUs−1Us=F;(b)
correctUs, see Remark 1;(c)
computeuPs and correct uMs, ∀M∈σ,σ∈Ε,P∈N, see Remark 2;(d)
computeAK,σs and BK,σs;(e)
compute residualAUsUs−F2;(f)
whetherAUsUs−F2≤εnonAU0U0−F2, if true, then go to (Step 4), otherwise, go to (Step 3).Step 4. Stop.Remark 1.
For∀uKs, if uKs>maxuMs−1,M∈σ,σ∈ΕK, let uKs=u¯Ks−1.Remark 2.
The value of∀uMs can be obtained by equation (22). If ∀uMs>maxuKs,uLs,uAs,uBs, where the common edge σ of two cells K and L is also denoted by AB (Figure 1). Let uMs=uMs−1.
It should be noted that the algorithm in Remark1 is important to avoid the numerical solution beyond the upper bound where the numerical results need to depend on nonnegative initial values of nonlinear iteration.Theorem 1.
LetF≥0,U0≥0, and linear systems in Picard iterations are solved exactly. Then,(35)Us≥0,s=1,2,3,…,.The detailed proof of positivity is given in the study by Yuan and Sheng [11].Now, we state our conclusion, which says that our scheme can avoid the numerical solution beyond the upper bound. Denoteumax=max0,uK,uM,∀K∈T,∀M∈ε.We assumeuK0=umax. Using Remark 1, we know that uK0≤maxuM,M∈σ,σ∈εK0.


## 4. Numerical Experiments

In order to demonstrate the accuracy and robustness of the scheme, we test several problems and takeεnon=1.0e−6 and εlin=1.0e−10. The convergence order can be obtained by the following formula:(36)Order=logErrorN1/ErrorN2logN2/N1,where N1 and N2 represent different number of cells, and ErrorN1 and ErrorN2 are the corresponding L2 errors.


### 4.1. The Problem with Anisotropic Diffusion Tensor

Consider the problems (1) and (2) with Dirichlet boundary condition on Ω=0,1×0,1, and take v=−1,1T. The exact solution is(37)ux,y=cosπx2ey,and the diffusion coefficient is(38)κ=100ε00ε.First, we test the accuracy of our scheme on random quadrilateral meshes shown in Figure4. Table 1 gives the corresponding L2 error and the numbers of nonlinear iteration numbers itnon# with a different parameter ε. We can see that our scheme obtains second-order accuracy for the solution and at least first-order accuracy for the flux. The average number of nonlinear iterations is 36.8 when ε=10−6. However, for ε=1.0, the corresponding number increases while the number of cell increases.Figure 4
The random quadrilateral meshes.Table 1
Accuracy for the problem with anisotropic coefficient.
Cells1445762304921636864ε=1.0ε2u3.11E−37.50E−42.09E−45.46E−51.42E−5Order—2.051.841.941.94ε2F1.06E+03.94E−11.76E−18.28E−23.96E−2Order—2.071.542.030.99itnon#5585116126153ε=10−6ε2u5.28E−31.24E−33.18E−48.06E−52.13E−5Order—2.091.961.981.92ε2F7.93E−31.99E−35.31E−41.34E−43.57E−5Order—2.071.542.030.99itnon#3640363339


### 4.2. The Problem with Discontinuous Coefficient

Consider the problems (1) and (2) with Dirichlet boundary condition on Ω=0,1×0,1, and take v=1,1T. The exact solution is(39)ux,y=sinπ2x+sinπ2y,x<12,2cπx−1/24+sinπ2y+22,x≥12,and the diffusion coefficient is(40)κ=c0×ε,x<12,ε,x≥12,where c0=40.We test this problem on random triangle meshes shown in Figure5. The numerical results with a different parameter ε=1.0,10−6 are given in Table 2. We can see that our scheme almost obtain second-order accuracy for the solution and at least first-order accuracy for the flux. The average numbers of nonlinear iterations are 25 and 23 for ε=1.0 and ε=1.0,10−6, respectively.Figure 5
The random triangle meshes.Table 2
Accuracy for the problem with discontinuous coefficient.
εCells2881152460818432737281.0ε2u4.87E−41.39E−43.82E−58.68E−62.13E−6Order—1.811.862.142.03ε2F7.83E−23.04E−21.39E−25.51E−32.50E−3Order—2.071.542.030.99itnon#232426262610−6ε2u4.54E−41.18E−43.13E−57.93E−62.01E−6Order—1.941.911.981.98ε2F6.60E−41.63E−44.15E−51.04E−52.63E−6Order—2.071.542.030.99itnon#2324232323Then, in order to illustrate the efficiency of our scheme, the comparison of accuracy between the studies by Lan et al. and Zhang et al. [19, 24] and present scheme is given in Table 3. As can be seen from the table, the accuracy is 1.5 in the study by Zhang et al. [24] and more than 2.0 in the study by Lan et al. [19] and our new scheme. However, the computed results in our scheme can approximate the exact solution more well.Table 3
Comparison of accuracy with discontinuous coefficientε=10−5.
MethodMeshCells2561024409616384Zhang et al. [24]Uniformε2u1.10E − 13.98E-21.42E − 25.02E − 3Order—1.471.491.50Randomε2u1.10E − 13.95E-21.38E − 25.11E − 3Order—1.481.521.43Lan et al. [19]Uniformε2u4.23E − 39.77E-42.31E − 45.50E − 5Order—2.122.082.07Randomε2u6.87E − 31.57E-33.78E − 48.96E − 5Order—2.132.052.08PresentUniformε2u9.34E − 41.82E − 43.63E − 58.43E − 6Order—2.362.322.11Randomε2u9.16E − 41.89E − 44.04E − 59.85E − 6Order—2.282.232.04


### 4.3. Positivity of Numerical Solutions

Now, we consider the problem (1)-(2) in the unit square Ω=0,12 with homogeneous Dirichlet boundary conditions. Take v⟶=1,1T, and set(41)κ=cosθsinθ−sinθcosθk100k2cosθ−sinθsinθcosθ,f=1,ifx,y∈38,582,0,otherwise.We takek1=1,k2=100, and θ=5π/6.The analytical solutionux,y is unknown, but the minimum principle states that it is nonnegative. It is a challenging task to solve it accurately because they can result in significant violation of the positivity and even produce a numerical solution with nonphysical oscillations. We will show that our new nonlinear scheme can also obtain the nonnegative solution.First, we test the new scheme on the random quadrilateral meshes with128×128 cells. The corresponding distribution is similarly shown in Figure 6), and the numerical solution is given in Figure 7. The minimum value is umin=0 and the maximum value is umax=6.7603×10−4, which show that our scheme preserves the positivity of the solution and does not produce any nonphysical oscillations. Then, we test it on the random triangular meshes. The computed results show that umin=0 and umax=6.7582×10−4.Figure 6
The random quadrilateral meshes.Figure 7
The numerical solutions on the random quadrilateral meshes.
(a)(b)These computed results illustrate that our new scheme preserves the positivity of numerical solutions and satisfies the discrete minimum principle.


### 4.4. Nonphysical Oscillations

We also consider the last nonsmooth anisotropic solution and compute it on the random quadrilateral meshes. Here, we resetf=0. The computational domain is a unit square with a hole, Ω=0,12∖4/9,5/92, so that the boundary ∂Ω is composed of two disjoint parts Γ0 and Γ1 as shown in Figure 8 where the number of cell is 72×72. Γ0 is the exterior boundary, and Γ1 is the interior boundary. We set g=0 on Γ0 and g=2 on Γ1.Figure 8
The random quadrilateral meshes with a hole.The numerical solutions on the random quadrilateral meshes are shown in Figure9. The computed results show that umin=6.94E−13 and umax=1.98. It means that the minimum 0 is attained on the Γ0, and the maximum 2 is attained on the Γ1. So, these computed results illustrate that our scheme can avoid the numerical solution beyond the upper bound and does not produce any nonphysical oscillations.Figure 9
The numerical results on the random quadrilateral meshes.
(a)(b)Then, the computed results without the correct method in Remarks1 and 2 are shown in Figure 10, and umin=3.87E−6 and umax=1.91. However, the numerical oscillations are produced in the computational domain.Figure 10
The numerical results without corrections on the random quadrilateral meshes.
(a)(b)


## 4.1. The Problem with Anisotropic Diffusion Tensor

Consider the problems (1) and (2) with Dirichlet boundary condition on Ω=0,1×0,1, and take v=−1,1T. The exact solution is(37)ux,y=cosπx2ey,and the diffusion coefficient is(38)κ=100ε00ε.First, we test the accuracy of our scheme on random quadrilateral meshes shown in Figure4. Table 1 gives the corresponding L2 error and the numbers of nonlinear iteration numbers itnon# with a different parameter ε. We can see that our scheme obtains second-order accuracy for the solution and at least first-order accuracy for the flux. The average number of nonlinear iterations is 36.8 when ε=10−6. However, for ε=1.0, the corresponding number increases while the number of cell increases.Figure 4
The random quadrilateral meshes.Table 1
Accuracy for the problem with anisotropic coefficient.
Cells1445762304921636864ε=1.0ε2u3.11E−37.50E−42.09E−45.46E−51.42E−5Order—2.051.841.941.94ε2F1.06E+03.94E−11.76E−18.28E−23.96E−2Order—2.071.542.030.99itnon#5585116126153ε=10−6ε2u5.28E−31.24E−33.18E−48.06E−52.13E−5Order—2.091.961.981.92ε2F7.93E−31.99E−35.31E−41.34E−43.57E−5Order—2.071.542.030.99itnon#3640363339


## 4.2. The Problem with Discontinuous Coefficient

Consider the problems (1) and (2) with Dirichlet boundary condition on Ω=0,1×0,1, and take v=1,1T. The exact solution is(39)ux,y=sinπ2x+sinπ2y,x<12,2cπx−1/24+sinπ2y+22,x≥12,and the diffusion coefficient is(40)κ=c0×ε,x<12,ε,x≥12,where c0=40.We test this problem on random triangle meshes shown in Figure5. The numerical results with a different parameter ε=1.0,10−6 are given in Table 2. We can see that our scheme almost obtain second-order accuracy for the solution and at least first-order accuracy for the flux. The average numbers of nonlinear iterations are 25 and 23 for ε=1.0 and ε=1.0,10−6, respectively.Figure 5
The random triangle meshes.Table 2
Accuracy for the problem with discontinuous coefficient.
εCells2881152460818432737281.0ε2u4.87E−41.39E−43.82E−58.68E−62.13E−6Order—1.811.862.142.03ε2F7.83E−23.04E−21.39E−25.51E−32.50E−3Order—2.071.542.030.99itnon#232426262610−6ε2u4.54E−41.18E−43.13E−57.93E−62.01E−6Order—1.941.911.981.98ε2F6.60E−41.63E−44.15E−51.04E−52.63E−6Order—2.071.542.030.99itnon#2324232323Then, in order to illustrate the efficiency of our scheme, the comparison of accuracy between the studies by Lan et al. and Zhang et al. [19, 24] and present scheme is given in Table 3. As can be seen from the table, the accuracy is 1.5 in the study by Zhang et al. [24] and more than 2.0 in the study by Lan et al. [19] and our new scheme. However, the computed results in our scheme can approximate the exact solution more well.Table 3
Comparison of accuracy with discontinuous coefficientε=10−5.
MethodMeshCells2561024409616384Zhang et al. [24]Uniformε2u1.10E − 13.98E-21.42E − 25.02E − 3Order—1.471.491.50Randomε2u1.10E − 13.95E-21.38E − 25.11E − 3Order—1.481.521.43Lan et al. [19]Uniformε2u4.23E − 39.77E-42.31E − 45.50E − 5Order—2.122.082.07Randomε2u6.87E − 31.57E-33.78E − 48.96E − 5Order—2.132.052.08PresentUniformε2u9.34E − 41.82E − 43.63E − 58.43E − 6Order—2.362.322.11Randomε2u9.16E − 41.89E − 44.04E − 59.85E − 6Order—2.282.232.04


## 4.3. Positivity of Numerical Solutions

Now, we consider the problem (1)-(2) in the unit square Ω=0,12 with homogeneous Dirichlet boundary conditions. Take v⟶=1,1T, and set(41)κ=cosθsinθ−sinθcosθk100k2cosθ−sinθsinθcosθ,f=1,ifx,y∈38,582,0,otherwise.We takek1=1,k2=100, and θ=5π/6.The analytical solutionux,y is unknown, but the minimum principle states that it is nonnegative. It is a challenging task to solve it accurately because they can result in significant violation of the positivity and even produce a numerical solution with nonphysical oscillations. We will show that our new nonlinear scheme can also obtain the nonnegative solution.First, we test the new scheme on the random quadrilateral meshes with128×128 cells. The corresponding distribution is similarly shown in Figure 6), and the numerical solution is given in Figure 7. The minimum value is umin=0 and the maximum value is umax=6.7603×10−4, which show that our scheme preserves the positivity of the solution and does not produce any nonphysical oscillations. Then, we test it on the random triangular meshes. The computed results show that umin=0 and umax=6.7582×10−4.Figure 6
The random quadrilateral meshes.Figure 7
The numerical solutions on the random quadrilateral meshes.
(a)(b)These computed results illustrate that our new scheme preserves the positivity of numerical solutions and satisfies the discrete minimum principle.


## 4.4. Nonphysical Oscillations

We also consider the last nonsmooth anisotropic solution and compute it on the random quadrilateral meshes. Here, we resetf=0. The computational domain is a unit square with a hole, Ω=0,12∖4/9,5/92, so that the boundary ∂Ω is composed of two disjoint parts Γ0 and Γ1 as shown in Figure 8 where the number of cell is 72×72. Γ0 is the exterior boundary, and Γ1 is the interior boundary. We set g=0 on Γ0 and g=2 on Γ1.Figure 8
The random quadrilateral meshes with a hole.The numerical solutions on the random quadrilateral meshes are shown in Figure9. The computed results show that umin=6.94E−13 and umax=1.98. It means that the minimum 0 is attained on the Γ0, and the maximum 2 is attained on the Γ1. So, these computed results illustrate that our scheme can avoid the numerical solution beyond the upper bound and does not produce any nonphysical oscillations.Figure 9
The numerical results on the random quadrilateral meshes.
(a)(b)Then, the computed results without the correct method in Remarks1 and 2 are shown in Figure 10, and umin=3.87E−6 and umax=1.91. However, the numerical oscillations are produced in the computational domain.Figure 10
The numerical results without corrections on the random quadrilateral meshes.
(a)(b)


## 5. Conclusion

The aim of this paper is to build a nonlinear finite volume scheme preserving positivity for solving the 2D convection-diffusion equation on arbitrary convex polygonal meshes. We first develop the nonlinear positive finite volume scheme. Then, a corrected method is proposed, and the numerical solution beyond the upper bound can be avoided.Our scheme includes only cell-centered unknowns. Numerical results show that our new scheme obtains second-order accuracy for the solution and first-order accuracy for the flux. In addition, it can not only keep the positivity but also do not produce any oscillation.



---

*Source: 7343716-2020-08-21.xml*