Elsevier

Journal of Computational Physics

Volume 344, 1 September 2017, Pages 419-436
Journal of Computational Physics

A vertex-centered and positivity-preserving scheme for anisotropic diffusion problems on arbitrary polygonal grids

https://doi.org/10.1016/j.jcp.2017.04.070Get rights and content

Abstract

We suggest a new positivity-preserving finite volume scheme for anisotropic diffusion problems on arbitrary polygonal grids. The scheme has vertex-centered, edge-midpoint and cell-centered unknowns. The vertex-centered unknowns are primary and have finite volume equations associated with them. The edge-midpoint and cell-centered unknowns are treated as auxiliary ones and are interpolated by the primary unknowns, which makes the final scheme a pure vertex-centered one. Unlike most existing positivity-preserving schemes, the construction of the scheme is based on a special nonlinear two-point flux approximation that has a fixed stencil and does not require the convex decomposition of the co-normal. In order to solve efficiently the nonlinear systems resulting from the nonlinear scheme, Picard method and its Anderson acceleration are discussed. Numerical experiments demonstrate the second-order accuracy and well positivity of the solution for heterogeneous and anisotropic problems on severely distorted grids. The high efficiency of the Anderson acceleration is also shown on reduction of the number of nonlinear iterations. Moreover, the proposed scheme does not have the so-called numerical heat-barrier issue suffered by most existing cell-centered and hybrid schemes. However, further improvements have to be made if the solution is very close to the machine precision and the mesh distortion is very severe.

Introduction

We consider the following diffusion problemdiv(Λu)=f in Ω,u=gD on ΓD,Λun=gN on ΓN, where

  • (1)

    Ω is an open bounded connected polygonal domain in R2;

  • (2)

    f is the source term, belonging to L2(Ω);

  • (3)

    Λ is a symmetric tensor such that (a) Λ is piecewise Lipschitz-continuous on Ω and (b) the set of eigenvalues of Λ is included in [λmin,λmax] with λmin>0;

  • (4)

    Ω=Γ¯DΓ¯N is the boundary of Ω, n denotes the exterior unit normal vector and gD (resp. gN) is the Dirichlet (resp. Neumann) boundary data defined on ΓD (resp. ΓN).

In addition, we assume that f0, gD0 and gN0, so that nonnegativity of the solution u is always assured.

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 first feature is the same as the one in [6] and the second and fourth are different, while the third one was not examined in [6]. The fourth feature makes it most attractive to extend the present scheme to arbitrary 3D grids. The main differences between the new scheme and the one in [13] are the first three features and moreover, the special fixed stencil here enable us to analyze the truncation error rigorously.

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 uh must be understood as their counterparts by replacing uh with u. To guarantee a second-order accuracy, a necessary condition is that the truncation term BK,σϵ,ϵ(uh) should satisfy the condition below|BK,σϵ,ϵ(u)|<Ch2,σEK,KM, where M denotes the set of primary cells, h=maxKMhK and hK 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 belowM(U)U=F(U), where U denotes the unknown vector and M(U) the coefficient matrix. The right-hand side vector F(U) 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:Eu=(xνΩ|Kν|(uνu(xν))2/xνΩ|Kν|u2(xν))1/2,Eq=(KMσEKSσ(F˜K,σ,νF˜K,σ,νext)2/KMσEKSσ(F˜K,σ,νext)2)1/2, where Sσ 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)

  • Z. Sheng et al.

    An improved monotone finite volume scheme for diffusion equation on polygonal meshes

    J. Comput. Phys.

    (2012)
  • C.D. Sijoy et al.

    TRHD: three-temperature radiation-hydrodynamics code with an implicit non-equilibrium radiation transport using a cell-centered monotonic finite volume scheme on unstructured-grids

    Comput. Phys. Commun.

    (2015)
  • J. Wu et al.

    Interpolation-based second-order monotone finite volume schemes for anisotropic diffusion equations on general grids

    J. Comput. Phys.

    (2014)
  • G. Yuan et al.

    Monotone finite volume schemes for diffusion equations on polygonal meshes

    J. Comput. Phys.

    (2008)
  • D.G. Anderson

    Iterative procedures for nonlinear integral equations

    J. Assoc. Comput. Mach.

    (1965)
  • X. Blanc et al.

    A positive scheme for diffusion problems on deformed meshes

    Z. Angew. Math. Mech.

    (2015)
  • F. Brezzi et al.

    A family of mimetic finite difference methods on polygonal and polyhedral meshes

    Math. Models Methods Appl. Sci.

    (2005)
  • Cited by (0)

    This work was supported by the National Natural Science Foundation of China of China (Nos. 91330205, 11671313, 11571226).

    View full text