1 Introduction

High order maximum principle preserving finite volume method for convection dominated problems

Pei Yang 111Department of Mathematics, University of Houston, Houston, 77204. E-mail: peiyang@math.uh.edu. Tao Xiong 222Department of Mathematics, University of Houston, Houston, 77204. E-mail: txiong@math.uh.edu. Jing-Mei Qiu 333Department of Mathematics, University of Houston, Houston, 77204. E-mail: jingqiu@math.uh.edu. The first, second and the third authors are supported by Air Force Office of Scientific Computing YIP grant FA9550-12-0318, NSF grant DMS-1217008. Zhengfu Xu 444Department of Mathematical Science, Michigan Technological University, Houghton, 49931. E-mail: zhengfux@mtu.edu. Supported by NSF grant DMS-1316662.

Abstract

In this paper, we investigate the application of the maximum principle preserving (MPP) parametrized flux limiters to the high order finite volume scheme with Runge-Kutta time discretization for solving convection dominated problems. Such flux limiter was originally proposed in [Xu, Math. Comp., 2013] and further developed in [Xiong et. al., J. Comp. Phys., 2013] for finite difference WENO schemes with Runge-Kutta time discretization for convection equations. The main idea is to limit the temporal integrated high order numerical flux toward a first order MPP monotone flux. In this paper, we generalize such flux limiter to high order finite volume methods solving convection-dominated problems, which is easy to implement and introduces little computational overhead. More importantly, for the first time in the finite volume setting, we provide a general proof that the proposed flux limiter maintains high order accuracy of the original WENO scheme for linear advection problems without any additional time step restriction. For general nonlinear convection-dominated problems, we prove that the proposed flux limiter introduces up to modification to the high order temporal integrated flux in the original WENO scheme without extra time step constraint. We also numerically investigate the preservation of up to ninth order accuracy of the proposed flux limiter in a general setting. The advantage of the proposed method is demonstrated through various numerical experiments.

Keywords: Convection diffusion equation; High order WENO scheme; Finite volume method; Maximum principle preserving; Flux limiters

## 1 Introduction

Recently, there is a growing interest in designing high order maximum principle preserving (MPP) schemes for solving scalar convection-dominated problems [17, 16, 15, 9, 10, 12], positivity preserving schemes for compressible Euler and Navier-Stokes equations [8, 13, 11, 18]. The motivation of this family of work arises from the observation that many existing high order conservative methods break down when simulating fluid dynamics in extreme cases such as near-vacuum state. To illustrate the purpose of the family of the MPP methods, we shall consider the solution to the following problem

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

with . The solution to (1.1) satisfies the maximum principle, i.e.

 if uM=maxxu0(x),um=minxu0(x), then u(x,t)∈[um,uM]. (1.2)

Within the high order finite volume (FV) Runge-Kutta (RK) weighted essentially non-oscillatory (WENO) framework, we would like to maintain a discrete form of (1.2):

 if uM=maxxu0(x),um=minxu0(x), then ¯unj∈[um,uM] for any n, j, (1.3)

where approximates the cell average of the exact solution with high order accuracy on a given th spatial interval at time .

Efforts for designing MPP high order schemes to solve (1.1) can be found in recent work by Zhang et al. [16, 19], as a continuous research effort to design high order FV and discontinuous Galerkin (DG) MPP schemes based on a polynomial rescaling limiter on the reconstructed (for FV) or representing (for DG) polynomials [17]. This approach requires the updated cell average to be written as a convex combination of some local quantities within the range . For convection-diffusion problems which do not have a finite speed of propagation, it is difficult to generalize such approach to design MPP schemes that are higher than third order accurate. In [9], an alternative approach via a parametrized flux limiter, developed earlier by Xu et al. [15, 12], is proposed for the finite difference (FD) RK WENO method in solving convection diffusion equations. The flux limiter is applied to convection and diffusion fluxes together to achieve (1.3) for the approximated point values in the finite difference framework. In this paper, we continue our effort in applying the MPP flux limiters to high order FV RK WENO methods to maintain (1.3) with efficiency. Furthermore, we provide some theoretical analysis on the preservation of high order accuracy for the proposed flux limiter in FV framework. Finally, we remark that our current focus is on convection-dominated diffusion problems for which explicit temporal integration proves to be efficient. For the regime of medium to large diffusion, where implicit temporal integration is needed for simulation efficiency, we refer to earlier work in [5, 3, 2, 4] and references therein for the construction of the MPP schemes with finite element framework. The generalization of the current flux limiter is not yet available and is subject to future investigation.

The MPP methods in [17, 15, 12] are designed base on the observation that first order monotone schemes in general satisfy MPP property (1.3) with proper Courant-Friedrichs-Lewy (CFL) numbers, while regular high order conservative schemes often fail to maintain (1.3). The MPP flux limiting approach is to seek a linear combination of the first order monotone flux with the high order flux, in the hope of that such combination can achieve both MPP property and high order accuracy under certain conditions, e.g. some mild time step constraint. This line of approach is proven to be successful in [12, 9] for the FD RK WENO schemes and it is later generalized to the high order semi-Lagrangian WENO method for solving the Vlasov-Poisson system [14]. A positivity preserving flux limiting approach is developed in [13] to ensure positivity of the computed density and pressure for compressible Euler simulations. Technically, the generalization of such MPP flux limiters from FD WENO [9] to FV WENO method is rather straightforward. Taking into the consideration that FV method offers a more natural framework for mass conservation and flexibility in handling irregular computational domain, we propose to apply the MPP flux limiters to the high order FV RK WENO method to solve (1.1). The proposed flux limiting procedure is rather easy to implement even with the complexity of the flux forms in multi-dimensional FV computation. Moreover, a general theoretical proof on preserving both MPP and high order accuracy without additional time step constraint can be done for FV methods when solving a linear advection equation; such result does not hold for high order FD schemes [12].

In this paper, for the first time, we establish a general proof that, there is no further time step restriction, besides the CFL condition under the linear stability requirement, to preserve high order accuracy when the high order flux is limited toward an upwind first order flux for solving linear advection problem, when the parametrized flux limiters are applied to FV RK WENO method. In other words, both the MPP property and high order accuracy of the original scheme can be maintained without additional time step constraint. For a general nonlinear convection problem, we prove that the flux limiter preserves up to third order accuracy and the discrete maximum principle with no further CFL restriction. This proof relies on tedious Taylor expansions, and it is difficult to generalize it to results with higher order accuracy (fourth order or higher). On the other hand, such analysis can be extended to a convection-dominated diffusion problem as done in [9]. Furthermore, numerical results indicate that mild CFL restriction is needed for the MPP flux limiting finite volume scheme without sacrificing accuracy. For more discussions, see Section 3.

The paper is organized as follows. In Section 2, we provide the numerical algorithm of the high order FV RK WENO schemes with MPP flux limiters. In Section 3, theoretical analysis is given for a linear advection problem and general nonlinear problems. Numerical experiments are demonstrated in Section 4. We give a brief conclusion in Section 5.

## 2 A MPP FV method

In this section, we propose a high order FV scheme for the convection-diffusion equation. In the proposed scheme, the high order WENO reconstruction of flux is used for the convection term, while a high order compact reconstruction of flux is proposed for the diffusion term.

For simplicity, we first consider a one dimensional (1D) case. The following uniform spatial discretization is used for a 1D bounded domain ,

 a=x12

with the computational cell and cell center defined as

 Ij=[xj−12,xj+12], xj=12(xj−12+xj+12), j=1,2,⋯,N. (2.2)

Let denote approximation to the cell average of over cell . The FV scheme is designed by integrating equation (1.1) over each computational cell and then dividing it by ,

 d¯ujdt=−1Δx(^HCj+12−^HCj−12)+1Δx(^HDj+12−^HDj−12), (2.3)

where and are the numerical fluxes for convection and diffusion terms respectively.

For the convection term, one can adopt any monotone flux. For example, in our simulations, we use the Lax-Friedrichs flux

 ^HCj+12(u−j+12,u+j+12)=12(f(u−j+12)+αu−j+12)+12(f(u+j+12)−αu+j+12), α=maxum≤u≤uM|f′(u)|. (2.4)

Here , where is obtained by reconstructing a order polynomial whose averages agree with those in a left-biased stencil ,

 1Δx∫IlP(x)dx=¯ul, l=j−k,⋯,j+k.

The reconstruction procedure for can be done similarly from a right-biased stencil. To suppress oscillation around discontinuities and maintain high order accuracy around smooth regions of the solution, the WENO mechanism can be incorporated in the reconstruction. Details of such procedure can be found in [1].

For the diffusion term, we propose the following compact reconstruction strategy for approximating fluxes at cell boundaries . Without loss of generality, we consider a fourth order reconstruction, while similar strategies can be extended to schemes with arbitrary high order. Below we let denote approximation to the point values of at .

1. Reconstruct from the cell averages by constructing a cubic polynomial , such that

 1Δx∫IlP(x)dx=¯ul, l=j−1,⋯,j+2.

Then , . We use to denote such reconstruction procedure,

 (uj−1,uj,uj+1,uj+2)=R1(¯uj−1,¯uj,¯uj+1,¯uj+2).

As a reference, the reconstruction formulas for are provided below,

 uj−1=1112¯uj−1+524¯uj−16¯uj+1+124¯uj+2, uj=−124¯uj−1+1312¯uj−124¯uj+1, uj+1=−124¯uj+1312¯uj+1−124¯uj+2, uj+2=124¯uj−1−16¯uj+524¯uj+1+1112¯uj+2.
2. Construct an interplant such that

 Q(xl)=a(ul), l=j−1,⋯,j+2.

Then let

 ^HDj+12=Q′(x)|xj+12.

Such procedure is denoted as

 ^HDj+12=R2(a(uj−1),a(uj),a(uj+1),a(uj+2)).

As a reference, we provide the formula for below

 ^HDj+12=124a(uj−1)−98a(uj)+98a(uj+1)−124a(uj+2).
###### Remark 2.1.

The reconstruction processes for and operators are designed such that is reconstructed from a compact stencil with a given order of accuracy. Because of such design, for the linear diffusion term , and can be combined and the strategy above turns out to be a classical fourth order central difference from a five-cell stencil with

 ^HDj+12=1Δx(12¯uj−1−1512¯uj+1512¯uj+1−112¯uj+2).

If each of () in Step 1 is reconstructed from symmetrical stencils (having the same number of cells from left and from right), the reconstruction of will depend on a much wider stencil . Such non-compact way of reconstructing numerical fluxes for diffusion terms will introduce some numerical instabilities when approximating nonlinear diffusion terms in our numerical tests, whereas the proposed compact strategy does not encounter such difficulty.

We use the following third order total variation diminishing (TVD) RK method [6] for the time discretization of (2.3), which reads

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

where denotes the right hand side of equation (2.3). Here and , denote the numerical solution of at time and corresponding RK stages. The fully discretized scheme (2) can be rewritten as

 ¯un+1j=¯unj−λ(^Hrkj+12−^Hrkj−12) (2.6)

with and

 ^Hrkj+12=16(^HC,nj+12−^HD,nj+12)+16(^HC,(1)j+12−^HD,(1)j+12)+23(^HC,(2)j+12−^HD,(2)j+12).

Here are the numerical fluxes at the intermediate stages in the RK scheme (2).

It has been known that the numerical solutions from schemes with a first order monotone flux for the convection term together with a first order flux for the diffusion term satisfy the maximum principle, if the time step is small enough [19]. However, if the numerical fluxes are of high order such as the one from the reconstruction process proposed above, the MPP property for the numerical solutions does not necessarily hold under the same time step constraint. Next we apply the parametrized flux limiters proposed in [12] to the scheme (2.6) to preserve the discrete maximum principle (1.3).

We modify the numerical flux in equation (2.6) with

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

by carefully seeking local parameters , such that the numerical solutions enjoy the MPP property yet is as close to as possible. In other words, is as close to the original high order flux as possible. Here denotes the first order flux for convection and diffusion terms, using which in the scheme (2.3) with a forward Euler time discretization guarantees the maximum principle of numerical solutions. For example, we can take

 ^hj+12=^hCj+12−^hDj+12=12(f(¯uj)+α¯uj)+12(f(¯uj+1)−α¯uj+1)−a(¯uj+1)−a(¯uj)Δx

with . The goal of the procedures outlined below is to adjust , so that with the modified flux , the numerical solutions satisfy the maximum principle,

 um≤¯unj−λ(~Hrkj+12−~Hrkj−12)≤uM,∀j. (2.8)

Detailed procedures in decoupling the above inequalities have been intensively discussed in our previous work, e.g. [12]. Below we only briefly describe the computational algorithm for the proposed limiter.

Let and

 ΓMj≐uM−(¯unj−λ(^hj+12−^hj−12)), Γmj≐um−(¯unj−λ(^hj+12−^hj−12)).

The MPP property is satisfied with the modified flux (2.7) when the following inequalities are hold,

 λθj−12Fj−12−λθj+12Fj+12−ΓMj≤0, (2.9) λθj−12Fj−12−λθj+12Fj+12−Γmj≥0. (2.10)

We first consider the inequality (2.9). We seek a local pair of numbers such that (1) and is as close to as possible, (2) for any , the inequality (2.9) holds. The inequality (2.9) can be decoupled based on the following four different cases:

1. If and , then

2. If and , then

3. If and , then

4. If and , then

 (ΛM−12,Ij,ΛM+12,Ij)=(min(1,ΓMjλFj−12−λFj+12),min(1,ΓMjλFj−12−λFj+12)).

Similarly, we can find a local pair of numbers such that for any

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

(2.10) holds. There are also four different cases:

1. If and , then

2. If and , then

3. If and , then

4. If and , then

 (Λm−12,Ij,Λm+12,Ij)=(min(1,ΓmjλFj−12−λFj+12),min(1,ΓmjλFj−12−λFj+12)).

Finally, the local limiter parameter at the cell boundary is defined as

 θj+12=min(ΛM+12,Ij,Λm+12,Ij,ΛM−12,Ij+1,Λm−12,Ij+1), (2.11)

so that the numerical solutions , satisfy the maximum principle.

The extension of the FV RK scheme and the MPP flux limiter from 1D case to two dimensional (2D) convection-diffusion problems is straightforward. For example, we consider a 2D problem on a rectangular domain ,

 ut+f(u)x+g(u)y=a(u)xx+b(u)yy. (2.12)

Without loss of generality, we consider a set of uniform mesh

 a=x12
 c=y12

with . A semi-discrete FV discretization of (2.12) gives

 ddt¯ui,j +1Δx(^fi+12,j−^fi−12,j)+1Δy(^gi,j+12−^gi,j−12) =1Δx(ˆ(ax)i+12,j−ˆ(ax)i−12,j)+1Δy(ˆ(by)i,j+12−ˆ(by)i,j−12), (2.13)

where and is the average of the flux over the right boundary of cell . , , can be defined similarly. The flux is evaluated by applying the Gaussian quadrature rule for integration,

 ^fi+12,j=12Σigωigf(ui+12,ig). (2.14)

Here represents the summation over the Gaussian quadratures with being quadrature weights and is the approximated value to with being the Gaussian quadrature points over . can be reconstructed from in the following two steps. Firstly, we reconstruct from . To do this, we construct a polynomial such that

 1Δy∫yj+12yj−12Q(y)dy=1ΔxΔy∫Ii,ju(x,y)dxdy=¯ui,j, (2.15)

with belongs to a reconstruction stencil in the -direction as in the one-dimensional case. Then is a high order approximation to . We let to denote such reconstruction process in -direction. Secondly, we construct a polynomial such that

 1Δx∫xi+12xi−12P(x)dx=1Δx∫xi+12xi−12u(x,yig)dx, (2.16)

with belongs to a reconstruction stencil in the -direction as in the one-dimensional case. Then . Such 1D reconstruction process is denoted as . The 2D reconstructing procedure can be summarized as the following flowchart

Detailed information on the 2D reconstruction procedure is similar to those described in [1]. The 2D MPP flux limiter is applied in a similar fashion as those in [10, 9, 12]. Thus details are omitted for brevity.

###### Remark 2.2.

The proposed generalization of the parametrized flux limiter to convection-diffusion problems is rather straightforward. In comparison, it is much more difficult to generalize the polynomial rescaling approach in [17] to schemes with higher than third order accuracy for convection diffusion problems. The approach there relies on rewriting the updated cell average as a convex combination of some local quantities within the range ; this is more difficult to achieve with the diffusion terms [16, 19]. Moreover, the proposed flux limiter introduces very mild time step constraint to preserve both MPP and high order accuracy of the original FV RK scheme, see the next section for more discussions.

## 3 Theoretical properties

In this section, we provide accuracy analysis for the MPP flux limiter applied to the high order FV RK scheme solving pure convection problems. Specifically, we will prove that the proposed parametrized flux limiter as in equation (2.7) introduces a high order modification in space and time to the temporal integrated flux of the original scheme, assuming that the solution is smooth enough. A general proof on preservation of arbitrary high order accuracy will be provided for linear problems. Then by performing Taylor expansions around extrema, we prove that the modification from the proposed flux limiter is of at least third order, for FV RK schemes that are third order or higher in solving general nonlinear problems.

The entropy solution to a scalar convection problem

 ut+f(u)x=0,u(x,0)=u0(x). (3.1)

satisfies

 ddt∫xj+12xj−12u(x,t)dx=f(u(xj+12,t))−f(u(xj−12,t)). (3.2)

Integrating (3.2) over the time period , we have

 ¯uj(tn+1)=¯uj(tn)−λ(ˇfj+12−ˇfj−12), (3.3)

where and

 ¯uj(t)=1Δx∫xj+1/2xj−1/2u(x,t)dx,ˇfj−1/2=1Δt∫tn+1tnf(u(xj−1/2,t))dt. (3.4)

The entropy solution satisfies the maximum principle in the form of

 um≤¯uj(tn)−λ(ˇfj+12−ˇfj−12)≤uM. (3.5)

For schemes with order finite volume spatial discretization (2.6) and order RK time discretization, we assume

 |ˇfj+12−^Hrkj+12|=O(Δx2k+1+Δtp),∀j. (3.6)

Our analysis is in the sense of local truncation analysis assuming the difference between and is of high order (). Under a corresponding order reconstruction, the difference between the point values and is also of high order. In the following, we use them interchangeably when such high order difference allows.

For the MPP flux limiter, we only consider the maximum value part as in equation (2.9). The proof of equation (2.10) for the minimum value would be similar. We would like to prove that the difference between and in (2.7) is of high order in both space and time, that is

 |^Hrkj+12−~Hrkj+12|=O(Δx2k+1+Δtp),∀j. (3.7)

There are four cases of the maximum value part (2.9) outlined in the previous section. The estimate (3.7) can be easily checked for case (a) and (d) under the assumption (3.6) and the fact (3.5), see arguments in [12]. Below we will only discuss case (b), as the argument for case (c) would be similar.

First we give the following lemma:

###### Lemma 3.1.

Consider applying the MPP flux limiter (2.7) for the maximum value part (2.9) with case (b), to prove (3.7), it suffices to have

 |uM−(¯uj−λ(ˇfj+12−^hj−12))|=O(Δx2k+1+Δtp), (3.8)

if .

###### Proof.

For case (b), we are considering the case when

 Λ+12,Ij=ΓMj−λFj+12<1,

which is equivalent to , and

 ~Hrkj+12−^Hrkj+12=ΓMj+λFj+12−λ=uM−(¯uj−λ(^Hrkj+12−^hj−12))−λ,

which indicates that it suffices to have (3.8) to obtain (3.7) with the assumption (3.6). ∎

###### Theorem 3.2.

Assuming and , we have

 ¯uj(tn)−λ(ˇfj+12−f(¯uj−1(tn)))≤uM (3.9)

if is the entropy solution to (3.1) subject to initial data .

###### Proof.

Consider the problem (3.1) with a different initial condition at time level ,

 ~u(x,tn)=⎧⎨⎩u(x,tn)x≥xj−12,¯uj−1(tn)x

here is the exact solution of (3.1) at time level . Assuming is its entropy solution corresponding to the initial data , instantly we have

 ¯~uj(tn)=¯uj(tn). (3.11)

Since , we have

 f(~u(xj−12,t))=f(¯uj−1(tn)), (3.12)

for . Since , the characteristic starting from would not hit the side , therefore

 ~u(xj+12,t)=u(xj+12,t) (3.13)

for . Also since satisfies the maximum principle , we have

 ¯~un+1j=¯~unj−λ(ˇ~fj+12−ˇ~fj−12)≤uM,

where

 ˇ~fj−1/2=1Δt∫tn+1tnf(~u(xj−1/2,t))dt. (3.14)

Substituting (3.11), (3.12) and (3.13) into the above inequality, it follows that

 ¯u