High order maximum principle preserving finite volume method for convection dominated problems
Pei Yang ^{1}^{1}1Department of Mathematics, University of Houston, Houston, 77204. Email: peiyang@math.uh.edu. Tao Xiong ^{2}^{2}2Department of Mathematics, University of Houston, Houston, 77204. Email: txiong@math.uh.edu. JingMei Qiu ^{3}^{3}3Department of Mathematics, University of Houston, Houston, 77204. Email: jingqiu@math.uh.edu. The first, second and the third authors are supported by Air Force Office of Scientific Computing YIP grant FA9550120318, NSF grant DMS1217008. Zhengfu Xu ^{4}^{4}4Department of Mathematical Science, Michigan Technological University, Houghton, 49931. Email: zhengfux@mtu.edu. Supported by NSF grant DMS1316662.
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 RungeKutta 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 RungeKutta 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 convectiondominated 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 convectiondominated 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 convectiondominated problems [17, 16, 15, 9, 10, 12], positivity preserving schemes for compressible Euler and NavierStokes 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 nearvacuum state. To illustrate the purpose of the family of the MPP methods, we shall consider the solution to the following problem
(1.1) 
with . The solution to (1.1) satisfies the maximum principle, i.e.
(1.2) 
Within the high order finite volume (FV) RungeKutta (RK) weighted essentially nonoscillatory (WENO) framework, we would like to maintain a discrete form of (1.2):
(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 convectiondiffusion 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 convectiondominated 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 CourantFriedrichsLewy (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 semiLagrangian WENO method for solving the VlasovPoisson 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 multidimensional 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 convectiondominated 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 convectiondiffusion 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 ,
(2.1) 
with the computational cell and cell center defined as
(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 ,
(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 LaxFriedrichs flux
(2.4) 
Here , where is obtained by reconstructing a order polynomial whose averages agree with those in a leftbiased stencil ,
The reconstruction procedure for can be done similarly from a rightbiased 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 .

Reconstruct from the cell averages by constructing a cubic polynomial , such that
Then , . We use to denote such reconstruction procedure,
As a reference, the reconstruction formulas for are provided below,

Construct an interplant such that
Then let
Such procedure is denoted as
As a reference, we provide the formula for below
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 fivecell stencil with
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 noncompact 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
(2.5)  
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
(2.6) 
with and
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
(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
with . The goal of the procedures outlined below is to adjust , so that with the modified flux , the numerical solutions satisfy the maximum principle,
(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
The MPP property is satisfied with the modified flux (2.7) when the following inequalities are hold,
(2.9)  
(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:

If and , then

If and , then

If and , then

If and , then
Similarly, we can find a local pair of numbers such that for any
(2.10) holds. There are also four different cases:

If and , then

If and , then

If and , then

If and , then
Finally, the local limiter parameter at the cell boundary is defined as
(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) convectiondiffusion problems is straightforward. For example, we consider a 2D problem on a rectangular domain ,
(2.12) 
Without loss of generality, we consider a set of uniform mesh
with . A semidiscrete FV discretization of (2.12) gives
(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,
(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
(2.15) 
with belongs to a reconstruction stencil in the direction as in the onedimensional case. Then is a high order approximation to . We let to denote such reconstruction process in direction. Secondly, we construct a polynomial such that
(2.16) 
with belongs to a reconstruction stencil in the direction as in the onedimensional case. Then . Such 1D reconstruction process is denoted as . The 2D reconstructing procedure can be summarized as the following flowchart
(2.17) 
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 convectiondiffusion 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
(3.1) 
satisfies
(3.2) 
Integrating (3.2) over the time period , we have
(3.3) 
where and
(3.4) 
The entropy solution satisfies the maximum principle in the form of
(3.5) 
For schemes with order finite volume spatial discretization (2.6) and order RK time discretization, we assume
(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
(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.
Proof.
Theorem 3.2.
Proof.
Consider the problem (3.1) with a different initial condition at time level ,
(3.10) 
here is the exact solution of (3.1) at time level . Assuming is its entropy solution corresponding to the initial data , instantly we have
(3.11) 
Since , we have
(3.12) 
for . Since , the characteristic starting from would not hit the side , therefore
(3.13) 
for . Also since satisfies the maximum principle , we have
where
(3.14) 
Substituting (3.11), (3.12) and (3.13) into the above inequality, it follows that