# Positive Definite Advection Transport Algorithm for Conservation Law Equations on Nonuniform Irregular Grids

Authors: Xinpeng Yuan; Chunguang Xiong
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The multidimensional positive definite advection transport algorithm (MPDATA) is an important numerical method for the computation of atmospheric dynamics. MPDATA is second-order accurate, positive definite, conservative, and computationally efficient. However, the method is problematic in which it results in a loss of precision when computing a nonuniform irregular grid. Furthermore, research revealed two reasons for this problem. On the one hand, numerical discretization of boundary derivatives of the finite-volume method is incompatible with nonuniform meshes (or grids); on the other hand, the up-wind scheme of staggered grids is not applicable to the calculation of irregular grids. We overcome these two problems by using the multipoint Taylor expansion method to obtain a boundary derivative numerical approximation scheme that does not depend on the grid structure. Furthermore, combined with the well-balance central-upwind scheme, a positive definite advection scheme for irregular meshes is proposed. Then, the positivity of the new numerical scheme is analyzed. Finally, the result of this study is verified by numerical simulation.


---

## Body

## 1. Introduction

Numerical modeling of atmospheric phenomena often necessitates solving the advection equation for positive definite scalar functions. Owing to the large spatial variation range, a large gradient and strong discontinuity exist in the distribution of these scalar fields [1]. These two factors lead to oscillation and negative values and cause unacceptable dispersion errors in the solution of ordinary advection transport schemes. These problems typically destroy the inherent spatial distribution and conservation of scalar substances such as water vapor. In addition, the positive definite scheme [2] is also of great significance in the calculation of advection transport of atmospheric aerosols and pollutants.Smolarkiewicz introduced a predictor-iterative-corrector scheme in upwind schemes and proposed the multidimensional positive definite advection transport algorithm (MPDATA) [3, 4] for approximating the advection terms in fluid equations. MPDATA is second-order accurate, positive definite, conservative, and computationally efficient. Since the first article in 1983 [3], the MPDATA has developed into a calculation method of multiple equations and global atmospheric model schemes [5]. However, its theoretical core is still a continuation of the ideas of 1983. Through numerical experiments, we found that the existing MPDATA for the nonuniform grid domains will be a loss of accuracy. In fact, for the nonuniform grid, many existing numerical schemes cannot achieve the design accuracy. In [6], Xie studies such precision problems. Numerical solutions on spheres are the key to developing global atmospheric models. In particular, uniform mesh distribution for spherical atmospheric dynamics simulation does not exist.Therefore, we hope to improve the numerical precision of the MPDATA on the nonuniform grid. The new MPDATA has two differences from the former one. On the one hand, it is the calculation of boundary derivative. Here, we use the multipoint Taylor expansion method to calculate the boundary derivative, which will be introduced in Section2. On the other hand, the upwind format of the former MPDATA scheme is replaced with a suitable form for the nonuniform grid. In recent years, Kurganov et al. integrated over nonuniform polygonal control volumes of different shapes and derived a fully discrete central-upwind scheme to present a family of central-upwind schemes on general triangulations. Those schemes are simpler than their upwind counterparts, have a genuinely multidimensional structure, and can be used to solve problems on irregular domains [7, 8]. Based on the work of Kurganov et al., we will improve the MPDATA for the nonuniform grid.In Section2, we will analyze the accuracy of the former MPDATA numerical scheme and propose our numerical precision solution. Furthermore, numerical simulation is employed to verify the result that forms the conclusion of this study.


## 2. Numerical Methods

### 2.1. Analysis of the Previous MPDATA

The irregular meshes discussed in this paper all refer to nonuniform irregular meshes. In this section, we analyze the reasons why previous MPDATA methods are not applicable to irregular grids. The one-dimensional advection equation is(1)∂∂tψ+∂∂xuψ=0,where ψ is a nondiffusive scalar quantity and u is the one-dimensional velocity. The former MPDATA scheme is proposed by Smolarkiewicz and Margolin [3], which is(2)ψi∗=ψin−Fψin,ψi+1n,ui+1/2n−Fψi−1n,ψin,ui−1/2n,(3)ψin+1=ψi∗−Fψi∗,ψi+1∗,u˜i+1/2−Fψi−1∗,ψi∗,u˜i−1/2,where(4)Fψin,ψi+1n,ui+1/2n=ui+1/2n+ui+1/2nψin+ui+1/2n−ui+1/2nψi+1nΔt2Δx.Here,ψin is the value of ψ at the ith frid point for the nth time step. u˜i+1/2 is the numerical approximation of ud at xi+1/2, and it is(5)u˜i+1/2=ui+1/2Δx−Δtui+1/22ψi+1∗+ψi∗+εψi+1∗−ψi∗Δx,where ε is a small value (e.g., 10−15) to ensure u˜=0 when ψi+1∗=ψi∗=0.For a regular grid, the numerical schemes (2) and (3) are an effective way to obtain positive definite solutions (see Section 3.1). However, neither (2) nor (3) are able to provide the second-order numerical accuracy solutions for irregular grids. There are two reasons for this problem. The first is the compatibility of numerical discretion of the “diffusion velocity” ud. The second is the compatibility of the first-order numerical scheme. To overcome the first problem, we have the following propositions.Proposition 1.
In (5), the numerical discretation of ∂ψ/∂xi+1/2 is(6)∂ψ∂xi+1/2=ψi+1∗−ψi∗Δx+OΔx,for regular grids, and discretation (6) is the first order. For an irregular grid, the discretation of (6) prevents step (3) from reaching first order. As a result, the MPDATA schemes (2) and (3) cannot reach the second order.Proof.
For the irregular grid, discretation (6) is(7)∂ψ∂xi+1/2=ψi+1∗−ψi∗Δx+C0,where C0 is a constant. This enables us to obtain the relationship between the numerical approximation values and the exact values:(8)ψi∗=ψxi+OΔx,(9)u˜=uxi+C0.
Substituting the abovementioned equations (8) and (9) into (4), we have(10)Fψi∗,ψi+1∗,u˜i+1/2=u˜xi+1/2+u˜xi+1/2ψxi+u˜xi+1/2−u˜xi+1/2ψxi+1+u˜xi+1/2OΔx+ψxiC1+ψxi+1C2+OΔx2Δt2Δx,where C1 and C2 are constant. Using the same method for ∂ψ/∂xi−1/2, we have(11)Fψi∗,ψi+1∗,u˜i+1/2−Fψi−1∗,ψi∗,u˜i−1/2=Fψxi,ψxi+1,u˜xi+1/2−Fψxi−1,ψxi,u˜xi−1/2+u˜xi+1/2−u˜xi−1/2Δt+C1′ψxi−1+C2′ψxi+C3′ψxi+1+OΔx2,where C1′, C2′, and C3′ are constant. When Δx⟶0,(12)C1′ψxi−1+C2′ψxi+C3′ψxi+1≠OΔx.
Thus, step (3) cannot become first order, and it cannot calculate the affect of diffusion effectively. Thus, the numerical schemes (2) and (3) cannot become the second order.
To obtain accurate numerical schemes on an irregular grid, a new approximation method of∂ψ/∂xi±1/2 is given here. The mesh structure around xi+1/2 is described in Figure 1.Figure 1
Schematic of the mesh structure aroundxi+1/2.Then, we expandψi−1, ψi, ψi+1,andψi+2 at the point ψi+1/2, and we have(13)ψs=ψi+1/2+xs−xi+1/2∂ψ∂xi+1/2+OΔx2.Here,s=i−1,i,i+1,i+2. A set of coefficients κ1i+1/2, i=1,2,3,4, need to be obtained, which satisfies(14)1111xi−1−xi+1/2xi−xi+1/2xi+1−xi+1/2xi+2−xi+1/2κ1i+1/2κ2i+1/2κ3i+1/2κ4i+1/2=01.It is easy to see that the rank of the coefficient matrix and the extended matrix for linear equations (11) are the same, and the linear equations (11) have infinite solutions. Here, the least-norm solutions are used (see Appendix). Then, we have the approximation for ∂ψ/∂x:(15)∂ψ∂xi+1/2=ψi−1⋅κ1i+1/2+ψiκ2i+1/2+ψi+1κ3i+1/2+ψi+2κ4i+1/2+OΔx.with the numerical approximation (), we have(16)ψi∗=ψxi+OΔx,u˜=uxi+OΔx.For the same reason as (11), we have(17)Fψi∗,ψi+1∗,u˜i+1/2−Fψi−1∗,ψi∗,u˜i−1/2=Fψxi,ψxi+1,u˜xi+1/2−Fψxi−1,ψxi,u˜xi−1/2+u˜xi+1/2−u˜xi−1/2Δt+C˜1ψxi−1Δt+C˜2ψxiΔt+C˜3ψxi+1Δt+OΔx2.Becauseu˜ is a continuous function, we have(18)C˜1ψxi−1Δt+C˜2ψxiΔt+C˜3ψxi+1Δt∼OΔxΔx⟶0.Thus, using (12) instead of (5) will enable step (3) to reach the first order, whereupon it would be able to calculate the affect of diffusion effectively. Furthermore, (2) and (3) can achieve second-order accuracy.Before we begin the discussion for a two-dimensional situation, we have another proposition to state.Proposition 2.
The compatibility of the first-order numerical scheme is important for the MPDATA method.
Generally speaking, the MPDATA schemes (2) and (3) belong to the class of staggered central-upwind schemes, which are excellent tools for solving various complex problems on regular domains. However, in practice complicated geometries exist, and the use of irregular triangular meshes could be advantageous or even unavoidable [7]. Staggered central-upwind schemes are sensitive to mesh effects. For example, Küther proved that Lax–Friedrichs-type finite-volume schemes estimate the error to be of the order of 0.25 [9, 10].
For the first-order scheme,(19)ψin+1=ψin−Fψin,ψi+1n,ui+1/2n−Fψi−1n,ψin,ui−1/2n.
With the Taylor expansionψi−1, ψi, ψi+1,andψi+2 at the point ψi+1/2, we have(20)Fψi,ψi+1,ui+1/2−Fψi−1,ψi,ui−1/2=ui+1/2−ui−1/2ψiΔtΔx+Δt4Δx∂ψ∂xiΔxi+Δxi+1ui+1/2−ui+1/2+Δxi−1+Δxiui−1/2+ui−1/2+OΔx2.
Whenu is a continuous constant, we have(21)Fψi,ψi+1,ui+1/2−Fψi−1,ψi,ui−1/2=u∂ψ∂xiΔxi−1+2Δxi+Δxi+14ΔxiΔt.
Based on the abovementioned equation, it is clear that if the mesh is distorted, a few points are incompatible. Whenu is not a discontinuous function, equation (21) implies that scheme (21) is largely incompatible.
We computed the irregular grid by carrying out a comparative simulation in Section 6.1. The results in Tables1–3 indicate that the standard MPDTA method reduces the order when calculated on an irregular grid.Table 1
Former MPDATA for regular grids.
N
L
1 error
L
1 order
L
2 error
L
2 order
L
n error
L
∞ order
16
4.333
E
−
02
—
3.330
E
−
02
—
3.313
E
−
02
—
32
1.115
E
−
02
1.96
8.513
E
−
03
1.97
9.022
E
−
03
1.87
64
2.717
E
−
03
2.04
2.132
E
−
03
2.00
2.545
E
−
03
1.82
128
6.648
E
−
04
2.03
5.280
E
−
04
2.01
6.412
E
−
04
1.98
256
1.640
E
−
04
2.02
1.308
E
−
04
2.01
1.596
E
−
04
2.00Table 2
Former MPDATA for irregular grids.
N
L
1 error
L
1 order
L
2 error
L
2 order
L
n error
L
∞ order
16
2.181
E
−
01
—
1.834
E
−
01
—
2.471
E
−
01
—
32
6.887
E
−
02
1.66
6.200
E
−
02
1.57
1.070
E
−
01
1.21
64
2.188
E
−
02
1.65
2.176
E
−
02
1.51
4.905
E
−
02
1.13
128
7.550
E
−
03
1.53
8.901
E
−
03
1.29
2.385
E
−
02
1.04
256
3.099
E
−
03
1.28
4.125
E
−
03
1.11
1.184
E
−
02
1.01Table 3
New MPDATA for irregular grids.
N
L
1 error
L
1 order
L
2 error
L
2 order
L
n error
L
∞ order
16
1.976
E
−
01
—
1.568
E
−
01
—
1.830
E
−
01
—
32
5.744
E
−
02
1.78
4.765
E
−
02
1.72
5.213
E
−
02
1.81
64
1.450
E
−
02
1.99
1.272
E
−
02
1.91
1.536
E
−
02
1.76
128
3.680
E
−
03
1.98
3.240
E
−
03
1.97
3.935
E
−
03
1.97
256
9.221
E
−
04
2.00
8.146
E
−
04
1.99
9.830
E
−
04
2.00


### 2.2. New MPDATA for Irregular Grids in One-Dimensional Problems

Base on the abovementioned analysis, we propose a new MPDATA method for irregular grids in the one-dimensional situation. One of the grid elements being studied is denoted asTj. After time step Δt, the discontinuity around xj±1/2 can be described as ajLinΔt and ajRinΔt for inward Tj, where ajLinΔt is around xj−1/2 and ajRinΔt is around xj+1/2. For outward, ajLout is around xj−1/2 and ajRout is around xj+1/2. Using this information, the control volume can be constructed as follows. As shown in Figure 2, the inside domain Dj; the boundary domain Dj±1/2, which can be divided into the inward domain Dj±1/2in (Dj−1/2in=ajLinΔt and Dj+1/2in=ajRinΔt) and the outward domain Dj±1/2out (Dj−1/2out=ajLoutΔt and Dj+1/2out=ajRoutΔt). All domains around element Tj are described as illustrated in Figure 2.Figure 2
Domains on control volumeTj.LetujL be the velocity in domain Dj−1/2in and ujR be the velocity in domain Dj+1/2in. For the advection equations (1), the discontinuity propagate with inward and outward Tj local speeds, which can be estimated by(22)ajLin=−min−uj−1R,−ujL,0,ajRin=−minuj+1L,ujR,0,(23)ajLout=max−uj−1R,−ujL,0,ajRout=maxuj+1L,ujR,0.Then, we can compute the averages over the corresponding domains. The domain average forDj×tn,tn+1 is(24)ϕ¯jn+1=1Dj∫Djψjnxdx−∫tntn+1uψxj−1/2+ajRinΔtxj+1/2−ajRinΔtdt=ψ¯jn−1Dj∫tntn+1uψxj−1/2+ajRinΔtxj+1/2−ajRinΔtdt+OΔx,and the domain average for Dj+1/2×tn,tn+1. is(25)ψ¯j+1/2n+1=1Dj+1/2∫Dj+1/2inψjnxdx+∫Dj+1/2outψj+1nxdx−∫tntn+1uj+1Lψj+1−ujRψjdt=1Dj+1/2Dj+1/2inψ¯jn+Dj+1/2outψ¯jn+1−∫tntn+1uj+1Lψj+1−ujRψjdt+OΔx.On the other hand, we can have the relationship(26)ψ¯jn=1Tj∫Tjψx,tndx=1Tj∫Tjϕnxdx.The time derivative ofψjt is expressed with the help of (24) as(27)ddtψ¯jt=limΔt⟶0ψ¯jn+1−ψ¯jnΔt=limΔt⟶01Δt1Tj∫Tjϕn+1xdx−ψ¯jn.Because the width of the Riemann fans is bounded byajLout+ajLinΔt or ajRout+ajRinΔt, we obtain(28)ϕj±1/2n+1x=ϕ¯j±1/2n+1+OΔt,where x∈Dj. From (28) and the conservation of ϕjn+1, we derive(29)ddtψ¯jt=limΔt⟶01Δt1Tj∫Dj−1/2inϕn+1xdx+limΔt⟶01Δt1Tj∫Djϕn+1xdx+limΔt⟶01Δt1Tj∫Dj+1/2inϕn+1xdx−limΔt⟶01Δtψ¯jn=ajLinTjlimΔt⟶0ϕ¯j−1/2n+1+limΔt⟶01ΔtDjTjϕ¯jn+1−ψ¯jn+ajRinTjlimΔt⟶0ϕ¯j+1/2n+1.With the help of (24) and (25), we obtain(30)limΔt⟶01ΔtDjTjϕ¯jn+1−ψ¯jn=−ajRinψjn+ajLinψjnTj−ujR−ujLψjnTj+OΔx,and(31)limΔt⟶0ϕ¯j+1/2n+1=ajRinψ¯j+ajRoutψ¯j+1ajRin+ajRout−−ujRψ¯j+uj+1Lψ¯j+1ajRin+ajRout+OΔx.With the help of the forward Euler method, the substitution of (30) and (31) in (29) results in the new first-order scheme:(32)ψ¯jn+1=ψ¯jn−ΔtΔxjHj+1/2ψ¯j,ψ¯j+1,u−Hj−1/2ψ¯j−1,ψ¯j,u,where(33)Hj+1/2=ajRoutujRψ¯jn+ajRinuj+1Lψ¯j+1najRout+ajRin−ajRoutajRinajRout+ajRinψ¯j+1n−ψ¯jn,Hj−1/2=ajLoutujLψ¯jn+ajLinuj−1Rψ¯j−1najLout+ajLin+ajLoutajLinajLout+ajLinψ¯j−1n−ψ¯jn.Using the Taylor expansionxj+1, xj at xj+1/2(34)ψj+1n=ψj+1/2n+12Δxj+1∂∂xψj+1/2n+OΔx2,(35)ψjn=ψj+1/2n−12Δxj∂∂xψj+1/2n+OΔx2,and xj, xj−1 at xj−1/2:(36)ψj−1n=ψj−1/2n−12Δxj−1∂∂xψj−1/2n+OΔx2,(37)ψjn=ψj−1/2n+12Δxj∂∂xψj−1/2n+OΔx2.Taking those Taylor expressions (34)–(37) into Hj+1/2 and Hj−1/2, we can obtain(38)Hj+1/2=ajRoutujR+ajRinuj+1LajRout+ajRinψj+1/2n+−ajRoutujRΔxj+ajRinuj+1LΔxj+1ajRout+ajRin12−Δxj+1+Δxj2ajRoutajRinajRout+ajRin∂∂xψj+1/2n+OΔx2,Hj−1/2=ajLoutujL+ajLinuj−1RajLout+ajLinψj−1/2n+ajLoutujLΔxj−ajLinuj−1RΔxj−1ajLout+ajLin12−Δxj−1+Δxj2ajLoutajLinajLout+ajLin∂∂xψj−1/2n+OΔx2.With the Taylor expansions on time,(39)ψjn+1=ψjn+1/2+12Δt∂ψ∂tjn+1/2+OΔt2,ψjn=ψjn+1/2−12Δt∂ψ∂tjn+1/2+OΔt2=ψjn+1/2+12Δt∂uψ∂xjn+1/2+OΔt2.Scheme (32) becomes(40)∂ψ∂tjn+1/2=1ΔxjajRoutujR+ajRinuj+1LajRout+ajRinψj+1/2n+1/2−ajLoutujL+ajLinuj−1RajLout+ajLinψj−1/2n+1/2+1ΔxjE˜,where(41)E˜:=u˜j+1/2ψj+1/2n+1/2−u˜j−1/2ψj−1/2n+1/2,u˜j+1/2:=1ψj+1/2n+1/2ajRoutujR+ajRinuj+1LajRout+ajRin12Δt∂∂xuψj+1/2n+1/2+1ψj+1/2n+1/2−ajRoutujRΔxj+ajRinuj+1LΔxj+1ajRout+ajRin12−Δxj+1+Δxj2ajRoutajRinajRout+ajRin∂∂xψj+1/2n,u˜j−1/2:=1ψj−1/2n+1/2ajLoutujL+ajLinuj−1RajLout+ajLin12Δt∂∂xuψj−1/2n+1/2+1ψj−1/2n+1/2ajLoutujLΔxj−ajLinuj−1RΔxj−1ajLout+ajLin12+Δxj−1+Δxj2ajLoutajLinajLout+ajLin∂∂xψj−1/2n.Whenu is continuous around xj−1/2 and xj+1/2, noticing one of aout or ain is zero. Then,(42)ajRoutujR+ajRinuj+1LajRout+ajRin⟹uj+1/2,ajLoutujL+ajLinuj−1RajLout+ajLin⟹uj−1/2.Thus, equation (40) can be viewed as(43)∂ψ∂t=∂uψ∂x+∂u˜ψ∂x.Here, the second term of (43) is the “pseudoadvection term.” (32) is the second-order scheme of equation (43). Therefore, we have the new MPDATA for irregular grids:(44)ψj∗=ψjn−Hj+1/2ψjn,ψj+1n,u−Hj−1/2ψj−1n,ψjn,uΔtΔxj,(45)ψjn+1=ψj∗−Hj+1/2ψj∗,ψj+1∗,u˜d−Hj−1/2ψj−1∗,ψj∗,u˜dΔtΔxj.Because the new MPDATA scheme depends on the two-step scheme, we can analyze the positivity of (32) instead of the schemes (44) and (45). Thus, we have the following claim for the new MPDATA schemes (44) and (45).Theorem 1.
Consider system (1) and the central-upwind discrete scheme (29), where ajLin, ajRin, ajLout, and ajRout are defined in (22) and (23). Then, for all j, ψ¯jn+1≥0, provided that Δt≤Δx/a, where Δx:=minjΔxj and a:=maxjmaxuj+1L,ujR.Proof.
For convenience of expression, let(46)λj=ΔtΔxj.
From (32), we can have(47)ψ¯jn+1=ψ¯jn−λjHj+1/2ψ¯jn,ψ¯j+1n,u−Hj−1/2ψ¯j−1n,ψ¯jn,u=12−λjajRoutajRin+ujRajRin+ajRoutψ¯jn+12+λjajLout−ajLin+ujLajLin+ajLoutψ¯jn+λjajRin−uj+1L+ajRoutajRin+ajRoutψ¯j+1n+λjajLinuj−1R+ajLoutajLin+ajLoutψ¯j−1n.
In the above equation (47), under ψ¯jn≥0, it is easy to see that(48)12+λjajLout−ajLin+ujLajLin+ajLoutψ¯jn≥0,λjajRin−uj+1L+ajRoutajRin+ajRoutψ¯j+1n≥0,λjajLinuj−1R+ajLoutajLin+ajLoutψ¯j−1n≥0,12−λjajRoutajRin+ujRajRin+ajRoutψ¯jn≥0,whenλj≤12ajRout.
Thus, we can claim that, for allj, ψ¯jn+1≥0, provided that Δt≤Δx/a, where Δx:=minjΔxj and a:=maxjmaxuj+1L,ujR.


### 2.3. New MPDATA for Irregular Grids in a Two-Dimensional Problem

For the two-dimensional advection problem,(49)∂ψ∂t=−∇⋅vψ,where ψx,t is a nondiffusive scalar field assumed to be nonnegative at t=0, and the v=vx is the prescribed flow. The finite-volume MPDATA is considered on equation (48). A triangular partition ∪jTj is used on the computational domain. The area of Tj is Tj. All domains around the triangular cell Tj are illustrated in Figure 3. The neighborhood cell near Tj is denoted as Tjk, k=1,2,3. The three edges of cell Tj are denoted as Sjk. The discontinuous long edges Sjk can move inward, in which case Tj is denoted as ajkinΔt, and outward such that Tj is represented as ajkoutΔt. The inside domain is Dj. The three nonconvex domains surrounding the corners are Ejk, k=1,2,3.Figure 3
Triangular mesh of central-upwind scheme.Similar to the one-dimensional problem,ajkin and ajkout can be interpreted as the possible discontinuities of ψj along the sides of Tj propagating with different inward, ajkin, and outward, ajkout, directional local speeds k=1,2,3. They can be estimated as(50)ajkin=−minminM∈Tj∩Tjkλ1VjkTj,minM∈Tj∩Tjkλ1VjkTjk,0,ajkout=maxmaxM∈Tj∩TjkλNVjkTj,maxM∈Tj∩TjkλNVjkTj,0,where λ1Vjk≤⋯≤λNVjk are the N eigenvalues of the matrix Vjk defined by(51)VjkT=njk⋅vjkT,and njk is the outer unit normals to the corresponding sides of Tj of length hjk, k=1,2,3. Following the similar method in Section 2 and the theory developed by Kurganov and Petrova [7], we obtain a consistent, accurate, and effective well-balanced numerical scheme:(52)ψ¯jn+1−ψ¯jnΔt=−1Tj∑k=13hjkajkinvjk−ψ¯jkn+ajkoutvjk+ψ¯jnajkin+ajkout⋅njk+1Tj∑k=13hjkajkinajkoutajkin+ajkoutψ¯jkn−ψ¯jn,where vjk+ and vjk− are advection velocities in domain Djk+ and domain Djk−. Here, (52) is a first-order scheme. Then, expanding ψj and ψjk at the middle point of the edge intersecting face Sj,(53)ψ¯jk=ψSj+ϑjkψSj+OΔx2+Δy2,where ϑjk=xjk−xSj∂/∂x+yjk−ySj∂/∂y.(54)ψ¯j=ψSj+ϑjψSj+OΔx2+Δy2,where ϑj=xj−xSj∂/∂x+yj−ySj∂/∂y. We substitute (53) and (54) into (52), and we have(55)ψ¯jn+1−ψ¯jnΔt=−1Tj∑k=13hjkajkinvjk−ψSjn+ajkoutvjk+ψSjnajkin+ajkout⋅njk−1Tj∑k=13hjkajkinvjk−ϑjkψSjn+ajkoutvjk+ϑjψSjnajkin+ajkout⋅njk+1Tj∑k=13hjkajkinajkoutajkin+ajkoutϑjkψSjn−ϑjψSjn.Expandingψn+1 and ψn at ψn+1/2,(56)ψn+1=ψn+1/2+12∂ψ∂tn+1/2Δt+OΔt2,ψn=ψn+1/2−12∂ψ∂tn+1/2Δt+OΔt2=ψn+1/2+12Δtψ∇vn+1/2+12Δtv∇ψn+1/2+OΔt2.Then, we have the relationship(57)∂ψ¯∂tn+1/2=−1Tj∑k=13hjkajkinvjk−+ajkoutvjk+ajkin+ajkoutψSjn+1/2⋅njk+Error,where(58)Error=−1Tj∑k=13hjkajkinvjk−+ajkoutvjk+ajkin+ajkout12Δtψ∇vSjn+1/2⋅njk−1Tj∑k=13hjkajkinvjk−+ajkoutvjk+ajkin+ajkout12Δtv∇ψSjn+1/2⋅njk−1Tj∑k=13hjkajkinvjkϑjkψSjn+ajkoutvjϑjψSjnajkin+ajkout⋅njk+1Tj∑k=13hjkajkinajkoutajkin+ajkout⋅njkϑjkψSjn−ϑjψSjn⋅njk.Noticing that(59)ajkinvjk−+ajkoutvjk+ajkin+ajkout⟹vSj,(59) can be can be understood as the average value of v in the small rectangle Djk (see Figure 3). Furthermore, we can assume t=nΔt, and (57) can be written as(60)∂ψ∂t=−∇⋅vψ+∇⋅v˜dψ,where(61)v˜jkd:=−ajkinvjk−+ajkoutvjk+ajkin+ajkout12Δt∇⋅vSjn+1/2−1ψajkinvjk−+ajkoutvjk+ajkin+ajkout12Δtv⋅∇ψSjn+1/2−1ψajkinvjk−ϑjkψSjn+ajkoutvjk+ϑjψSjnajkin+ajkout+1ψajkinajkoutajkin+ajkout⋅njkϑjkψSjn−ϑjψSjn.Here, the gradient∇ψ at the edges can be obtained by the combination of the physical values around the edge. The mesh structure, which is used to obtain the gradient ∇ψj1, around the midpoint of an edge Sj1 is described in Figure 4.Figure 4
Schematic of the mesh structure aroundSj1.Here, we useψj1 to represent the physical value at point Mj1. Using the Taylor expression, expand ψ¯j, ψ¯j1, ψ¯j2, ψ¯j3, ψ¯j11, and ψ¯j12 at point Mj1 to second-order, and we have(62)ψ¯j=ψSj1+xj−xSj1∂∂xψSj1+yj−ySj1∂∂yψSj1+OΔx2+Δy2,ψ¯j1=ψSj1+xj1−xSj1∂∂xψSj1+yj1−ySj1∂∂yψSj1+OΔx2+Δy2,ψ¯j2=ψSj1+xj2−xSj1∂∂xψSj1+yj2−ySj1∂∂yψSj1+OΔx2+Δy2,ψ¯j3=ψSj1+xj3−xSj1∂∂xψSj1+yj3−ySj1∂∂yψSj1+OΔx2+Δy2,ψ¯j11=ψSj1+xj11−xSj1∂∂xψSj1+yj11−ySj1∂∂yψSj1+OΔx2+Δy2,ψ¯j12=ψSj1+xj12−xSj1∂∂xψSj1+yj12−ySj1∂∂yψSj1+OΔx2+Δy2.Next, use the relationship(63)∑i=05αiφSi=0,∑i=05αiyi−ySj1=0,∑i=05αixj1−xSj1=1,to obtain ∂/∂xφSj1, and relationship(64)∑i=05αiφSi=0,∑i=05αiyi−ySj1=0,∑i=05αixj1−xSj1=1,for ∂/∂yφSj1. Both relationship (63) and (64) are linear equations. We can use the method showed in A1 to obtain the least-norm solutions such as in the one-dimensional problem. Then, the new two-step MPDATA for two-dimensional advection equations is(65)ψj∗=ψjn−ΔtTj∑k=13Hψjn,ψjkn,vjk,njkhjk,ψjn+1=ψj∗−ΔtTj∑k=13Hψj∗,ψjk∗,v˜jkd,njkhjk,where(66)Hψjkn,ψjn,vjk,njk=ajkoutvjk+ψ¯jn+ajkinvjk−ψ¯jknajkin+ajkout⋅njk−ajkinajkoutajkin+ajkoutψ¯jkn−ψ¯jn.Next, we focus our attention on the positivity. Here, we have the following results for the one-step scheme (52).Theorem 2.
For finite-volume scheme (52), solving the advection system (49), if ψ¯jn≥0 for all j, then ψ¯jn+1≥0, under the conditions(67)Δt≤minj,kTjhjk4a∑k=13hjk,where a=maxj,kvjk+,vjk−.Proof.
Then, from scheme (52), we can have(68)ψjn+1=−ΔtTj∑k=13Hψjn,ψjkn,vjk,njkhjk=14ψjn+14M1+14M2+14M3,where(69)M1=ψjn−4ΔtTjHψjn,ψj1n,vj1,nj1hj1+Hψjn,ψjn,vj1+vj2+,−nj1hj1,(70)M2=ψjn−4ΔtTjHψjn,ψjn,vj1+vj2+,nj1hj1+Hψj2n,ψjn,vj2,nj2hj2+Hψjn,ψjn,vj3+vj2+,nj3hj3,(71)M3=ψjn−4ΔtTjHψjn,ψjn,vj3+vj2+,−nj3hj3+Hψj3n,ψjn,vj3,nj3hj3.
Then, we cxan say thatψjn+1 is a convex combination of ψjn, M1, M2, and M3, which implies if ψjn, M1, M2, and M3 are nonnegative, ψjn+1 is nonnegative. Noticing the identical relation,(72)∑k=13hjkHψjn,ψjn,vj2+vj2+,njk=∑k=13hjkajkoutvj2+ψ¯jn+ajkinvj2+ψ¯jnajkin+ajkout⋅njk−ajkinajkoutajkin+ajkoutψ¯jn−ψ¯jn=0.
Incorporating (72) into (70), we obtain(73)M2=hj1∑k=13jkM21+hj2∑k=13hjkM22+hj3∑k=13hjkM23,where(74)M21=ψjn−4Δt∑k=13hjkTjHψjn,ψjn,vj1+vj2+,nj1−Hψjn,ψjn,vj2+vj2+,nj1,M22=ψjn−4Δt∑k=13hjkTjHψj2n,ψjn,vj2,nj2−Hψjn,ψjn,vj2+vj2+,nj2,M23=ψjn−4Δt∑k=13hjkTjHψjn,ψjn,vj3+vj2+,nj3−Hψjn,ψjn,vj2+vj2+,nj3.
As eachM1, M3, M21, M22, and M23 can be viewed as a type of one-dimensional scheme, we can say that under the conditions(75)Δt≤minj,kTjhjk4a∑k=13hjk,for finite-volume scheme (49), solving the advection system (49) if ψ¯jn≥0 for all j, then ψ¯jn+1≥0.


## 2.1. Analysis of the Previous MPDATA

The irregular meshes discussed in this paper all refer to nonuniform irregular meshes. In this section, we analyze the reasons why previous MPDATA methods are not applicable to irregular grids. The one-dimensional advection equation is(1)∂∂tψ+∂∂xuψ=0,where ψ is a nondiffusive scalar quantity and u is the one-dimensional velocity. The former MPDATA scheme is proposed by Smolarkiewicz and Margolin [3], which is(2)ψi∗=ψin−Fψin,ψi+1n,ui+1/2n−Fψi−1n,ψin,ui−1/2n,(3)ψin+1=ψi∗−Fψi∗,ψi+1∗,u˜i+1/2−Fψi−1∗,ψi∗,u˜i−1/2,where(4)Fψin,ψi+1n,ui+1/2n=ui+1/2n+ui+1/2nψin+ui+1/2n−ui+1/2nψi+1nΔt2Δx.Here,ψin is the value of ψ at the ith frid point for the nth time step. u˜i+1/2 is the numerical approximation of ud at xi+1/2, and it is(5)u˜i+1/2=ui+1/2Δx−Δtui+1/22ψi+1∗+ψi∗+εψi+1∗−ψi∗Δx,where ε is a small value (e.g., 10−15) to ensure u˜=0 when ψi+1∗=ψi∗=0.For a regular grid, the numerical schemes (2) and (3) are an effective way to obtain positive definite solutions (see Section 3.1). However, neither (2) nor (3) are able to provide the second-order numerical accuracy solutions for irregular grids. There are two reasons for this problem. The first is the compatibility of numerical discretion of the “diffusion velocity” ud. The second is the compatibility of the first-order numerical scheme. To overcome the first problem, we have the following propositions.Proposition 1.
In (5), the numerical discretation of ∂ψ/∂xi+1/2 is(6)∂ψ∂xi+1/2=ψi+1∗−ψi∗Δx+OΔx,for regular grids, and discretation (6) is the first order. For an irregular grid, the discretation of (6) prevents step (3) from reaching first order. As a result, the MPDATA schemes (2) and (3) cannot reach the second order.Proof.
For the irregular grid, discretation (6) is(7)∂ψ∂xi+1/2=ψi+1∗−ψi∗Δx+C0,where C0 is a constant. This enables us to obtain the relationship between the numerical approximation values and the exact values:(8)ψi∗=ψxi+OΔx,(9)u˜=uxi+C0.
Substituting the abovementioned equations (8) and (9) into (4), we have(10)Fψi∗,ψi+1∗,u˜i+1/2=u˜xi+1/2+u˜xi+1/2ψxi+u˜xi+1/2−u˜xi+1/2ψxi+1+u˜xi+1/2OΔx+ψxiC1+ψxi+1C2+OΔx2Δt2Δx,where C1 and C2 are constant. Using the same method for ∂ψ/∂xi−1/2, we have(11)Fψi∗,ψi+1∗,u˜i+1/2−Fψi−1∗,ψi∗,u˜i−1/2=Fψxi,ψxi+1,u˜xi+1/2−Fψxi−1,ψxi,u˜xi−1/2+u˜xi+1/2−u˜xi−1/2Δt+C1′ψxi−1+C2′ψxi+C3′ψxi+1+OΔx2,where C1′, C2′, and C3′ are constant. When Δx⟶0,(12)C1′ψxi−1+C2′ψxi+C3′ψxi+1≠OΔx.
Thus, step (3) cannot become first order, and it cannot calculate the affect of diffusion effectively. Thus, the numerical schemes (2) and (3) cannot become the second order.
To obtain accurate numerical schemes on an irregular grid, a new approximation method of∂ψ/∂xi±1/2 is given here. The mesh structure around xi+1/2 is described in Figure 1.Figure 1
Schematic of the mesh structure aroundxi+1/2.Then, we expandψi−1, ψi, ψi+1,andψi+2 at the point ψi+1/2, and we have(13)ψs=ψi+1/2+xs−xi+1/2∂ψ∂xi+1/2+OΔx2.Here,s=i−1,i,i+1,i+2. A set of coefficients κ1i+1/2, i=1,2,3,4, need to be obtained, which satisfies(14)1111xi−1−xi+1/2xi−xi+1/2xi+1−xi+1/2xi+2−xi+1/2κ1i+1/2κ2i+1/2κ3i+1/2κ4i+1/2=01.It is easy to see that the rank of the coefficient matrix and the extended matrix for linear equations (11) are the same, and the linear equations (11) have infinite solutions. Here, the least-norm solutions are used (see Appendix). Then, we have the approximation for ∂ψ/∂x:(15)∂ψ∂xi+1/2=ψi−1⋅κ1i+1/2+ψiκ2i+1/2+ψi+1κ3i+1/2+ψi+2κ4i+1/2+OΔx.with the numerical approximation (), we have(16)ψi∗=ψxi+OΔx,u˜=uxi+OΔx.For the same reason as (11), we have(17)Fψi∗,ψi+1∗,u˜i+1/2−Fψi−1∗,ψi∗,u˜i−1/2=Fψxi,ψxi+1,u˜xi+1/2−Fψxi−1,ψxi,u˜xi−1/2+u˜xi+1/2−u˜xi−1/2Δt+C˜1ψxi−1Δt+C˜2ψxiΔt+C˜3ψxi+1Δt+OΔx2.Becauseu˜ is a continuous function, we have(18)C˜1ψxi−1Δt+C˜2ψxiΔt+C˜3ψxi+1Δt∼OΔxΔx⟶0.Thus, using (12) instead of (5) will enable step (3) to reach the first order, whereupon it would be able to calculate the affect of diffusion effectively. Furthermore, (2) and (3) can achieve second-order accuracy.Before we begin the discussion for a two-dimensional situation, we have another proposition to state.Proposition 2.
The compatibility of the first-order numerical scheme is important for the MPDATA method.
Generally speaking, the MPDATA schemes (2) and (3) belong to the class of staggered central-upwind schemes, which are excellent tools for solving various complex problems on regular domains. However, in practice complicated geometries exist, and the use of irregular triangular meshes could be advantageous or even unavoidable [7]. Staggered central-upwind schemes are sensitive to mesh effects. For example, Küther proved that Lax–Friedrichs-type finite-volume schemes estimate the error to be of the order of 0.25 [9, 10].
For the first-order scheme,(19)ψin+1=ψin−Fψin,ψi+1n,ui+1/2n−Fψi−1n,ψin,ui−1/2n.
With the Taylor expansionψi−1, ψi, ψi+1,andψi+2 at the point ψi+1/2, we have(20)Fψi,ψi+1,ui+1/2−Fψi−1,ψi,ui−1/2=ui+1/2−ui−1/2ψiΔtΔx+Δt4Δx∂ψ∂xiΔxi+Δxi+1ui+1/2−ui+1/2+Δxi−1+Δxiui−1/2+ui−1/2+OΔx2.
Whenu is a continuous constant, we have(21)Fψi,ψi+1,ui+1/2−Fψi−1,ψi,ui−1/2=u∂ψ∂xiΔxi−1+2Δxi+Δxi+14ΔxiΔt.
Based on the abovementioned equation, it is clear that if the mesh is distorted, a few points are incompatible. Whenu is not a discontinuous function, equation (21) implies that scheme (21) is largely incompatible.
We computed the irregular grid by carrying out a comparative simulation in Section 6.1. The results in Tables1–3 indicate that the standard MPDTA method reduces the order when calculated on an irregular grid.Table 1
Former MPDATA for regular grids.
N
L
1 error
L
1 order
L
2 error
L
2 order
L
n error
L
∞ order
16
4.333
E
−
02
—
3.330
E
−
02
—
3.313
E
−
02
—
32
1.115
E
−
02
1.96
8.513
E
−
03
1.97
9.022
E
−
03
1.87
64
2.717
E
−
03
2.04
2.132
E
−
03
2.00
2.545
E
−
03
1.82
128
6.648
E
−
04
2.03
5.280
E
−
04
2.01
6.412
E
−
04
1.98
256
1.640
E
−
04
2.02
1.308
E
−
04
2.01
1.596
E
−
04
2.00Table 2
Former MPDATA for irregular grids.
N
L
1 error
L
1 order
L
2 error
L
2 order
L
n error
L
∞ order
16
2.181
E
−
01
—
1.834
E
−
01
—
2.471
E
−
01
—
32
6.887
E
−
02
1.66
6.200
E
−
02
1.57
1.070
E
−
01
1.21
64
2.188
E
−
02
1.65
2.176
E
−
02
1.51
4.905
E
−
02
1.13
128
7.550
E
−
03
1.53
8.901
E
−
03
1.29
2.385
E
−
02
1.04
256
3.099
E
−
03
1.28
4.125
E
−
03
1.11
1.184
E
−
02
1.01Table 3
New MPDATA for irregular grids.
N
L
1 error
L
1 order
L
2 error
L
2 order
L
n error
L
∞ order
16
1.976
E
−
01
—
1.568
E
−
01
—
1.830
E
−
01
—
32
5.744
E
−
02
1.78
4.765
E
−
02
1.72
5.213
E
−
02
1.81
64
1.450
E
−
02
1.99
1.272
E
−
02
1.91
1.536
E
−
02
1.76
128
3.680
E
−
03
1.98
3.240
E
−
03
1.97
3.935
E
−
03
1.97
256
9.221
E
−
04
2.00
8.146
E
−
04
1.99
9.830
E
−
04
2.00


## 2.2. New MPDATA for Irregular Grids in One-Dimensional Problems

Base on the abovementioned analysis, we propose a new MPDATA method for irregular grids in the one-dimensional situation. One of the grid elements being studied is denoted asTj. After time step Δt, the discontinuity around xj±1/2 can be described as ajLinΔt and ajRinΔt for inward Tj, where ajLinΔt is around xj−1/2 and ajRinΔt is around xj+1/2. For outward, ajLout is around xj−1/2 and ajRout is around xj+1/2. Using this information, the control volume can be constructed as follows. As shown in Figure 2, the inside domain Dj; the boundary domain Dj±1/2, which can be divided into the inward domain Dj±1/2in (Dj−1/2in=ajLinΔt and Dj+1/2in=ajRinΔt) and the outward domain Dj±1/2out (Dj−1/2out=ajLoutΔt and Dj+1/2out=ajRoutΔt). All domains around element Tj are described as illustrated in Figure 2.Figure 2
Domains on control volumeTj.LetujL be the velocity in domain Dj−1/2in and ujR be the velocity in domain Dj+1/2in. For the advection equations (1), the discontinuity propagate with inward and outward Tj local speeds, which can be estimated by(22)ajLin=−min−uj−1R,−ujL,0,ajRin=−minuj+1L,ujR,0,(23)ajLout=max−uj−1R,−ujL,0,ajRout=maxuj+1L,ujR,0.Then, we can compute the averages over the corresponding domains. The domain average forDj×tn,tn+1 is(24)ϕ¯jn+1=1Dj∫Djψjnxdx−∫tntn+1uψxj−1/2+ajRinΔtxj+1/2−ajRinΔtdt=ψ¯jn−1Dj∫tntn+1uψxj−1/2+ajRinΔtxj+1/2−ajRinΔtdt+OΔx,and the domain average for Dj+1/2×tn,tn+1. is(25)ψ¯j+1/2n+1=1Dj+1/2∫Dj+1/2inψjnxdx+∫Dj+1/2outψj+1nxdx−∫tntn+1uj+1Lψj+1−ujRψjdt=1Dj+1/2Dj+1/2inψ¯jn+Dj+1/2outψ¯jn+1−∫tntn+1uj+1Lψj+1−ujRψjdt+OΔx.On the other hand, we can have the relationship(26)ψ¯jn=1Tj∫Tjψx,tndx=1Tj∫Tjϕnxdx.The time derivative ofψjt is expressed with the help of (24) as(27)ddtψ¯jt=limΔt⟶0ψ¯jn+1−ψ¯jnΔt=limΔt⟶01Δt1Tj∫Tjϕn+1xdx−ψ¯jn.Because the width of the Riemann fans is bounded byajLout+ajLinΔt or ajRout+ajRinΔt, we obtain(28)ϕj±1/2n+1x=ϕ¯j±1/2n+1+OΔt,where x∈Dj. From (28) and the conservation of ϕjn+1, we derive(29)ddtψ¯jt=limΔt⟶01Δt1Tj∫Dj−1/2inϕn+1xdx+limΔt⟶01Δt1Tj∫Djϕn+1xdx+limΔt⟶01Δt1Tj∫Dj+1/2inϕn+1xdx−limΔt⟶01Δtψ¯jn=ajLinTjlimΔt⟶0ϕ¯j−1/2n+1+limΔt⟶01ΔtDjTjϕ¯jn+1−ψ¯jn+ajRinTjlimΔt⟶0ϕ¯j+1/2n+1.With the help of (24) and (25), we obtain(30)limΔt⟶01ΔtDjTjϕ¯jn+1−ψ¯jn=−ajRinψjn+ajLinψjnTj−ujR−ujLψjnTj+OΔx,and(31)limΔt⟶0ϕ¯j+1/2n+1=ajRinψ¯j+ajRoutψ¯j+1ajRin+ajRout−−ujRψ¯j+uj+1Lψ¯j+1ajRin+ajRout+OΔx.With the help of the forward Euler method, the substitution of (30) and (31) in (29) results in the new first-order scheme:(32)ψ¯jn+1=ψ¯jn−ΔtΔxjHj+1/2ψ¯j,ψ¯j+1,u−Hj−1/2ψ¯j−1,ψ¯j,u,where(33)Hj+1/2=ajRoutujRψ¯jn+ajRinuj+1Lψ¯j+1najRout+ajRin−ajRoutajRinajRout+ajRinψ¯j+1n−ψ¯jn,Hj−1/2=ajLoutujLψ¯jn+ajLinuj−1Rψ¯j−1najLout+ajLin+ajLoutajLinajLout+ajLinψ¯j−1n−ψ¯jn.Using the Taylor expansionxj+1, xj at xj+1/2(34)ψj+1n=ψj+1/2n+12Δxj+1∂∂xψj+1/2n+OΔx2,(35)ψjn=ψj+1/2n−12Δxj∂∂xψj+1/2n+OΔx2,and xj, xj−1 at xj−1/2:(36)ψj−1n=ψj−1/2n−12Δxj−1∂∂xψj−1/2n+OΔx2,(37)ψjn=ψj−1/2n+12Δxj∂∂xψj−1/2n+OΔx2.Taking those Taylor expressions (34)–(37) into Hj+1/2 and Hj−1/2, we can obtain(38)Hj+1/2=ajRoutujR+ajRinuj+1LajRout+ajRinψj+1/2n+−ajRoutujRΔxj+ajRinuj+1LΔxj+1ajRout+ajRin12−Δxj+1+Δxj2ajRoutajRinajRout+ajRin∂∂xψj+1/2n+OΔx2,Hj−1/2=ajLoutujL+ajLinuj−1RajLout+ajLinψj−1/2n+ajLoutujLΔxj−ajLinuj−1RΔxj−1ajLout+ajLin12−Δxj−1+Δxj2ajLoutajLinajLout+ajLin∂∂xψj−1/2n+OΔx2.With the Taylor expansions on time,(39)ψjn+1=ψjn+1/2+12Δt∂ψ∂tjn+1/2+OΔt2,ψjn=ψjn+1/2−12Δt∂ψ∂tjn+1/2+OΔt2=ψjn+1/2+12Δt∂uψ∂xjn+1/2+OΔt2.Scheme (32) becomes(40)∂ψ∂tjn+1/2=1ΔxjajRoutujR+ajRinuj+1LajRout+ajRinψj+1/2n+1/2−ajLoutujL+ajLinuj−1RajLout+ajLinψj−1/2n+1/2+1ΔxjE˜,where(41)E˜:=u˜j+1/2ψj+1/2n+1/2−u˜j−1/2ψj−1/2n+1/2,u˜j+1/2:=1ψj+1/2n+1/2ajRoutujR+ajRinuj+1LajRout+ajRin12Δt∂∂xuψj+1/2n+1/2+1ψj+1/2n+1/2−ajRoutujRΔxj+ajRinuj+1LΔxj+1ajRout+ajRin12−Δxj+1+Δxj2ajRoutajRinajRout+ajRin∂∂xψj+1/2n,u˜j−1/2:=1ψj−1/2n+1/2ajLoutujL+ajLinuj−1RajLout+ajLin12Δt∂∂xuψj−1/2n+1/2+1ψj−1/2n+1/2ajLoutujLΔxj−ajLinuj−1RΔxj−1ajLout+ajLin12+Δxj−1+Δxj2ajLoutajLinajLout+ajLin∂∂xψj−1/2n.Whenu is continuous around xj−1/2 and xj+1/2, noticing one of aout or ain is zero. Then,(42)ajRoutujR+ajRinuj+1LajRout+ajRin⟹uj+1/2,ajLoutujL+ajLinuj−1RajLout+ajLin⟹uj−1/2.Thus, equation (40) can be viewed as(43)∂ψ∂t=∂uψ∂x+∂u˜ψ∂x.Here, the second term of (43) is the “pseudoadvection term.” (32) is the second-order scheme of equation (43). Therefore, we have the new MPDATA for irregular grids:(44)ψj∗=ψjn−Hj+1/2ψjn,ψj+1n,u−Hj−1/2ψj−1n,ψjn,uΔtΔxj,(45)ψjn+1=ψj∗−Hj+1/2ψj∗,ψj+1∗,u˜d−Hj−1/2ψj−1∗,ψj∗,u˜dΔtΔxj.Because the new MPDATA scheme depends on the two-step scheme, we can analyze the positivity of (32) instead of the schemes (44) and (45). Thus, we have the following claim for the new MPDATA schemes (44) and (45).Theorem 1.
Consider system (1) and the central-upwind discrete scheme (29), where ajLin, ajRin, ajLout, and ajRout are defined in (22) and (23). Then, for all j, ψ¯jn+1≥0, provided that Δt≤Δx/a, where Δx:=minjΔxj and a:=maxjmaxuj+1L,ujR.Proof.
For convenience of expression, let(46)λj=ΔtΔxj.
From (32), we can have(47)ψ¯jn+1=ψ¯jn−λjHj+1/2ψ¯jn,ψ¯j+1n,u−Hj−1/2ψ¯j−1n,ψ¯jn,u=12−λjajRoutajRin+ujRajRin+ajRoutψ¯jn+12+λjajLout−ajLin+ujLajLin+ajLoutψ¯jn+λjajRin−uj+1L+ajRoutajRin+ajRoutψ¯j+1n+λjajLinuj−1R+ajLoutajLin+ajLoutψ¯j−1n.
In the above equation (47), under ψ¯jn≥0, it is easy to see that(48)12+λjajLout−ajLin+ujLajLin+ajLoutψ¯jn≥0,λjajRin−uj+1L+ajRoutajRin+ajRoutψ¯j+1n≥0,λjajLinuj−1R+ajLoutajLin+ajLoutψ¯j−1n≥0,12−λjajRoutajRin+ujRajRin+ajRoutψ¯jn≥0,whenλj≤12ajRout.
Thus, we can claim that, for allj, ψ¯jn+1≥0, provided that Δt≤Δx/a, where Δx:=minjΔxj and a:=maxjmaxuj+1L,ujR.


## 2.3. New MPDATA for Irregular Grids in a Two-Dimensional Problem

For the two-dimensional advection problem,(49)∂ψ∂t=−∇⋅vψ,where ψx,t is a nondiffusive scalar field assumed to be nonnegative at t=0, and the v=vx is the prescribed flow. The finite-volume MPDATA is considered on equation (48). A triangular partition ∪jTj is used on the computational domain. The area of Tj is Tj. All domains around the triangular cell Tj are illustrated in Figure 3. The neighborhood cell near Tj is denoted as Tjk, k=1,2,3. The three edges of cell Tj are denoted as Sjk. The discontinuous long edges Sjk can move inward, in which case Tj is denoted as ajkinΔt, and outward such that Tj is represented as ajkoutΔt. The inside domain is Dj. The three nonconvex domains surrounding the corners are Ejk, k=1,2,3.Figure 3
Triangular mesh of central-upwind scheme.Similar to the one-dimensional problem,ajkin and ajkout can be interpreted as the possible discontinuities of ψj along the sides of Tj propagating with different inward, ajkin, and outward, ajkout, directional local speeds k=1,2,3. They can be estimated as(50)ajkin=−minminM∈Tj∩Tjkλ1VjkTj,minM∈Tj∩Tjkλ1VjkTjk,0,ajkout=maxmaxM∈Tj∩TjkλNVjkTj,maxM∈Tj∩TjkλNVjkTj,0,where λ1Vjk≤⋯≤λNVjk are the N eigenvalues of the matrix Vjk defined by(51)VjkT=njk⋅vjkT,and njk is the outer unit normals to the corresponding sides of Tj of length hjk, k=1,2,3. Following the similar method in Section 2 and the theory developed by Kurganov and Petrova [7], we obtain a consistent, accurate, and effective well-balanced numerical scheme:(52)ψ¯jn+1−ψ¯jnΔt=−1Tj∑k=13hjkajkinvjk−ψ¯jkn+ajkoutvjk+ψ¯jnajkin+ajkout⋅njk+1Tj∑k=13hjkajkinajkoutajkin+ajkoutψ¯jkn−ψ¯jn,where vjk+ and vjk− are advection velocities in domain Djk+ and domain Djk−. Here, (52) is a first-order scheme. Then, expanding ψj and ψjk at the middle point of the edge intersecting face Sj,(53)ψ¯jk=ψSj+ϑjkψSj+OΔx2+Δy2,where ϑjk=xjk−xSj∂/∂x+yjk−ySj∂/∂y.(54)ψ¯j=ψSj+ϑjψSj+OΔx2+Δy2,where ϑj=xj−xSj∂/∂x+yj−ySj∂/∂y. We substitute (53) and (54) into (52), and we have(55)ψ¯jn+1−ψ¯jnΔt=−1Tj∑k=13hjkajkinvjk−ψSjn+ajkoutvjk+ψSjnajkin+ajkout⋅njk−1Tj∑k=13hjkajkinvjk−ϑjkψSjn+ajkoutvjk+ϑjψSjnajkin+ajkout⋅njk+1Tj∑k=13hjkajkinajkoutajkin+ajkoutϑjkψSjn−ϑjψSjn.Expandingψn+1 and ψn at ψn+1/2,(56)ψn+1=ψn+1/2+12∂ψ∂tn+1/2Δt+OΔt2,ψn=ψn+1/2−12∂ψ∂tn+1/2Δt+OΔt2=ψn+1/2+12Δtψ∇vn+1/2+12Δtv∇ψn+1/2+OΔt2.Then, we have the relationship(57)∂ψ¯∂tn+1/2=−1Tj∑k=13hjkajkinvjk−+ajkoutvjk+ajkin+ajkoutψSjn+1/2⋅njk+Error,where(58)Error=−1Tj∑k=13hjkajkinvjk−+ajkoutvjk+ajkin+ajkout12Δtψ∇vSjn+1/2⋅njk−1Tj∑k=13hjkajkinvjk−+ajkoutvjk+ajkin+ajkout12Δtv∇ψSjn+1/2⋅njk−1Tj∑k=13hjkajkinvjkϑjkψSjn+ajkoutvjϑjψSjnajkin+ajkout⋅njk+1Tj∑k=13hjkajkinajkoutajkin+ajkout⋅njkϑjkψSjn−ϑjψSjn⋅njk.Noticing that(59)ajkinvjk−+ajkoutvjk+ajkin+ajkout⟹vSj,(59) can be can be understood as the average value of v in the small rectangle Djk (see Figure 3). Furthermore, we can assume t=nΔt, and (57) can be written as(60)∂ψ∂t=−∇⋅vψ+∇⋅v˜dψ,where(61)v˜jkd:=−ajkinvjk−+ajkoutvjk+ajkin+ajkout12Δt∇⋅vSjn+1/2−1ψajkinvjk−+ajkoutvjk+ajkin+ajkout12Δtv⋅∇ψSjn+1/2−1ψajkinvjk−ϑjkψSjn+ajkoutvjk+ϑjψSjnajkin+ajkout+1ψajkinajkoutajkin+ajkout⋅njkϑjkψSjn−ϑjψSjn.Here, the gradient∇ψ at the edges can be obtained by the combination of the physical values around the edge. The mesh structure, which is used to obtain the gradient ∇ψj1, around the midpoint of an edge Sj1 is described in Figure 4.Figure 4
Schematic of the mesh structure aroundSj1.Here, we useψj1 to represent the physical value at point Mj1. Using the Taylor expression, expand ψ¯j, ψ¯j1, ψ¯j2, ψ¯j3, ψ¯j11, and ψ¯j12 at point Mj1 to second-order, and we have(62)ψ¯j=ψSj1+xj−xSj1∂∂xψSj1+yj−ySj1∂∂yψSj1+OΔx2+Δy2,ψ¯j1=ψSj1+xj1−xSj1∂∂xψSj1+yj1−ySj1∂∂yψSj1+OΔx2+Δy2,ψ¯j2=ψSj1+xj2−xSj1∂∂xψSj1+yj2−ySj1∂∂yψSj1+OΔx2+Δy2,ψ¯j3=ψSj1+xj3−xSj1∂∂xψSj1+yj3−ySj1∂∂yψSj1+OΔx2+Δy2,ψ¯j11=ψSj1+xj11−xSj1∂∂xψSj1+yj11−ySj1∂∂yψSj1+OΔx2+Δy2,ψ¯j12=ψSj1+xj12−xSj1∂∂xψSj1+yj12−ySj1∂∂yψSj1+OΔx2+Δy2.Next, use the relationship(63)∑i=05αiφSi=0,∑i=05αiyi−ySj1=0,∑i=05αixj1−xSj1=1,to obtain ∂/∂xφSj1, and relationship(64)∑i=05αiφSi=0,∑i=05αiyi−ySj1=0,∑i=05αixj1−xSj1=1,for ∂/∂yφSj1. Both relationship (63) and (64) are linear equations. We can use the method showed in A1 to obtain the least-norm solutions such as in the one-dimensional problem. Then, the new two-step MPDATA for two-dimensional advection equations is(65)ψj∗=ψjn−ΔtTj∑k=13Hψjn,ψjkn,vjk,njkhjk,ψjn+1=ψj∗−ΔtTj∑k=13Hψj∗,ψjk∗,v˜jkd,njkhjk,where(66)Hψjkn,ψjn,vjk,njk=ajkoutvjk+ψ¯jn+ajkinvjk−ψ¯jknajkin+ajkout⋅njk−ajkinajkoutajkin+ajkoutψ¯jkn−ψ¯jn.Next, we focus our attention on the positivity. Here, we have the following results for the one-step scheme (52).Theorem 2.
For finite-volume scheme (52), solving the advection system (49), if ψ¯jn≥0 for all j, then ψ¯jn+1≥0, under the conditions(67)Δt≤minj,kTjhjk4a∑k=13hjk,where a=maxj,kvjk+,vjk−.Proof.
Then, from scheme (52), we can have(68)ψjn+1=−ΔtTj∑k=13Hψjn,ψjkn,vjk,njkhjk=14ψjn+14M1+14M2+14M3,where(69)M1=ψjn−4ΔtTjHψjn,ψj1n,vj1,nj1hj1+Hψjn,ψjn,vj1+vj2+,−nj1hj1,(70)M2=ψjn−4ΔtTjHψjn,ψjn,vj1+vj2+,nj1hj1+Hψj2n,ψjn,vj2,nj2hj2+Hψjn,ψjn,vj3+vj2+,nj3hj3,(71)M3=ψjn−4ΔtTjHψjn,ψjn,vj3+vj2+,−nj3hj3+Hψj3n,ψjn,vj3,nj3hj3.
Then, we cxan say thatψjn+1 is a convex combination of ψjn, M1, M2, and M3, which implies if ψjn, M1, M2, and M3 are nonnegative, ψjn+1 is nonnegative. Noticing the identical relation,(72)∑k=13hjkHψjn,ψjn,vj2+vj2+,njk=∑k=13hjkajkoutvj2+ψ¯jn+ajkinvj2+ψ¯jnajkin+ajkout⋅njk−ajkinajkoutajkin+ajkoutψ¯jn−ψ¯jn=0.
Incorporating (72) into (70), we obtain(73)M2=hj1∑k=13jkM21+hj2∑k=13hjkM22+hj3∑k=13hjkM23,where(74)M21=ψjn−4Δt∑k=13hjkTjHψjn,ψjn,vj1+vj2+,nj1−Hψjn,ψjn,vj2+vj2+,nj1,M22=ψjn−4Δt∑k=13hjkTjHψj2n,ψjn,vj2,nj2−Hψjn,ψjn,vj2+vj2+,nj2,M23=ψjn−4Δt∑k=13hjkTjHψjn,ψjn,vj3+vj2+,nj3−Hψjn,ψjn,vj2+vj2+,nj3.
As eachM1, M3, M21, M22, and M23 can be viewed as a type of one-dimensional scheme, we can say that under the conditions(75)Δt≤minj,kTjhjk4a∑k=13hjk,for finite-volume scheme (49), solving the advection system (49) if ψ¯jn≥0 for all j, then ψ¯jn+1≥0.


## 3. Numerical Results

### 3.1. One-Dimensional Accuracy Test

We consider the equation that expresses the linear scalar conservation law:(76)∂∂tψ+∂∂xψ=0.Here, the periodic boundary condition is used, and the exact solution is(77)ψx,t=1+0.5sinπx−t.The computational domain is0,2, and the final time T=2.0. The initial condition is ux,0. To solve this problem, we use three approaches for simulation for the purpose of comparison. First, the simulation is conducted on regular grids with the aforementioned former MPDATA proposed by Smolarkiewicz and Margolin [3]. Second, the same algorithm, MPDATA by Smolarkiewicz, is used on irregular grids. Lastly, the simulation is carried out on irregular grids with the new MPDATA method. The construction of these irregular grids is illustrated in Figure 5. Here, the left-most cell size is assumed to be Δx, in which case the second left-most cell size is 2Δx, the third left-most cell size is 3Δx, and so on. In Figure 6, the convergence order of the former MPDATA (Figure 6(a)) and the new MPDATA (Figure 6(b)) is compared with different number of grids. For the sake of distinction, a set of parallel dashed lines with a slope of −2 is drawn in Figure 6. It can be seen from the figure that the new MPDATA can achieve second-order convergence, while the former scheme cannot reach the second-order accuracy. The simulation results of all three methods with multiple group cells are presented in Tables 1–3. N is the number of cells in the simulation. Table 1 contains the results of the simulation of the former MPDATA for regular grids, indicating that the former MPDATA scheme for regular grid simulation is a useful tool and that L1, L2, and L∞ tend to the second-order norm as the mesh increases. The simulation results obtained with the former MPDATA on irregular grids are listed in Table 2. A comparison of Tables 1 and 2 reveals that the former MPDATA is useful for regular grids but that it fails for irregular grids. Table 3 is devoted to the simulation results of the new MPDATA on irregular grids, showing that L1, L2, and L∞ tend to second-order as the mesh increases.Figure 5
Schematic of construction of irregular grids.Figure 6
MPDATA simulation with 64 cells. (a) Former MPDATA with regular grids; (b) new MPDATA with irregular grid.
(a)
(b)


### 3.2. Two-Dimensional Accuracy Test

The two-dimensional accuracy test is taken on the linear convection equation:(78)∂ψ∂t+v0⋅∇ψ=0,v0=0.5,0.5.The periodic boundary condition is used as the one-dimensional situation. The exact solution is taken as(79)ψx,y,t=1+0.2sinx+y−t.The computational domain is0,2π×0,2π. The computational time is T=2.0. The periodic boundary condition is used. The computational domain with 32×32 is shown in Figure 7 and is surrounded by ghost cells. The simulation results are shown in Figure 8. Figure 8(a) shows the simulation result with 128×128 cells, which implies that our scheme can accommodate the situation of an unstructured grid well. The results in Figure 8(b) were simulated at a different resolution. The results confirm that our schemes can attain second-order accuracy as the mesh increases. In Table 4, we give the test results of the accuracy of the new MPDATA for two-dimensional situation. It can be seen that the new MPDATA can reached second order.Figure 7
Computational domain with32×32 mesh.Figure 8
Simulation results. (a) Simulation with128×128 mesh. (b) Simulation at a different resolution.
(a)
(b)Table 4
New MPDATA for irregular grids.
N
L
1 error
L
1 order
L
2 error
L
2 order
L
n error
L
∞ order
16
×
16
6.623
E
−
01
—
1.242
E
−
01
—
3.892
E
−
02
—
32
×
32
1.968
E
−
01
1.75
3.799
E
−
02
1.70
1.657
E
−
02
1.23
64
×
64
4.870
E
−
02
2.01
9.651
E
−
03
1.97
5.945
E
−
03
1.47
128
×
128
1.152
E
−
02
2.07
2.255
E
−
03
2.09
1.504
E
−
03
1.98
256
×
256
2.721
E
−
03
2.08
5.382
E
−
04
2.06
3.735
E
−
04
2.01


### 3.3. Cosine Bell Test for Two-Dimensional Situations

The cosine bell test was suggested by Huang et al. [11]. The cosine bell with an unsmooth boundary and all nonnegative values was used to examine the spurious oscillations as well as the damping effects of both schemes in the 2D case. The initial distribution of the cosine bell is given by(80)ψx,y,0=121.0+cosπrr0,r≤r0,0,r>r0,where r=x−x02+y−y02, r=0.5, and x0=y0. Irregular grids were used in the simulations. The initial values are synoptically shown in Figure 9 for a 32×32 area.Figure 9
Initial computational area.The finial time of the simulation isT=2. In our test, the cosine bell passes through areas of dense grid and areas of sparse grid. A mesh of 300×300 is used in the simulation in Figure 10. Figure 10(a) shows the results at time T=1, whereas those in Figure 10(b) were obtained at time T=2. These results enabled us to conclude that our schemes are able to solve this problem without obtaining negative values.Figure 10
Computational results at time (a)T=1; (b)T=2.
(a)
(b)


## 3.1. One-Dimensional Accuracy Test

We consider the equation that expresses the linear scalar conservation law:(76)∂∂tψ+∂∂xψ=0.Here, the periodic boundary condition is used, and the exact solution is(77)ψx,t=1+0.5sinπx−t.The computational domain is0,2, and the final time T=2.0. The initial condition is ux,0. To solve this problem, we use three approaches for simulation for the purpose of comparison. First, the simulation is conducted on regular grids with the aforementioned former MPDATA proposed by Smolarkiewicz and Margolin [3]. Second, the same algorithm, MPDATA by Smolarkiewicz, is used on irregular grids. Lastly, the simulation is carried out on irregular grids with the new MPDATA method. The construction of these irregular grids is illustrated in Figure 5. Here, the left-most cell size is assumed to be Δx, in which case the second left-most cell size is 2Δx, the third left-most cell size is 3Δx, and so on. In Figure 6, the convergence order of the former MPDATA (Figure 6(a)) and the new MPDATA (Figure 6(b)) is compared with different number of grids. For the sake of distinction, a set of parallel dashed lines with a slope of −2 is drawn in Figure 6. It can be seen from the figure that the new MPDATA can achieve second-order convergence, while the former scheme cannot reach the second-order accuracy. The simulation results of all three methods with multiple group cells are presented in Tables 1–3. N is the number of cells in the simulation. Table 1 contains the results of the simulation of the former MPDATA for regular grids, indicating that the former MPDATA scheme for regular grid simulation is a useful tool and that L1, L2, and L∞ tend to the second-order norm as the mesh increases. The simulation results obtained with the former MPDATA on irregular grids are listed in Table 2. A comparison of Tables 1 and 2 reveals that the former MPDATA is useful for regular grids but that it fails for irregular grids. Table 3 is devoted to the simulation results of the new MPDATA on irregular grids, showing that L1, L2, and L∞ tend to second-order as the mesh increases.Figure 5
Schematic of construction of irregular grids.Figure 6
MPDATA simulation with 64 cells. (a) Former MPDATA with regular grids; (b) new MPDATA with irregular grid.
(a)
(b)


## 3.2. Two-Dimensional Accuracy Test

The two-dimensional accuracy test is taken on the linear convection equation:(78)∂ψ∂t+v0⋅∇ψ=0,v0=0.5,0.5.The periodic boundary condition is used as the one-dimensional situation. The exact solution is taken as(79)ψx,y,t=1+0.2sinx+y−t.The computational domain is0,2π×0,2π. The computational time is T=2.0. The periodic boundary condition is used. The computational domain with 32×32 is shown in Figure 7 and is surrounded by ghost cells. The simulation results are shown in Figure 8. Figure 8(a) shows the simulation result with 128×128 cells, which implies that our scheme can accommodate the situation of an unstructured grid well. The results in Figure 8(b) were simulated at a different resolution. The results confirm that our schemes can attain second-order accuracy as the mesh increases. In Table 4, we give the test results of the accuracy of the new MPDATA for two-dimensional situation. It can be seen that the new MPDATA can reached second order.Figure 7
Computational domain with32×32 mesh.Figure 8
Simulation results. (a) Simulation with128×128 mesh. (b) Simulation at a different resolution.
(a)
(b)Table 4
New MPDATA for irregular grids.
N
L
1 error
L
1 order
L
2 error
L
2 order
L
n error
L
∞ order
16
×
16
6.623
E
−
01
—
1.242
E
−
01
—
3.892
E
−
02
—
32
×
32
1.968
E
−
01
1.75
3.799
E
−
02
1.70
1.657
E
−
02
1.23
64
×
64
4.870
E
−
02
2.01
9.651
E
−
03
1.97
5.945
E
−
03
1.47
128
×
128
1.152
E
−
02
2.07
2.255
E
−
03
2.09
1.504
E
−
03
1.98
256
×
256
2.721
E
−
03
2.08
5.382
E
−
04
2.06
3.735
E
−
04
2.01


## 3.3. Cosine Bell Test for Two-Dimensional Situations

The cosine bell test was suggested by Huang et al. [11]. The cosine bell with an unsmooth boundary and all nonnegative values was used to examine the spurious oscillations as well as the damping effects of both schemes in the 2D case. The initial distribution of the cosine bell is given by(80)ψx,y,0=121.0+cosπrr0,r≤r0,0,r>r0,where r=x−x02+y−y02, r=0.5, and x0=y0. Irregular grids were used in the simulations. The initial values are synoptically shown in Figure 9 for a 32×32 area.Figure 9
Initial computational area.The finial time of the simulation isT=2. In our test, the cosine bell passes through areas of dense grid and areas of sparse grid. A mesh of 300×300 is used in the simulation in Figure 10. Figure 10(a) shows the results at time T=1, whereas those in Figure 10(b) were obtained at time T=2. These results enabled us to conclude that our schemes are able to solve this problem without obtaining negative values.Figure 10
Computational results at time (a)T=1; (b)T=2.
(a)
(b)


## 4. Conclusion and Discussion

We combined the MPDATA method with the semidiscrete central difference strategy. We presented a positive definite shape-preserving scheme for complex grid structures. The MPDATA method proposed by Smolarkiewicz et al. is a useful tool for atmospheric dynamics simulation. However, their method is problematic as it experiences loss of precision when computing on an irregular grid. The reasons for this problem are two fold: the first is the numerical discretization of the boundary derivative in the finite-volume method; the other is that the numerical format is developed based on a staggered grid. Our approach overcomes these two problems in combination with the multipoint Taylor expansion method. We thus obtain the boundary derivative approximation method that does not depend on the grid structure. Combined with the central-upwind scheme, we propose a positive definite advection scheme for complex grid structures. On this basis, we carried out the positive definite analysis and numerical test of the new numerical scheme.In this paper, starting from the one-dimensional MPDATA [3], we found that the existing MPDATA method [3–5] could not reach the second-order calculation accuracy required by the design for nonuniform grids, and then through improvement, we obtained the MPDATA method that could reach the second-order accuracy for the irregular grid. Here, we need to explain that the irregular grid studied in this paper is different from the unstructured grid. The irregular grid we studied specifically refers to the grid with the nonuniform grid structure. Traditional MPDATA unstructured mesh algorithms can achieve two precision for uniform unstructured grids but cannot achieve second precision for nonuniform irregular grids.In addition, in order to explain that the existing MPDATA methods cannot achieve second-order accuracy for nonuniform irregular grids, the one-dimensional case of traditional MPDATA methods is mainly discussed here. However, in the case of two-dimensional MPDATA method, there are similar accuracy problems. The most direct point is in the calculation of boundary derivatives. The unified MPDTA scheme still adopts the calculation of the difference between physical quantities of adjacent grids which is similar to one-dimensional MPDATA [4, 5]:(81)∂Φ∂xIj=1Vj+Vi∑m=1liΦ¯i,mSmI+∑m′=1liΦ¯i,m′Sm′I,where Φ¯i,m=0.5Φi+Φm.



---

*Source: 5283139-2020-04-10.xml*