1 Introduction

High Order Maximum Principle Preserving Semi-Lagrangian Finite Difference WENO schemes for the Vlasov Equation

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 and 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. Andrew Christlieb 444 Department of Mathematics, Michigan State University, East Lansing, MI, 48824. E-mail: andrewch@math.msu.edu

Abstract

In this paper, we propose the parametrized maximum principle preserving (MPP) flux limiter, originally developed in [Z. Xu, Math. Comp., (2013), in press], to the semi-Lagrangian finite difference weighted essentially non-oscillatory scheme for solving the Vlasov equation. The MPP flux limiter is proved to maintain up to fourth order accuracy for the semi-Lagrangian finite difference scheme without any time step restriction. Numerical studies on the Vlasov-Poisson system demonstrate the performance of the proposed method and its ability in preserving the positivity of the probability distribution function while maintaining the high order accuracy.

Keywords: Semi-Lagrangian method; Finite difference WENO scheme; Maximum principle preserving; Parametrized flux limiter; Vlasov equation

1 Introduction

In this paper, we will consider the Vlasov-Poisson (VP) system

 ∂tf+v⋅∇xf+E(t,x)⋅∇vf=0, (1.1) E(t,x)=−∇ϕ(t,x),−Δϕ(t,x)=ρ(t,x). (1.2)

on the domain , where and . It is an important system for modelling the collisionless plasma. The Vlasov equation (1.1) is a kinetic equation that describes the time evolution of the probability distribution function (PDF) of finding an electron at position with velocity at time . is the electric field and is a self-consistent electrostatic potential function described by the Poisson equation (1.2). The probability distribution function is coupled with the long range fields via the charge density , with uniformly distributed infinitely massive ions in the background. The total charge neutrality condition is imposed. The Vlasov equation (1.1) in the hyperbolic form enjoys the maximum principle preserving (MPP) and positivity preserving (PP) properties. That is, let

 fm=minx,v(f(x,v,0))≥0,fM=maxx,v(f(x,v,0)),

the solution at later time is still within the interval and stays positive.

The VP system has been extensively studied numerically. Popular numerical methods, besides the semi-Lagrangian (SL) approaches which will be reviewed in the next paragraph, include particle-in-cell (PIC) methods [12, 5, 20], Lagrangian particle methods [4, 13], weighted essentially non-oscillatory (WENO) coupled with Fourier collocation [42], Fourier-Fourier spectral methods [23, 22], finite volume methods [14, 16], continuous finite element methods [39, 38], Runge-Kutta discontinuous Galerkin methods [19, 3, 11, 9, 8], and methods in references therein.

In this paper, we focus on the SL method, which has been proposed and applied for a wide range of applications, e.g. in atmospheric modeling and simulations [32, 17], in capturing the moving interface via solving level set equations [33], in fluid dynamics [36] and in kinetic simulations [15, 34]. Compared to the Eulerian approach, the SL approach also has a fixed set of numerical meshes, but in each time step evolution, the information propagates along the characteristics in the Lagrangian fashion. The SL method can be designed as accurate as the Eulerian approach. In addition, the SL method is free of time step restriction by utilizing the characteristic method in the temporal direction. The design of SL approach requires the tracking of characteristics forward or backward in time, which could be challenging for nonlinear problems. To avoid such difficulty in solving the VP system, the dimensional splitting approach was proposed in [7]. Along this line, SL methods were developed in various settings. For example, in the finite difference framework, different interpolation strategies, such as the cubic spline interpolation [31], the cubic interpolated propagation [25], the high order WENO interpolation in a non-conservative form [6] and in a conservative form [26, 27, 28] are proposed. There have also been many work in designing the SL methods in the finite volume framework for the VP system [16] and for the guiding center Vlasov model [10], in the finite element discontinuous Galerkin (DG) framework [30, 29] and with a hybrid finite difference-finite element approach [18].

Our main interest in this paper is to develop a MPP, and in particular PP, SL finite difference WENO scheme for the VP system (1.1). The main challenge of designing MPP (and/or PP) schemes within the SL WENO framework is to maintain the designed high order of accuracy from the conservative WENO approximation. Meanwhile, it is desired that no additional restrictive CFL constraint will be introduced. There have been some research efforts on designing MPP and PP high order SL schemes. For example, in [16, 10], PP property is preserved in the finite volume framework. However, the PP preservation of the PDF is accompanied with sacrificing high order spatial accuracy. Recently, in [30, 29], the SL discontinuous Galerkin (DG) methods to solve the VP system is coupled with the MPP limiters, that were originally proposed by Zhang and Shu [40]. In these approaches, the limiters are applied to the reconstructed polynomials (for finite volume) or the representing polynomials (for DG). In general, the MPP (or PP) property together with the maintenance of high order accuracy is much more challenging to achieve in the finite difference framework than in the finite volume and finite element DG framework via limiting polynomials, see [41].

We propose to generalize the recently developed parametrized MPP flux limiter [37] to a conservative SL finite difference WENO method solving the VP system. The original parametrized flux limiter for 1D scalar conservation laws was later extended to the two-dimensional case [24]. Xiong et al. [35] proposed to apply the parametrized MPP flux limiter to the final Runge-Kutta (RK) stage only, with significantly improved time step restriction for maintenance of high order accuracy, leading to much reduced computational cost. It has also been proved in [35] that the parametrized MPP flux limiter can maintain up to third order accuracy both in space and in time for nonlinear scalar conservation laws.

In order to apply the parametrized MPP flux limiter in [37] for solving the VP system, we start with the dimensional splitting approach. The parametrized MPP flux limiter is proposed for the conservative high order SL finite difference WENO scheme [27, 26]. We mimic the proof in [35] to prove that the parametrized MPP flux limiter for the SL finite difference scheme solving linear advection equations can maintain up to fourth order accuracy without any time step constraint. We also apply the parametrized flux limiter proposed in [35] to the RK finite difference WENO approximation of the VP system without operator splitting. Through numerical studies on weak and strong Landau damping, two stream instabilities, and KEEN waves, we show that both methods perform very well with the designed MPP properties, while maintaining high order accuracy and mass conservation.

The rest of the paper is organized as follows. In Section 2, the SL finite difference WENO scheme is reviewed and the parametrized MPP flux limiting procedure is proposed for the SL WENO scheme. We also prove that the parametrized MPP flux limiter maintains up to fourth order accuracy without additional time step restriction. Numerical studies of the scheme are presented in Section 3. Conclusions are made in Section 4.

2 The parametrized MPP flux limiter for the SL WENO scheme

In this section, we will briefly describe the SL finite difference WENO scheme for solving Vlasov equation (1.1) in one dimension (both physical and phase spaces). We adopt the method in [28], which is in a conservative flux difference form and is suitable to be coupled with the newly developed parametrized MPP flux limiter [37].

For the one-dimensional problem, periodic boundary conditions are imposed in x-direction on a finite domain and in v-direction with a cut-off domain with chosen to be large enough to guarantee for . The domain is discretized by the the computational grid

 0=x12

Let and to be the middle point of each cell. The uniform mesh sizes in each direction are and .

Following [7], the Vlasov equation is dimensionally split to the following form

 ∂tf+v⋅∇xf=0, (2.1) ∂tf+E(t,x)⋅∇vf=0. (2.2)

Using the second order Strang splitting strategy, the numerical solution is updated from time level to time level by solving equation (2.1) for half a time step, then solving equation (2.2) by a full time step, followed by solving equation (2.1) for another half a time step. Note that each split equation (2.1) or (2.2) still preserves the MPP (and/or PP) property.

In the following, we will take the prototype 1D linear advection equation

 ut+aux=0 (2.3)

with constant , to present the SL finite difference scheme with the parametrized MPP flux limiters.

2.1 Review of SL finite difference WENO scheme

The SL finite difference WENO scheme proposed in [28] is based on integrating equation (2.3) in time over ,

 u(x,tn+1)=u(x,tn)−F(x)x, (2.4)

where

 F(x)=∫tn+1tnau(x,τ)dτ. (2.5)

By introducing a sliding average function

 F(x)=1Δx∫x+Δx2x−Δx2H(ζ)dζ, (2.6)

with

 F(x)x=1Δx(H(x+Δx2)−H(x−Δx2)), (2.7)

the evaluation of equation (2.4) at the grid point can be written in a conservative form

 un+1i=uni−1Δx(H(xi+12)−H(xi−12)), (2.8)

where is called the flux function. A SL finite difference WENO reconstruction is used to approximate the numerical flux function based on its cell averages

 ¯Hj=1Δx∫xj+12xj−12H(ζ)dζ,j=i−p,⋯,i+q, (2.9)

with

 ¯Hj=F(xj)=∫xjx⋆ju(ζ,tn)dζ,

where is the point tracing from the grid point along characteristics back to the time level . The last equality above is essential for the SL scheme, and is obtained via following characteristics. can be reconstructed from for each . In summary, the SL finite difference WENO scheme procedure in evolving equation (2.3) from to is as follows:

1. At each of the grid points at time level , say , trace the characteristic back to time level at .

2. Reconstruct from . We use to denote this reconstruction procedure

 R1[x⋆i,xi](uni−p1,⋯,uni+q1), (2.10)

to approximate , where indicates the stencil used in the reconstruction. indicates the reconstruction of .

3. Reconstruct from . We use to denote this reconstruction procedure

 Hi+12≐R2(¯Hi−p2,⋯,¯Hi+q2), (2.11)

to approximate . Here indicates the stencil used in the reconstruction.

4. Update the solution by

 un+1i=uni−1Δx(Hi+12−Hi−12), (2.12)

with numerical fluxes computed in the previous step.

When the reconstruction stencils in and above only involve one neighboring point value of the solution, then the scheme reduces to a first order monotone scheme when the time step is within CFL restriction. We let denote the first order flux. The proposed SL finite difference scheme can be designed to be of high order accuracy by including more points in the stencil for (the composition of and ), to reconstruct the numerical flux

 Hi+12=R2∘R1(uni−p,⋯,uni+q), (2.13)

where indicates the stencil used in the reconstruction process. The WENO mechanism can be introduced in the reconstruction procedures in order to realize a stable and non-oscillatory capture of fine scale structures. In the Appendix, we provide formulas to obtain the high order fluxes in (2.13) for the fifth order SL finite difference WENO scheme. For more details, we refer to [28].

For the case of large time step , if , is no longer inside , let to be the index such that and , we have

 hi+12 = i∑j=i⋆+1Δxunj+(xi⋆−x⋆i)uni⋆, (2.14) Hi+12 = i∑j=i⋆+1Δxunj+Hi⋆+12, (2.15)

where is reconstructed in the same fashion as (2.13), but replacing by . Similarly, if , let to be the index such that and , we have

 hi+12 = −i⋆∑j=i+1Δxunj+(xi⋆−x⋆i)uni⋆+1, (2.16) Hi+12 = −i⋆∑j=i+1Δxunj+Hi⋆+12. (2.17)

It is numerically demonstrated in [28] that the proposed high order SL WENO method works very well in Vlasov simulations with extra large time step evolution.

2.2 Parametrized MPP flux limiters

In this subsection, we propose a parametrized MPP flux limiter, as proposed in [37, 35], for the high order SL finite difference WENO scheme (2.12).

For simplicity, Let and as the minimum and maximum values of the initial condition. It has been known that the numerical solutions updated by (2.12) with the first order monotone flux satisfy the maximum principle. The MPP flux limiters are designed as modifying the high order numerical flux towards the first order monotone flux in the following way,

 ~Hi+12=θi+12(Hi+12−hi+12)+hi+12 (2.18)

where is the parameter to be designed to take advantage of the first order monotone flux in the MPP property and to take advantage of the high order flux in the high order accuracy.

Below is a detailed procedure of designing , in order to guarantee the MPP property of the numerical solutions, yet to choose ’s to be as close to as possible for high order accuracy. For each in limiting the numerical flux as in equation (2.18), we are looking for the upper bounds and , such that, for all

 θi+12∈[0,Λ+12,Ii]∩[0,Λ−12,Ii+1], (2.19)

the updated numerical solutions and by the SL WENO scheme (2.12) with the modified numerical fluxes (2.18) are within . Let

 ΓMi=uM−uni+1Δx(hi+12−hi−12),Γmi=um−uni+1Δx(hi+12−hi−12),

the MPP property of the first order monotone flux guarantees

 ΓMi≥0,Γmi≤0. (2.20)

To ensure with as in equation (2.18), it is sufficient to have

 1Δxθi−12(Hi−12−hi−12)−1Δxθi+12(Hi+12−hi+12)−ΓMi ≤ 0, (2.21) 1Δxθi−12(Hi−12−hi−12)−1Δxθi+12(Hi+12−hi+12)−Γmi ≥ 0. (2.22)

The linear decoupling of , subject to the constraints (2.21) and (2.22), is achieved via a case-by-case discussion based on the sign of

 Fi±12≐1Δx(Hi±12−hi±12),

as outlined below.

1. Assume

 θi−12∈[0,ΛM−12,Ii],θi+12∈[0,ΛM+12,Ii],

where and are the upper bounds of , subject to the constraint (2.21).

1. If and ,

 (ΛM−12,Ii,ΛM+12,Ii)=(1,1).
2. If and ,

 (ΛM−12,Ii,ΛM+12,Ii)=(1,min(1,ΓMi−Fi+12)).
3. If and ,

 (ΛM−12,Ii,ΛM+12,Ii)=(min(1,ΓMiFi−12),1).
4. If and ,

• If equation (2.21) is satisfied with , then

 (ΛM−12,Ii,ΛM+12,Ii)=(1,1).
• Otherwise,

 (ΛM−12,Ii,ΛM+12,Ii)=(ΓMiFi−12−Fi+12,ΓMiFi−12−Fi+12).
2. Similarly assume

 θi−12∈[0,Λm−12,Ii],θi+12∈[0,Λm+12,Ii],

where and are the upper bounds of , subject to the constraint (2.22).

1. If and ,

 (Λm−12,Ii,Λm+12,Ii)=(1,1).
2. If and ,

 (Λm−12,Ii,Λm+12,Ii)=(1,min(1,Γmi−Fi+12)).
3. If and ,

 (Λm−12,Ii,Λm+12,Ii)=(min(1,ΓmiFi−12),1).
4. If and ,

• If equation (2.22) is satisfied with , then

 (Λm−12,Ii,Λm+12,Ii)=(1,1).
• Otherwise,

 (Λm−12,Ii,Λm+12,Ii)=(ΓmiFi−12−Fi+12,ΓmiFi−12−Fi+12).

Notice that the range of (2.19) is determined by the need to ensure both the upper bound (2.21) and the lower bound (2.22) of numerical solutions in both cell and . Thus the locally defined limiting parameter is given as

 θi+12=min(Λ+12,Ii,Λ−12,Ii+1), (2.23)

with and . The modified flux in equation (2.18) with the designed above ensures the maximum principle. Such modified flux is consistent and preserves the maximum principle by its design, since it is a convex combination () of a high order flux and the first order flux . The mass conservation property is preserved, due to the flux difference form of the scheme.

Remark 2.1.

(Machine zero) Due to rounding floating point arithmetic, equations (2.20), (2.21) and (2.22) are enforced only at the level of machine precision. In some extreme cases, numerical solutions may go out of bound, but only at the level of machine precision.

In [35], it was proved that the MPP flux limiter for the RK finite difference scheme can maintain up to third order accuracy both in space and in time for the nonlinear scalar equation . In Theorem 2.2 below, we prove the proposed MPP limiter for the SL finite difference scheme can maintain up to fourth order accuracy for solving the linear advection equation (2.3) without any time step restriction. The generalization of the proof to maintain up to fifth order accuracy would be much more technical and will be investigated in our future work.

Theorem 2.2.

Consider solving the linear advection equation (2.3) using the finite difference SL method (2.12) with a fourth order reconstruction procedure. Assume the global error,

 enj=|unj−u(xj,tn)|=O(Δx4),∀n,j. (2.24)

Consider applying the parametrized MPP flux limiter to the numerical fluxes, then

 ∣∣∣1Δx(Hj+12−~Hj+12)∣∣∣=O(Δx4),∀j, (2.25)

without any time step restriction.

Proof.

Without loss of generality, we consider with . The case of can be reduced to by shifting the numerical solution with several whole grid points. Without specifying, in the following, we use instead of and use instead of . Since the difference between and is of high order due to assumption (2.24), we use and interchangeably when such high order difference allows. The high order flux is taken to be (2.13) with a 4th order finite difference reconstruction, the first order monotone flux is .

We mimic the proof in [35] and only consider the most challenging case: case (b) for the maximum value part. The proof for other cases would be similar to that in [35].

We consider case (b) when

 Λ+12,Ij=ΓMj−Fj+12<1. (2.26)

To prove (2.25), it suffices to prove

 uM−(uj−ξ(^Hj+12−uj−1))=O(Δx4), (2.27)

if . For the SL method, we have

 ^Hj+12=1ΔtHj+12=1Δt∫tn+1tnH(xj+12,τ)dτ, (2.28)

where is defined by

 u(x,t)=1Δx∫x+Δx2x−Δx2H(ξ,t)dξ. (2.29)

follows the same characteristics as the linear advection equation, hence

 ^Hj+12=1Δt∫Δt0H(xj+12−τ,tn)dτ. (2.30)

Let . We approximate by a 4th order reconstructed polynomial from the solution based on the stencil , with which

 I1 = uj−ξ(1Δt∫Δt0H(xj+12−τ,tn)dτ−uj−1) (2.31) = 124(26uj−10uj−1+2uj−2+6uj+1)+ξ24(9uj+3uj−1−uj−2−11uj+1) +ξ224(−14uj+10uj−1−2uj−2+6uj+1)+ξ324(3uj−3uj−1+uj−2−uj+1).

We first consider the case when the maximum point at , with , and . is the dependent domain for the exact solution when . We perform Taylor expansions around up to th order and obtain

 uj+1 = uM+u′M(xj−xM+Δx)+u′′M(xj−xM+Δx)22+u′′′M(xj−xM+Δx)36+O(Δx4), uj = uM+u′M(xj−xM)+u′′M(xj−xM)22+u′′′M(xj−xM)36+O(Δx4), uj−1 = uM+u′M(xj−xM−Δx)+u′′M(xj−xM−Δx)22+u′′′M(xj−xM−Δx)36+O(Δx4), uj−2 = uM+u′M(xj−xM−2Δx)+u′′M(xj−xM−2Δx)22+u′′′M(xj−xM−2Δx)36+O(Δx4).

Since , we can write to be

 I1=uM+u′′MΔx22R2+u′′′MΔx36R3+O(Δx4), (2.32)

where

 R2 = 56ξ+12ξ2−13ξ3+(ξ2−3ξ)z+z2, R3 = 14(−4ξ+ξ2−2ξ3+ξ4)+14(10ξ+6ξ2−4ξ3)z+14(−18ξ+6ξ2)z2+z3,

and .

We consider a quantity denoted by , which can be written as follows

 I2 = β1u(xM)+(1−β1)u(xM+c1Δx) (2.33) = uM+u′′MΔx22(1−β1)c21+u′′′MΔx36(1−β1)c31+O(Δx4),

with parameters and to be determined. If , .

Since , if we can find , such that

 (1−β1)c21≤R2,(1−β1)c31=R3, (2.34)

we would have

 |I1−uM|=O(Δx4), (2.35)

when (2.32) is compared with (2.33) under the assumption (2.26), thus proving (2.27) in this case. In the following, to determine the parameters and satisfying (2.34), we first need , that is

 1−β1≤min0≤z≤1(R32R23),with0≤ξ≤1. (2.36)

Using Mathematica, we have

 (2.37)

and

 6481≤54−108ξ+72ξ2−16ξ354−81ξ+27ξ2≤1, if0≤ξ≤12, (2.38) 0≤12−32√32√2−7ξ+9ξ2−5ξ3+ξ427−32ξ+9ξ2≤1, if0≤ξ≤1. (2.39)

Thus, we can choose . Let . Below, we verify that is bounded.

If , from (2.38), we have

 |c1| ≤ ⎛⎜⎝max0≤ξ≤1,0≤z≤1|R3|min0≤ξ≤12(1−β1)⎞⎟⎠1/3≤(8164max0≤ξ≤,0≤z≤1|R3|)1/3≤3.

Otherwise if , we have

 R31−β1=Λ(ξ)(14(−4+ξ−2ξ2+ξ3)+14(10+6ξ−4ξ2)z+14(−18+6ξ)z2+z3ξ), (2.40)

with

 216125≤Λ(ξ)=ξ12−32√32√2−7ξ+9ξ2−5ξ3+ξ427−32ξ+9ξ2≤2, (2.41) 14≤|14(−4+ξ−2ξ2+ξ3)+14(10+