High order maximum principle preserving discontinuous Galerkin method for convection-diffusion equations
Tao Xiong 111Department of Mathematics, University of Houston, Houston, 77204. E-mail: email@example.com Jing-Mei Qiu 222Department of Mathematics, University of Houston, Houston, 77204. E-mail: firstname.lastname@example.org. The first and second authors are supported by Air Force Office of Scientific Computing YIP grant FA9550-12-0318, NSF grant DMS-1217008, University of Houston. Zhengfu Xu 333Department of Mathematical Science, Michigan Technological University, Houghton, 49931. E-mail: email@example.com. Supported by NSF grant DMS-1316662.
In this paper, we propose to apply the parametrized maximum-principle-preserving (MPP) flux limiter in [Xiong et. al., JCP, 2013] to the discontinuous Galerkin (DG) method for solving the convection-diffusion equations. The feasibility of applying the MPP flux limiters to the DG solution of convection-diffusion 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 Navier-Stokes 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, convection-diffusion equations
In this paper, we propose a parametrized maximum-principle-preserving (MPP) flux limiter for the high order discontinuous Galerkin (DG) finite element method in order to solve the nonlinear convection-diffusion equation
The exact solution of (1.1) satisfies the maximum principle, that is, if
When describes the density of a particular species, , the problem is generally addressed as positivity preserving.
High order shock-capturing numerical methods for convection-dominated problems include the high resolution finite volume (FV) and finite difference (FD) essentially non-oscillatory (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  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 well-known for its flexibility, h-p adaptivity, compactness and high parallel efficiency . Later the DG method was generalized to the convection-dominated diffusion equations. Different types of DG approaches for solving the convection diffusion equations include the local DG (LDG) method , the DG formulation of Cheng and Shu . When a convection-dominated 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
is no longer satisfied.
Two kinds of high order discrete maximum-principle-preserving 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 convection-diffusion equations based on a twice-integrated FV formulation of (1.1) within the FV high order WENO framework . The same technique under the DG framework for hyperbolic conservation laws has been applied to the convection-diffusion equations on a triangular mesh in , however the approach does not work for the Runge-Kutta 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.  by applying the limiter only at the final stage of a high order RK method. The MPP flux limiter has been generalized to convection-diffusion equations under the FD WENO framework in  and the FV WENO framework in . Early discussion of the discrete maximum principle for the convection diffusion equations includes the linear finite element solutions for parabolic equations  with recent developments in [10, 9, 11, 24] and the Petrov-Galerkin finite element method for convection-dominated problems . However, they are under a different framework.
In this paper, we propose to apply the parametrized MPP flux limiter in  to the RKDG method, for solving convection-diffusion 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 , 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  is proposed for the DG schemes with and solution spaces only. Our proposed approach can be viewed as a low-cost easy-to-implement post-processing 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  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 Navier-Stokes 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  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  for directly solving convection-diffusion equations. Then we will apply the parametrized MPP flux limiters developed in  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,
with the cell, cell center to be
The mesh size and let . In the following, for simplicity, we assume uniform mesh sizes .
The DG method in  is defined as follows: find , such that
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
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  for the semi-discrete scheme (2.3) is given as
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
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
where is a first order monotone flux with which the scheme is maximum principle preserving, e.g., the global Lax-Friedrichs flux . The parameter is defined to ensure , for which sufficient inequalities to have are
Let , the parameter can be obtained as follows, for details see :
The local parameter is determined to be
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 multi-stage RK and multi-step 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  in the final time step only to ensure the numerical solution (piecewise polynomials) to be within bounds.
The proposed flux limiter is different from the polynomial rescaling techniques introduced in [31, 33] in several aspects. First of all, in , the entire polynomial (or at least the Gaussian-Lobatto 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 multi-stage RK method may be affected, e.g. see discussions in . Secondly, in  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 .
(accuracy) For the DG method (2.3) with the numerical fluxes (2.4), it has been proved in  that the method is stable with sub-optimal 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 convection-dominated 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].
where and . The LDG method is defined to be: find , such that for any test functions , we have
is the monotone flux for the convection part. For the diffusion part, the numerical fluxes are
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.
For convection-diffusion equations with source terms , the technique in  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 convection-diffusion equation
on a bounded domain of , where is a symmetric semi-positive-definite 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  for two dimensions with rectangular mesh is defined as: find , such that for any test function ,
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 , e.g., the global Lax-Friedrichs flux. Other numerical fluxes are defined by 
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 .
Taking in (2.18), for the cell average, simply we have
where and are with and respectively.
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,
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
when a monotone numerical flux is used under a suitable CFL constraint, which will be specified in the numerical part. Denote
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
For the minimum value part, let
A similar procedure would be applied:
Identify negative values out of the four locally defined numbers ,