High order maximum principle preserving discontinuous Galerkin method for convectiondiffusion equations
Tao Xiong ^{1}^{1}1Department of Mathematics, University of Houston, Houston, 77204. Email: txiong@math.uh.edu JingMei Qiu ^{2}^{2}2Department of Mathematics, University of Houston, Houston, 77204. Email: jingqiu@math.uh.edu. The first and second authors are supported by Air Force Office of Scientific Computing YIP grant FA9550120318, NSF grant DMS1217008, University of Houston. Zhengfu Xu ^{3}^{3}3Department of Mathematical Science, Michigan Technological University, Houghton, 49931. Email: zhengfux@mtu.edu. Supported by NSF grant DMS1316662.
Abstract
In this paper, we propose to apply the parametrized maximumprinciplepreserving (MPP) flux limiter in [Xiong et. al., JCP, 2013] to the discontinuous Galerkin (DG) method for solving the convectiondiffusion equations. The feasibility of applying the MPP flux limiters to the DG solution of convectiondiffusion problem is based on the fact that the cell averages for the DG solutions are updated in a conservative fashion (by using flux difference) even in the presence of diffusion terms. The main purpose of this paper is to address the difficulty of obtaining higher than second order accuracy while maintaining a discrete maximum principle for the DG method solving convection diffusion equations. We found that the proposed MPP flux limiter can be applied to arbitrarily high order DG method. Numerical evidence is presented to show that the proposed MPP flux limiter method does not adversely affect the desired high order accuracy, nor does it require restrictive time steps. Numerical experiments including incompressible NavierStokes equations demonstrate the high order accuracy preserving, the MPP performance, and the robustness of the proposed method.
Keywords: Discoutinuous Galerkin method, high order, maximum principle preserving, flux limiter, convectiondiffusion equations
1 Introduction
In this paper, we propose a parametrized maximumprinciplepreserving (MPP) flux limiter for the high order discontinuous Galerkin (DG) finite element method in order to solve the nonlinear convectiondiffusion equation
(1.1) 
The exact solution of (1.1) satisfies the maximum principle, that is, if
(1.2) 
we have
(1.3) 
When describes the density of a particular species, , the problem is generally addressed as positivity preserving.
High order shockcapturing numerical methods for convectiondominated problems include the high resolution finite volume (FV) and finite difference (FD) essentially nonoscillatory (ENO) and weighted ENO (WENO) methods for the convection part, which are capable of producing solutions with fidelity without spurious oscillations. In this framework, high order central difference is generally used to approximate the second order derivative terms. See the lecture notes [20, 21] and the review paper [22] of Shu and reference therein for more discussion of numerical methods in this aspect. The discontinuous Galerkin (DG) method in the finite element framework is another type method; it was developed by Cockburn et. al. in a series of papers [5, 4, 2, 7] for hyperbolic conservation laws and systems. The DG method has been wellknown for its flexibility, hp adaptivity, compactness and high parallel efficiency [3]. Later the DG method was generalized to the convectiondominated diffusion equations. Different types of DG approaches for solving the convection diffusion equations include the local DG (LDG) method [7], the DG formulation of Cheng and Shu [1]. When a convectiondominated diffusion problem is solved within any of the two previously mentioned frameworks, numerical solutions may exhibit overshoots or undershoots, i.e, a discrete version of the maximum principle
(1.4) 
is no longer satisfied.
Two kinds of high order discrete maximumprinciplepreserving limiters are newly developed for convection dominated problems. One is the polynomial rescaling MPP limiter proposed by Zhang and Shu in [31, 32] for hyperbolic conservation laws. It has been extended to the convectiondiffusion equations based on a twiceintegrated FV formulation of (1.1) within the FV high order WENO framework [30]. The same technique under the DG framework for hyperbolic conservation laws has been applied to the convectiondiffusion equations on a triangular mesh in [33], however the approach does not work for the RungeKutta DG (RKDG) method with order higher than . The high order parametrized MPP flux limiter is developed in [28, 17] for hyperbolic conservation laws, which was later improved by Xiong et. al. [26] by applying the limiter only at the final stage of a high order RK method. The MPP flux limiter has been generalized to convectiondiffusion equations under the FD WENO framework in [14] and the FV WENO framework in [29]. Early discussion of the discrete maximum principle for the convection diffusion equations includes the linear finite element solutions for parabolic equations [12] with recent developments in [10, 9, 11, 24] and the PetrovGalerkin finite element method for convectiondominated problems [18]. However, they are under a different framework.
In this paper, we propose to apply the parametrized MPP flux limiter in [26] to the RKDG method, for solving convectiondiffusion equations. For the convection part, the parametrized MPP flux limiter is proposed to preserve the MPP property for the cell averages and at the final RK stage only. This is different from the polynomial rescaling limiter proposed by Zhang and Shu that preserves the MPP property for the entire polynomial (or at Gaussian quadrature points) per element and at each of the RK stages. For the diffusion part, the parametrized MPP flux limiter is proposed for the DG formulations in [1, 7] with general piecewise () polynomial solution spaces. By the design, the parametrized MPP flux limiter preserves the MPP property of cell averages, thus avoids the main difficulty in the approach of polynomial rescaling limiter. Specifically, in [33], great effort is made to rewrite the updated solution as a convex combination of the solution point values in the current time step via approximating the second derivative term by point values; because of such complications, the limiter in [33] is proposed for the DG schemes with and solution spaces only. Our proposed approach can be viewed as a lowcost easytoimplement postprocessing procedure that modifies the high order numerical fluxes towards a first order one only at the final RK stage for the evolution of cell averages (not for higher moments), in order to preserve the solution cell averages’ MPP property. We remark that the proposed DG solutions (piecewise polynomials) with the parametrized flux limiters might be out of the bound of with the cell averages well bounded by . One can apply the polynomial rescaling limiters as in [31] in the final time step only to ensure the numerical solution (piecewise polynomials) to be within bounds.
The proposed MPP flux limiter in the DG framework is mass conservative due to the flux difference form. It is very efficient due to the fact that the limiter is only applied at the final RK stage and for the cell averages per time step. The parametrized MPP flux limiter for the DG method can be proved to maintain up to 3rd order accuracy under the time step constraint of the original DG method, by following a similar analysis as in [29, 14] for the finite volume and finite difference WENO method. The proof for higher than third order case is very technical and algebraically complicated, thus we rely on extensive numerical tests to showcase that up to fourth order accuracy can be preserved. Extensive numerical tests, including the incompressible NavierStokes system, are presented to demonstrate the robust performance of the proposed approach in preserving the high order accuracy as well as the MPP property.
The outline of the paper is as follows. In Section 2, the DG formulation of Cheng and Shu [1] for one and two dimensions are described. The application of the parametrized MPP flux limiters on the cell averages of the DG solution is presented. Extension to the LDG method will be discussed. Numerical results are provided in Section 3. Finally conclusions are made in Section 4.
2 The MPP flux limiter for the RKDG method
In this section, we will first briefly describe the DG method developed by Cheng and Shu [1] for directly solving convectiondiffusion equations. Then we will apply the parametrized MPP flux limiters developed in [26] to the RKDG method. Both one and two dimensional cases will be presented.
2.1 One dimensional case
Without loss of generality, we assume periodic boundary condition or zero boundary condition with compact support for 1D cases. The spatial domain is discretized by cells,
(2.1) 
with the cell, cell center to be
(2.2) 
The mesh size and let . In the following, for simplicity, we assume uniform mesh sizes .
The DG method in [1] is defined as follows: find , such that
(2.3)  
for any test function and , where and is a piecewise polynomial space with degree up to on the cell . is a monotone flux for the convection term, and are numerical fluxes chosen to be
(2.4) 
here is a positive constant chosen for stability. and are the left and right limit from cell and respectively. is the jump of at the cell interface . A third order strong stability preserving (SSP) RK time discretization [23] for the semidiscrete scheme (2.3) is given as
(2.5)  
The parametrized MPP flux limiters are applied to keep only the cell averages of within the range of . If we take in (2.3) and divided by on both sides, we have
(2.6) 
where the flux . With the third order SPP RK method (2.5), the update of cell averages in equation (2.6) can be written as
(2.7) 
where
(2.8) 
with being the numerical flux obtained by at each RK stage for , respectively.
The MPP flux limiter is proposed to replace the numerical flux by a modified one
(2.9) 
where is a first order monotone flux with which the scheme is maximum principle preserving, e.g., the global LaxFriedrichs flux [8]. The parameter is defined to ensure , for which sufficient inequalities to have are
(2.10)  
(2.11) 
with and
Let , the parameter can be obtained as follows, for details see [26]:
The local parameter is determined to be
(2.12) 
by the consideration to ensure both the upper bound (2.10) and lower bound (2.11) of the cell averages in both cell and .
Remark 2.1.
In the proposed approach, the convection and diffusion terms are treated together for the MPP property of the cell averages. The parametrized flux limiters are applied only for cell averages (not for higher moments) and at the final RK stage as in (2.7). Hence, the proposed flux limiting procedure has low computational cost and is easy to implement. The approach can also be generalized to other multistage RK and multistep methods. We remark that the DG solutions (piecewise polynomials) with the parametrized flux limiters might be out of the bound of with the cell averages well bounded by . One can apply the polynomial rescaling limiters as in [31] in the final time step only to ensure the numerical solution (piecewise polynomials) to be within bounds.
Remark 2.2.
The proposed flux limiter is different from the polynomial rescaling techniques introduced in [31, 33] in several aspects. First of all, in [31], the entire polynomial (or at least the GaussianLobatto quadrature points) over each cell at each of the RK stage are rescaled to satisfy the MPP property. As a result, the temporal accuracy for a multistage RK method may be affected, e.g. see discussions in [31]. Secondly, in [33] the convection and diffusion terms are treated separately; for both terms, great effort has been made to rewrite the updated cell average as a convex combination of point values in the current time step. Such approach introduces extra CFL time step constraint on the DG method; moreover, it is difficult to generalize such approach for the diffusion term with higher than second order accuracy, see Remark 2.2 in [33].
Remark 2.3.
(accuracy) For the DG method (2.3) with the numerical fluxes (2.4), it has been proved in [1] that the method is stable with suboptimal error estimate (th order for polynomial space) for the norm. Numerically, both and norms are of the optimal th order. Regarding the preservation of high order accuracy of the original RK DG method with the proposed MPP flux limiter, similar conclusions can be established following the same line of proof in [29, 14], i.e. without any additional time step condition, (1) for the special case of linear advection problem, the high order accuracy of the original RK DG solutions is maintained with the proposed flux limiter; (2) for the general convectiondominated diffusion problem, up to third order accuracy will be maintained with the flux limiter. In fact, numerically, it can be shown that arbitrary high order accuracy is preserved, under the time step constraint from the linear stability analysis for the DG method [8, 25].
Remark 2.4.
The parametrized MPP flux limiter can also be applied to the local DG (LDG) method for the convectiondiffusion equations [6]. To obtain an LDG formulation for (1.1), first we rewrite it as
(2.13) 
where and . The LDG method is defined to be: find , such that for any test functions , we have
(2.14)  
is the monotone flux for the convection part. For the diffusion part, the numerical fluxes are
(2.16) 
If we take in (2.14), we have the same equation (2.6) for the cell average of , the only difference is the flux given by . The rest of applying the flux limiter would be the same as described above. Similar arguments hold for the two dimensional case in the following subsection.
Remark 2.5.
For convectiondiffusion equations with source terms , the technique in [27] can be used for the source term to ensure the MPP property.
2.2 Two dimensional case
In this subsection, we consider the generalization of the parametrized flux limiter to the two dimensional convectiondiffusion equation
(2.17) 
on a bounded domain of , where is a symmetric semipositivedefinite matrix. Similar observation of (1.3) also holds for the two dimensional case.
For simplicity, in the following, we assume periodic boundary conditions or zero boundary conditions with compact support in each direction. A spatial discretization with rectangular meshes is defined as
where the cell, cell centers and cell sizes are defined by
and , , . For simplicity, in the following, we assume and .
The DG scheme in [1] for two dimensions with rectangular mesh is defined as: find , such that for any test function ,
(2.18)  
Here and is the two dimensional polynomial space with degree up to on the cell , is the outward unit normal vector on the edges. is a monotone numerical flux for the convection part [2], e.g., the global LaxFriedrichs flux. Other numerical fluxes are defined by [33]
(2.19) 
here is the maximum absolute eigenvalue of the symmetric matrix , is a parameter large enough to ensure the stability of the scheme, which will be specified later. are the left and right limit values from the cells adjacent to the edges respectively. On the left boundary of , we have and , while on the right boundary, and . Similarly for and on the top and bottom boundaries of .
With the third order RK time discretization (2.5), the last stage of (2.20) can be written as
(2.21) 
where and . is the integral of the numerical flux along the cell interface , which could be approximated by a numerical quadrature. At each fixed quadrature point , is defined the same as (2.8). Similarly for .
Let and , numerically to preserve the cell averages within the range , we are looking for the type of limiters,
(2.22)  
(2.23) 
such that
(2.24) 
where and are first order monotone fluxes which can form a maximum principle preserving first order scheme similarly as the one dimensional case. (2.22)(2.24) form coupled inequalities for the limiting parameters . In each cell , as the 1D case, the MPP flux limiters can be parametrized in the sense that we can find a group of numbers , such that the numerical solutions of (2.21) satisfy the MPP property (2.24) with
For the maximum value case, let
(2.25) 
when a monotone numerical flux is used under a suitable CFL constraint, which will be specified in the numerical part. Denote
(2.26) 
The coupled inequalities (2.22)(2.24) can be rewritten as
(2.27) 
To decouple the inequality (2.27), for the specific cell , two steps are followed:

Identify positive values out of the four locally defined numbers , , , ;

Corresponding to those positive values, collectively, the limiting parameters can be defined. For example, if and , then
(2.28)
For the minimum value part, let
(2.29) 
The coupled inequalities (2.22)(2.24) can be rewritten as
(2.30) 
A similar procedure would be applied:

Identify negative values out of the four locally defined numbers ,