A vertex-centered and positivity-preserving scheme for anisotropic diffusion problems on arbitrary polygonal grids☆
Introduction
We consider the following diffusion problem where
- (1)
Ω is an open bounded connected polygonal domain in ;
- (2)
f is the source term, belonging to ;
- (3)
Λ is a symmetric tensor such that (a) Λ is piecewise Lipschitz-continuous on Ω and (b) the set of eigenvalues of Λ is included in with ;
- (4)
is the boundary of Ω, n denotes the exterior unit normal vector and (resp. ) is the Dirichlet (resp. Neumann) boundary data defined on (resp. ).
Diffusion equations of the form (1) are at the core of complex models such as those in radiation hydrodynamics (RHD), magnetohydrodynamics (MHD), reservoir modeling, and so on. Finite volume (FV) discretizations in these applications are well received due to some important characteristics such as simplicity and local conservation. The readers are referred to, e.g., [6], [8], [9], for some recent developments. In many situations, one of the significant requirements of the FV schemes is that the discrete solution should be nonnegative. A scheme with such a property is usually referred to as positivity-preserving or monotone. In some applications such as RHD and reservoir simulations, the meshes are usually distorted and the media are heterogeneous and anisotropic, which may cause a scheme to violate the positivity-preserving property and exhibit nonphysical oscillations. Many scientists have made great efforts in the development of FV schemes that preserve the positivity-preserving property and at the same time, have approximately a second-order accuracy on severely distorted meshes in the case that the diffusion tensor is taken to be highly anisotropic, heterogeneous, and/or discontinuous.
As has been shown in [5], [15], [21] that, no linear consistent and conservative nine-point scheme can unconditionally respect the positivity-preserving property. One must therefore look for nonlinear schemes. In [22], C. Le Potier suggested a nonlinear two-point flux approximation to construct a cell-centered scheme on triangular meshes. Since then, this approach has been developed to obtain positivity-preserving schemes [3], [14], [16], [26], [35] or extreme-preserving ones [12], [25] on general grids with general diffusion tensors. Two key ingredients of these works are the convex decomposition of the co-normal and the positivity-preserving interpolation of the auxiliary unknowns. By comparison, the later is quite difficult, even more difficult than the design of positivity-preserving FV scheme itself. The methods in [3], [16], [26], [35] generally introduce auxiliary unknowns at mesh vertices or edge midpoints, whose values are interpolated from the neighboring cell-centered unknowns. When some of the weights in the interpolation formula are negative, the authors suggested that the corresponding auxiliary unknowns must be computed by some lower-order but positivity-preserving interpolation methods, such as the inverse distance weighting method [16]. Unfortunately, this technique usually leads to a great loss of accuracy on largely distorted meshes even in the case of a constant diffusion coefficient with a smooth solution. An interpolation-free approach based on local repositioning of cell centers was proposed in [17], but it can be used only for meshes with only one discontinuous interface per cell. A nonlinear two-point flux approximation a little different from that of Potier's was first proposed in [32] and further developed in [13] where both the positive interpolation of the auxiliary unknowns and the convex decomposition of the co-normal are not required. Therefore, it is possible for us to use arbitrary second-order interpolation algorithms and a relatively arbitrary co-normal decomposition. However, the truncation error in this approach is a little difficult to be analyzed rigorously.
To our knowledge, all the nonlinear positivity-preserving FV schemes are cell-centered except those in [23] and [6]. In [23], a nonlinear correction technique is applied to a class of hybrid schemes to insure the discrete maximum principle. Here we are more concerned about the nonlinear scheme in [6] where both the cell-centered and vertex-centered unknowns are treated as primary ones. This scheme is constructed on the primary mesh and its dual counterpart. Like most existing positivity-preserving schemes, it requires the convex decomposition of the co-normal. Moreover, it has FV equations for the two sets of primary unknowns so that no interpolation is involved. However, this nonlinear scheme can not reach the second-order accuracy in the case of discontinuous diffusion coefficients and the authors did not know how to improve it. We show in this work that the convex decomposition of the co-normal may result in the violation of the linearity-preserving property. To construct a vertex-centered nonlinear FV scheme with a second-order accuracy is among the main motivations of this paper.
Recently, the authors in [20] pointed out that many existing cell-centered and hybrid FV schemes for nonlinear parabolic equations, including the mimetic finite difference schemes [4], [19], suffer the so-called heat-barrier issue. More explicitly, any schemes based on the harmonic averaging of cell-centered diffusion coefficients will break down when some of these coefficients go to zero or their ratio grows, which results in totally wrong numerical solution profiles in some strongly nonlinear problems such as the propagation of a nonlinear heat wave in a cold media. To our knowledge, the report of this issue can be traced back to [2]. Our numerical experiments indicate that all the cell-centered and hybrid linearity-preserving schemes we studied before, including those nonlinear ones in [12], [13], [32], suffer the same problem. The vertex-centered schemes suggested in [31], [33] do not have the numerical heat-barrier issue but they are not positivity-preserving. To design a positivity-preserving scheme that can overcome the heat-barrier issue is another main motivation here.
The efficiency of nonlinear solvers is another important issue. At the present, in solving the nonlinear systems resulting from the aforementioned nonlinear FV schemes, Picard method or fixed-point iteration is frequently used because it preserves the M-matrix structure of the coefficient matrix, which is required to guarantee the positivity of the discrete solutions. In some extreme cases such as largely distorted meshes, Picard method can be very slow, which motivates us to seek an efficient nonlinear solver that does not spoil the M-matrix structure of the coefficient matrix. Anderson acceleration (AA) has been provided for nonlinear algebraic or transcendental equations in [1], and has been successfully used in electronic-structure computations (see, e.g. [10], [34] and the references therein). It has been shown to be potentially useful in a broader variety of applications [29], [30]. Recently, a convergence result of AA has been obtained in [28]. More interesting is that AA algorithm employs the same matrix structure of Picard method. Based on AA and the adaptive inexact solution of linear systems, an efficient strategy was proposed in [18] for accelerating the Picard method. We mention that, some authors [24] investigated the non-fixed-point iteration such as Newton method, however, methods of this kind usually spoil the M-matrix structure.
In the present work, we further develop the idea in [13] to construct a vertex-centered and positivity-preserving scheme on general meshes, and discuss also its AA algorithm. Except for some standard features such as the local conservation, the second-order accuracy and the validness for arbitrary star-shaped meshes with arbitrary diffusion tensors, the present scheme has the following special ones:
- 1.
Its construction is based on the primary mesh and its dual counterpart;
- 2.
It is a pure vertex-centered one. The cell-centered and edge-midpoint unknowns are computed by simple interpolation algorithms;
- 3.
It does not involve the harmonic averaging of cell-centered diffusion coefficients so that it does not suffer the numerical heat-barrier issue;
- 4.
The convex decomposition of the co-normal is not required and the nonlinear two-point flux approximation has a fixed stencil.
The rest contents of the paper is organized as follows. In Section 2, the construction of the new scheme is described based on a new nonlinear two-point flux approximation. Then the truncation error is analyzed in Section 3. In Section 4, a simple interpolation algorithm for the auxiliary unknowns is suggested and in Section 5, the nonlinear iterations, including the Picard method and its Anderson acceleration, are discussed. Finally, numerical experiments are presented in Section 6 and some concluding remarks are given in the last section.
Section snippets
Notations and assumptions
Suppose that Ω is partitioned into a number of non-overlapped polygonal cells that form the so-called primary mesh, see the mesh with solid line segments in Fig. 1. The cell center is defined as the geometric center whose coordinates are the simple averages of the cell vertices. Other definitions for cell centers are also allowed, which may result in more complicated interpolation algorithms. Assume that all primary cells are star-shaped with respect to their cell centers. Then, each primary
Analysis of the truncation error
Throughout this section, without special mention, all the quantities defined in the previous section and involved must be understood as their counterparts by replacing with u. To guarantee a second-order accuracy, a necessary condition is that the truncation term should satisfy the condition below where denotes the set of primary cells, and is the diameter of cell K. In the following discussion, notation (u) will be dropped for
The interpolation of the auxiliary unknowns
In design of both vertex-centered and cell-centered FV schemes, the auxiliary unknowns must be interpolated by primary unknowns. A desirable interpolation algorithm is usually required to be second-order, simple, positivity-preserving, topology-independent, discontinuity-independent and so on, which is quite difficult for cell-centered schemes. As for the present scheme, since the auxiliary unknowns are defined at the centers and edge midpoints of the primary cells, they can be easily
Monotonicity and iteration of the nonlinear scheme
The nonlinear scheme (28) can be formulated as the matrix form below where U denotes the unknown vector and the coefficient matrix. The right-hand side vector is generated by the source term and the boundary data.
Numerical experiments
In this section, we examine the numerical performance of the vertex-centered and positive-preserving scheme (VPPS) discussed in the previous sections. The scheme with the Picard (resp. Anderson acceleration) method is denoted as VPPS-P (resp. VPPS-AA) for simplicity. We investigate the discrete relative errors for both the solution and the flux: where is an area associated with σ
Conclusion
In this paper, a new nonlinear positivity-preserving FV scheme is suggested for the numerical solution of anisotropic diffusion problems on general polygonal grids with star-shaped cells. This scheme differs from existing nonlinear FV schemes in the way that it is of a pure vertex-centered type and relies on a new nonlinear two-point flux approximation that has a fixed stencil. The key ingredient is to abandon the traditional convex decomposition of the co-normal vector, which seems to a unique
Acknowledgements
The authors thank the anonymous reviewers for their careful readings and useful suggestions.
References (35)
- et al.
An efficient cell-centered diffusion scheme for quadrilateral grids
J. Comput. Phys.
(2009) - et al.
On the nonexistence of monotone linear schema for some linear parabolic equations
C. R. Acad. Sci. Paris, Ser. I
(2005) - et al.
A small stencil and extremum-preserving scheme for anisotropic diffusion problems on arbitrary 2D and 3D meshes
J. Comput. Phys.
(2013) - et al.
Sufficient criteria are necessary for monotone control volume methods
Appl. Math. Lett.
(2009) - et al.
Monotone finite volume schemes for diffusion equations on unstructured triangular and shape-regular polygonal meshes
J. Comput. Phys.
(2007) - et al.
Interpolation-free monotone finite volume method for diffusion equations on polygonal meshes
J. Comput. Phys.
(2009) - et al.
The mimetic finite difference method for elliptic and parabolic problems with a staggered discretization of diffusion coefficient
J. Comput. Phys.
(2016) Schéma volumes finis monotones pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés
C. R. Acad. Sci. Paris, Ser. I
(2005)- et al.
A nonlinear correction and maximum principle for diffusion operators discretized using hybrid schemes
C. R. Acad. Sci. Paris, Ser. I
(2012) - et al.
The finite volume scheme preserving extremum principle for diffusion equations on polygonal meshes
J. Comput. Phys.
(2011)