1 Introduction

High order maximum principle preserving discontinuous Galerkin method for convection-diffusion equations

Tao Xiong 111Department of Mathematics, University of Houston, Houston, 77204. E-mail: txiong@math.uh.edu Jing-Mei Qiu 222Department of Mathematics, University of Houston, Houston, 77204. E-mail: jingqiu@math.uh.edu. 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: zhengfux@mtu.edu. Supported by NSF grant DMS-1316662.

Abstract

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

## 1 Introduction

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

 ut+f(u)x=a(u)xx,u(x,0)=u0(x). (1.1)

The exact solution of (1.1) satisfies the maximum principle, that is, if

 uM=maxxu0(x),um=minxu0(x), (1.2)

we have

 u(x,t)∈[um,uM],∀t>0. (1.3)

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 [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 well-known for its flexibility, h-p adaptivity, compactness and high parallel efficiency [3]. 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 [7], the DG formulation of Cheng and Shu [1]. 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

 um≤unj≤uM,∀n,j, (1.4)

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 [30]. The same technique under the DG framework for hyperbolic conservation laws has been applied to the convection-diffusion equations on a triangular mesh in [33], 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. [26] 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 [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 Petrov-Galerkin finite element method for convection-dominated 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 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 [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 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 [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 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 [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 convection-diffusion 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,

 a=x12

with the cell, cell center to be

 Ij=(xj−12,xj+12),xj=12(xj−12+xj+12),∀j=1,⋯N. (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

 ∫Ij(uh)tvhdx − ∫Ijf(uh)(vh)xdx−∫Ija(uh)(vh)xxdx (2.3) + (^f(u−h,u+h)v−h)j+12−(ˆf(u−h,u+h)v+h)j−12−(˜a(uh)xv−h)j+12 + (˜a(uh)xv+h)j−12+(ˆa(uh)(vh)−x)j+12−(ˆa(uh)(vh)+x)j−12 = 0,

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

 ˜a(uh)x=[a(uh)][uh]((uh)−x+αh[a(uh)]),ˆa(uh)=a(u+h), (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 semi-discrete scheme (2.3) is given as

 u(1) = un+ΔtL(un), u(2) = un+Δt(14L(un)+14L(u(1))), (2.5) un+1 = un+Δt(16L(un)+23L(u(2))+16L(u(1))).

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

 ddt¯uh+1h(^Hj+12−^Hj−12)=0, (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

 (¯uh)n+1j=(¯uh)nj−λ(^Hrkj+12−^Hrkj−12), (2.7)

where

 ^Hrkj+12≐16^Hnj+12+23^H(2)j+12+16^H(1)j+12, (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

 ~Hrkj+12=θj+12(^Hrkj+12−^hj+12)+^hj+12, (2.9)

where is a first order monotone flux with which the scheme is maximum principle preserving, e.g., the global Lax-Friedrichs flux [8]. The parameter is defined to ensure , for which sufficient inequalities to have are

 λθj−12(^Hrkj−12−^hj−12)−λθj+12(^Hrkj+12−^hj+12)−ΓMj ≤ 0, (2.10) λθj−12(^Hrkj−12−^hj−12)−λθj+12(^Hrkj+12−^hj+12)−Γmj ≥ 0, (2.11)

with and

 ΓMj=uM−(¯uh)nj+λ(^hj+12−^hj−12)≥0,Γmj=um−(¯uh)nj+λ(^hj+12−^hj−12)≤0.

Let , the parameter can be obtained as follows, for details see [26]:

1. Assume where and are designed to preserve the upper bound by equation (2.10),

1. If and ,

2. If and ,

3. If and ,

4. If and ,

• If equation (2.10) is satisfied with , then

• If equation (2.10) is not satisfied with , then

2. Similarly assume

 θj−12∈[0,Λm−12,Ij],θj+12∈[0,Λm+12,Ij],

where and are designed to preserve the lower bound by equation (2.11),

1. If and ,

2. If and ,

3. If and ,

4. If and ,

• If equation (2.11) is satisfied with , then

• If equation (2.11) is not satisfied with , then

The local parameter is determined to be

 θj+12=min(ΛM+12,Ij,ΛM−12,Ij+1,Λm+12,Ij,Λm−12,Ij+1), (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 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 [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 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 [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 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].

###### Remark 2.4.

The parametrized MPP flux limiter can also be applied to the local DG (LDG) method for the convection-diffusion equations [6]. To obtain an LDG formulation for (1.1), first we rewrite it as

 ut+f(u)x=(γ(u)q)x,q−Γ(u)x=0, (2.13)

where and . The LDG method is defined to be: find , such that for any test functions , we have

 ∫Ij(uh)tvhdx−∫Ij(f(uh)−γ(uh)qh)(vh)xdx+(^f−^γ^q)j+12(vh)−j+12−(^f−^γ^q)j−12(vh)+j−12=0, (2.14) ∫Ijqhwhdx+∫IjΓ(uh)(wh)xdx−^Γj+12(wh)−j+12+^Γj−12(wh)+j−12=0.

is the monotone flux for the convection part. For the diffusion part, the numerical fluxes are

 ^γ=Γ(u+h)−Γ(u−h)u+h−u−h,^q=q−h,^Γ=Γ(u+h). (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 convection-diffusion 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 convection-diffusion equation

 ut+∇⋅F(u)=∇⋅(A∇u),F(u)=(f(u),g(u)), (2.17)

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

 a=x12

where the cell, cell centers and cell sizes are defined by

 Kij=Ii×Jj,Ii=(xi+12,xi−12),Jj=(yj−12,yj+12),
 hxi=xi+12−xi−12,xi=12(xi+12+xi−12),hyj=yj+12−yj−12,yj=12(yj−12+yj+12),

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 ,

 ∫Kij (uh)t vhdx−∫KijF(uh)⋅∇vhdx−∫Kijuh∇⋅(A∇vh)dx (2.18) + ∫∂Kij(˜uhA∇vh)⋅nds+∫∂Kij(^F⋅n−ˆA∇uh⋅n)vhds=0.

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 Lax-Friedrichs flux. Other numerical fluxes are defined by [33]

 ˆA∇uh⋅n=A(u−h)∇u−h⋅n+αΛh(uouth−uinh),˜uhA=u+hA(u+h), (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 .

Taking in (2.18), for the cell average, simply we have

 ddt¯uh + 1hx(1hy∫Jj^H(xi+12,y)dy−1hy∫Jj^H(xi−12,y)dy) (2.20) + 1hy(1hx∫Ii^G(x,yj+12)dx−1hx∫Ii^G(x,yj−12)dx)=0.

where and are with and respectively.

With the third order RK time discretization (2.5), the last stage of (2.20) can be written as

 ¯un+1h=¯unh−λx(^Hrki+12,j−^Hrki−12,j)−λy(^Grki,j+12−^Grki,j−12), (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,

 ~Hi+12,j = θi+12,j(^Hrki+12,j−^hi+12,j)+^hi+12,j, (2.22) ~Gi,j+12 = θi,j+12(^Grki,j+12−^gi,j+12)+^gi,j+12, (2.23)

such that

 um≤(¯uh)ni,j−λx(~Hi+12,j−~Hi−12,j)−λy(~Gi,j+12−~Gi,j−12)≤uM, (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

 (θi−12,j,θi+12,j,θi,j−12,θi,j+12)∈[0,ΛL,i,j]×[0,ΛR,i,j]×[0,ΛD,i,j]×[0,ΛU,i,j].

For the maximum value case, let

 ΓMi,j=uM−((¯uh)ni,j−λx(^hi+12,j−^hi−12,j)−λy(^gi,j+12−^gi,j−12))≥0, (2.25)

when a monotone numerical flux is used under a suitable CFL constraint, which will be specified in the numerical part. Denote

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Fi−12,j=λx(^Hrki−12,j−^hi−12,j),Fi+12,j=−λx(^Hrki+12,j−^hi+12,j),Fi,j−12=λy(^Grki,j−12−^gi,j−12),Fi,j+12=−λy(^Grki,j+12−^gi,j+12). (2.26)

The coupled inequalities (2.22)-(2.24) can be rewritten as

 θi+12,jFi+12,j+θi−12,jFi−12,j+θi,j+12Fi,j+12+θi,j−12Fi,j−12≤ΓMi,j, (2.27)

To decouple the inequality (2.27), for the specific cell , two steps are followed:

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

2. Corresponding to those positive values, collectively, the limiting parameters can be defined. For example, if and , then

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩ΛMi+12,j,ΛMi−12,j=min(ΓMi,jFi+12,j+Fi−12,j,1),ΛMi,j−12,ΛMi,j+12=1. (2.28)

For the minimum value part, let

 Γmi,j=um−((¯uh)ni,j−λx(^hi+12,j−^hi−12,j)−λy(^gi,j+12−^gi,j−12))≤0. (2.29)

The coupled inequalities (2.22)-(2.24) can be rewritten as

 Γmi,j≤θi+12,jFi+12,j+θi−12,jFi−12,j+θi,j+12Fi,j+12+θi,j−12Fi,j−12. (2.30)

A similar procedure would be applied:

1. Identify negative values out of the four locally defined numbers ,