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.
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 Section 2. 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 Section 2, 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 iswhere is a nondiffusive scalar quantity and is the one-dimensional velocity. The former MPDATA scheme is proposed by Smolarkiewicz and Margolin [3], which iswhere
Here, is the value of at the th frid point for the th time step. is the numerical approximation of at , and it iswhere is a small value (e.g., ) to ensure when .
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” . 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 isfor 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) iswhere is a constant. This enables us to obtain the relationship between the numerical approximation values and the exact values:Substituting the abovementioned equations (8) and (9) into (4), we havewhere and are constant. Using the same method for , we havewhere , , and are constant. When ,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 is given here. The mesh structure around is described in Figure 1.
Then, we expand , , at the point , and we have
Here, . A set of coefficients , , need to be obtained, which satisfies
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 :with the numerical approximation (), we have
For the same reason as (11), we have
Because is a continuous function, we have
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,With the Taylor expansion , , at the point , we haveWhen is a continuous constant, we haveBased on the abovementioned equation, it is clear that if the mesh is distorted, a few points are incompatible. When 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 Tables 1–3 indicate that the standard MPDTA method reduces the order when calculated on an irregular grid.
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 as . After time step , the discontinuity around can be described as and for inward , where is around and is around . For outward, is around and is around . Using this information, the control volume can be constructed as follows. As shown in Figure 2, the inside domain ; the boundary domain , which can be divided into the inward domain ( and ) and the outward domain ( and ). All domains around element are described as illustrated in Figure 2.
Let be the velocity in domain and be the velocity in domain . For the advection equations (1), the discontinuity propagate with inward and outward local speeds, which can be estimated by
Then, we can compute the averages over the corresponding domains. The domain average for isand the domain average for is
On the other hand, we can have the relationship
The time derivative of is expressed with the help of (24) as
Because the width of the Riemann fans is bounded by or , we obtainwhere . From (28) and the conservation of , we derive
With the help of (24) and (25), we obtainand
With the help of the forward Euler method, the substitution of (30) and (31) in (29) results in the new first-order scheme:where
Using the Taylor expansion , at and , at :
Taking those Taylor expressions (34)–(37) into and , we can obtainWith the Taylor expansions on time,
Scheme (32) becomeswhere
When is continuous around and , noticing one of or is zero. Then,
Thus, equation (40) can be viewed as
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:
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 , , , and are defined in (22) and (23). Then, for all , , provided that , where and .
Proof. For convenience of expression, letFrom (32), we can haveIn the above equation (47), under , it is easy to see thatThus, we can claim that, for all j, , provided that , where and .
2.3. New MPDATA for Irregular Grids in a Two-Dimensional Problem
For the two-dimensional advection problem,where is a nondiffusive scalar field assumed to be nonnegative at , and the is the prescribed flow. The finite-volume MPDATA is considered on equation (48). A triangular partition is used on the computational domain. The area of is . All domains around the triangular cell are illustrated in Figure 3. The neighborhood cell near is denoted as , . The three edges of cell are denoted as . The discontinuous long edges can move inward, in which case is denoted as , and outward such that is represented as . The inside domain is . The three nonconvex domains surrounding the corners are , .
Similar to the one-dimensional problem, and can be interpreted as the possible discontinuities of along the sides of propagating with different inward, , and outward, , directional local speeds . They can be estimated aswhere are the eigenvalues of the matrix defined byand is the outer unit normals to the corresponding sides of of length , . 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:where and are advection velocities in domain and domain . Here, (52) is a first-order scheme. Then, expanding and at the middle point of the edge intersecting face ,where .where . We substitute (53) and (54) into (52), and we have
Expanding and at ,
Then, we have the relationshipwhere
Noticing that(59) can be can be understood as the average value of in the small rectangle (see Figure 3). Furthermore, we can assume , and (57) can be written aswhere
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 , around the midpoint of an edge is described in Figure 4.
Here, we use to represent the physical value at point . Using the Taylor expression, expand , , , , , and at point to second-order, and we have
Next, use the relationshipto obtain , and relationshipfor . 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 iswhere
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 for all , then , under the conditionswhere .
Proof. Then, from scheme (52), we can havewhereThen, we cxan say that is a convex combination of , , , and , which implies if , , , and are nonnegative, is nonnegative. Noticing the identical relation,Incorporating (72) into (70), we obtainwhereAs each , , , , and can be viewed as a type of one-dimensional scheme, we can say that under the conditionsfor finite-volume scheme (49), solving the advection system (49) if for all , then .
3. Numerical Results
3.1. One-Dimensional Accuracy Test
We consider the equation that expresses the linear scalar conservation law:
Here, the periodic boundary condition is used, and the exact solution is
The computational domain is , and the final time . The initial condition is . 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 , in which case the second left-most cell size is , the third left-most cell size is , 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 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. 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 , , and 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 , , and tend to second-order as the mesh increases.
(a)
(b)
3.2. Two-Dimensional Accuracy Test
The two-dimensional accuracy test is taken on the linear convection equation:
The periodic boundary condition is used as the one-dimensional situation. The exact solution is taken as
The computational domain is . The computational time is . The periodic boundary condition is used. The computational domain with 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 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.
(a)
(b)
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 bywhere , , and . Irregular grids were used in the simulations. The initial values are synoptically shown in Figure 9 for a area.
The finial time of the simulation is . In our test, the cosine bell passes through areas of dense grid and areas of sparse grid. A mesh of is used in the simulation in Figure 10. Figure 10(a) shows the results at time , whereas those in Figure 10(b) were obtained at time . These results enabled us to conclude that our schemes are able to solve this problem without obtaining negative values.
(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]:where .
Appendix
Least-Norm Solution
The linear equationsatisfy the following conditions:(i)The rank of matrix is the same as the rank of extend matrix (ii)The rank of matrix is less than
Thus, there are infinite solutions to linear equations (A.1). The least-norm solution is the minimum value ofwhich also satisfies equation (A.1) Without loss of generality, we assume the rank of is . Then, the matrix can be decomposed as
Then, can be written as
Thus, (A.1) becomeswhich leads
In addition, we have
Then, the derivative of is obtained to determine the extreme values:
In the extreme value point, the least-norm solution is found:
Data Availability
The data in this paper are all open, and this paper is mainly engaged in theoretical analysis; numerical simulation is an auxiliary means of verification.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Basic Research and Development Program of China (Grant no. 2017YFC1502201), Basic Scientific Research and Operation Foundation of the Chinese Academy of Meteorological Sciences (CAMS) (Grant nos. 2018Y007 and 2017Z017), and development fund of State Key Laboratory of Severe Weather (Grant nos. 2018LASW-B12). Thanks to Yuanfu Xie for his suggestions on this article.