Study of conservation and recurrence of Runge-Kutta discontinuous Galerkin schemes for Vlasov-Poisson systems

Study of conservation and recurrence of Runge-Kutta discontinuous Galerkin schemes for Vlasov-Poisson systems

Yingda Cheng Department of Mathematics, Michigan State University, East Lansing, MI 48824 U.S.A. ycheng@math.msu.edu    Irene M. Gamba Department of Mathematics and ICES, University of Texas at Austin, Austin, TX 78712 U.S.A. gamba@math.utexas.edu    Philip J. Morrison Department of Physics and Institute for Fusion Studies, University of Texas at Austin, Austin, TX 78712 U.S.A. morrison@physics.utexas.edu
July 6, 2019
Abstract

In this paper we consider Runge-Kutta discontinuous Galerkin (RKDG) schemes for Vlasov-Poisson systems that model collisionless plasmas. One-dimensional systems are emphasized. The RKDG method, originally devised to solve conservation laws, is seen to have excellent conservation properties, be readily designed for arbitrary order of accuracy, and capable of being used with a positivity-preserving limiter that guarantees positivity of the distribution functions. The RKDG solver for the Vlasov equation is the main focus, while the electric field is obtained through the classical representation by Green’s function for the Poisson equation. A rigorous study of recurrence of the DG methods is presented by Fourier analysis, and the impact of different polynomial spaces and the positivity-preserving limiters on the quality of the solutions is ascertained. Several benchmark test problems, such as Landau damping, the two-stream instability, and the KEEN (Kinetic Electro static Electron Nonlinear) wave, are given.

Keywords: Vlasov-Poisson, discontinuous Galerkin methods, recurrence, positivity-preserving, BGK mode, KEEN wave.

1 Introduction

The Vlasov-Poisson (VP) system is an important equation for modeling collisionless plasmas, one that possesses computational difficulties of more complete kinetic theories. Thus, it serves as an important test bed for algorithm development. The VP system describes the evolution of , the probability distribution function () for finding an electron (at position with velocity at time ) with a uniform background of fixed ions under a self-consistent electrostatic field. In particular, the non-dimensionalized VP system (with time scaled by the inverse plasma frequency and length scaled by the Debye length ) is given by

 ∂tf+v⋅∇xf−E⋅∇vf=0 Ω×(0,T] −ΔxΦ=1−∫Rnfdv Ωx×(0,T] (1) E=−∇xΦ Ωx×(0,T].

Here the domain , where can be either a finite domain or . The boundary conditions for the above systems are summarized as follows: as or . If is finite, then we can impose either inflow boundary conditions with on , where is the outward normal vector, or more simply impose periodic boundary conditions. For simplicity of discussion, in this paper, we will always assume periodicity in . Also, we add that when the VP system is applied to plasmas the total charge neutrality condition, , is imposed.

The following physical quantities associated with this system are related to its conservation properties:

 charge density ρ(x,t)=∫Rnf(x,v,t)dv, momentum density j(x,t)=∫Rnvf(x,v,t)dv, (2) kinetic energy density ξk(x,t)=12∫Rn|v|2f(x,v,t)dv, electrostatic energy density ξe(x,t)=12|E(x,t)|2.

Indeed, it is well-known that the VP system conserves the total electron charge , momentum , and energy . Moreover, any functional of the form is a constant of motion. In particular, this includes the -th order invariant and the entropy . Sometimes the functional is also called the enstrophy, and all of these invariants are called Casimir invariants (see, e.g., [34]).

The VP system has been studied extensively for the simulation of collisionless plasmas. Popular numerical approaches include Particle-In-Cell (PIC) methods [6, 24], Lagrangian particle methods [4, 17], semi-Lagrangian methods [8, 45], the WENO method coupled with Fourier collocation [58], finite volume (flux balance) methods [7, 18, 19], Fourier-Fourier spectral methods [27, 28], continuous finite element methods [49, 50], among many others. In this paper, we will focus on the discontinuous Galerkin (DG) method to solve the VP system. The original DG method was introduced by Reed and Hill [42] for neutron transport. Lesaint and Raviart [32] performed the first convergence study for the original DG method. Cockburn and Shu in a series of papers [14, 13, 12, 11, 15] developed the Runge-Kutta DG (RKDG) method for hyperbolic equations. The RKDG methods have been used to simulate the VP system in plasmas by Heath, Gamba, Morrison and Michler [23, 22] and for the gravitational infinite homogeneous stellar system by Cheng and Gamba [9]. Theoretical aspects about stability, accuracy and conservation of those methods are discussed in [22, 23] and more recently in [3, 2] for energy conserving schemes. Such methods have excellent conservation properties, can be readily designed for arbitrary order of accuracy, and have the potential for implementation on unstructured meshes with adaptivity. To ensure the positivity of the solution, one can use a maximum-principle-satisfying limiter that has been recently proposed by Zhang and Shu in [52] for conservation laws on cartesian meshes, and later extended on triangular meshes [56]. This limiter has been used to develop positivity-preserving schemes for compressible Euler [53, 55], shallow water equations [48], and Vlasov-Boltzmann transport equations [10]. It has also been employed recently in the framework of semi-Lagrangian DG methods [43, 41] for the VP system.

The scope of the present paper is as follows: we focus on a detailed study of the RKDG scheme for the Vlasov equation from both the numerical and analytical points of view. Since we are only considering one-dimensional problems, we use the classical representation of the solution by Green’s function to compute the Poisson equation; therefore, the electric field is explicitly given as a function of the numerical density. This removes all discretization error from the Poisson equation and lets us more accurately investigate our DG solver for the Vlasov equations. We rigorously study recurrence, which is an important numerical phenomenon that commonly appears with many solvers. We use Fourier analysis and obtain eigenvalues of the amplification matrix, and then investigate the impact of different polynomial spaces on the quality of the solution by examining conserved quantities as well as convergence to BGK states [5] for some choices of initial states. We consider benchmark test problems such as simulations of Landau damping phenomena for the linearized and nonlinear Vlasov Poisson systems, two-stream instability, and their long time BGK states and the formation of KEEN waves, both for the nonlinear system as well.

The remaining part of the paper is organized as follows: in Section 2, we describe the numerical algorithm and summarize its conservation properties. In Section 3, we study the recurrence phenomena that occurs for linear Vlasov type transport equations discretized by DG methods with various polynomial orders. Sections 4.1 and 4.2 are devoted to discussions of simulation results for the linearized and nonlinear VP system, respectively, for diverse choices of initial data and external drive potentials. We conclude with a few remarks in Section 5.

2 Numerical methods

In this section we first describe the proposed DG numerical algorithm and then discuss some of its basic conservation properties related to the quantities of (1). This is done for both the fully nonlinear VP system of (1) and the linearized VP system obtained by linearizing about the homogeneous equilibrium , with corresponding vanishing electric equilibrium field. Periodic boundary conditions in are assumed.

Thus, setting and expanding the system to first order approximation, the perturbation satisfies the Linear Vlasov-Poisson (LVP) system,

 ∂tf+v⋅∇xf=E⋅∇vfeq Ω×(0,T] ΔxΦ=∫Rnfdv Ωx×(0,T] (3) E=−∇xΦ Ωx×(0,T],

where has been replaced by to ease the notation. We find it convenient and efficient to intertwine the discussion of our algorithms for the VP and LVP systems. To avoid confusion in Section 2.1 we underline the words linear and nonlinear, signaling where discussions specific to each apply.

2.1 Numerical algorithm

For one-dimensional problems, we use a mesh that is a tensor product of grids in the and directions, because this simplifies the definitions of the mesh and polynomial space for the Poisson equation. Specifically, the domain is partitioned as follows:

 0=x12

where is chosen appropriately large to guarantee for . This is a reasonable assumption, because of the well-posedness of the one-dimensional Vlasov-Poisson system as indicated in [21]. The grid is defined as

 Ii,j=[xi−12,xi+12]×[vj−12,vj+12], Ji=[xi−1/2,xi+1/2],Kj=[vj−1/2,vj+1/2], i=1,…Nx,j=1,…Nv,

where and are center points of the cells.

We will make use of several approximation spaces. For the -domain, we consider the piecewise polynomial space of functions ,

 Zlh={ξ:ξ|Ji∈Pl(Ji),i=1,…Nx},

where is the space of polynomials in one dimension of degree up to . For the space, we consider two approximation spaces of functions ,

 Vlh={ϕ:ϕ|Ii,j∈Ql(Ii,j),i=1,…Nx,j=1,…Nv}

and

 Wlh={φ:φ|Ii,j∈Pl(Ii,j),i=1,…Nx,j=1,…Nv},

where denotes all polynomials of degree at most in and on , and . These two spaces are widely considered in the DG literature for multi-dimensional problems. A simple calculation shows that the number of degrees of freedom of is . For this is larger than the number of degrees of freedom of , which is .

First we describe the RKDG scheme for the linear Vlasov equation. We seek (or ), such that

 ∫Ii,j(fh)tφhdxdv−∫Ii,jvfh(φh)xdxdv+∫Kj(ˆvfhφ−h)i+12,vdv −∫Kj(ˆvfhφ+h)i−12,vdv=∫Ii,jEhf′eqφhdxdv (4)

holds for any test function (or ). Here and below, we use the following notations: is the discrete electric field, which is to be computed from Poisson’s equation, , , and is a numerical flux. We can assume that in each , has a single sign by properly partitioning the mesh. Then, the upwind flux is defined as

 (5)

The scheme for the nonlinear Vlasov equation is similar. We seek (or ), such that

 ∫Ii,j(fh)tφhdxdv−∫Ii,jvfh(φh)xdxdv+∫Kj(ˆvfhφ−h)i+12,vdv−∫Kj(ˆvfhφ+h)i−12,vdv +∫Ii,jEhfh(φh)vdxdv−∫Ji(ˆEhfhφ−h)x,j+12dx+∫Ji(ˆEhfhφ+h)x,j−12dx=0 (6)

holds for any test function (or ). The upwind flux for has been defined in (5) and the new flux needed for the nonlinear case is given by

 ˆEhfh={Ehf−hif∫JiEhdx≤0Ehf+hif∫JiEhdx>0. (7)

The above descriptions coupled with a suitable time discretization scheme, such as the TVD Runge-Kutta method [44], completes the RKDG methods. For example, suppose the semi-discrete schemes in (2.1) and (2.1) are written in the compact form

 ∫Ii,j(fh)tφhdxdv=Hi,j(fh,Eh,φh)

where for the linear Vlasov of (2.1)

 Hlini,j(fh,Eh,φh) = ∫Ii,jvfh(φh)xdxdv−∫Kj(ˆvfhφ−h)i+12,vdv + ∫Kj(ˆvfhφ+h)i−12,vdv+∫Ii,jEhf′eqφhdxdv,

while for the nonlinear Vlasov of (2.1)

 Hnonlini,j(fh,Eh,φh) = ∫Ii,jvfh(φh)xdxdv−∫Kj(ˆvfhφ−h)i+12,vdv+∫Kj(ˆvfhφ+h)i−12,vdv − ∫Ii,jEhfh(φh)vdxdv+∫Ji(ˆEhfhφ−h)x,j+12dx−∫Ji(ˆEhfhφ+h)x,j−12dx.

The third order TVD-RK method implements the following procedure for going from to :

 ∫Ii,jf(1)hφhdxdv=∫Ii,jfnhφhdxdv+△tHi,j(fnh,Enh,φh), ∫Ii,jf(2)hφhdxdv=34∫Ii,jfnhφhdxdv+14∫Ii,jf(1)hφhdxdv+△t4Hi,j(f(1)h,E(1)h,φh), (8) ∫Ii,jfn+1hφhdxdv=13∫Ii,jfnhφhdxdv+23∫Ii,jf(2)hφhdxdv+2△t3Hi,j(f(2)h,E(2)h,φh).

Poisson’s equation is used to obtain , , and . Beyond periodicity, we need to enforce some additional conditions to uniquely determine . For example, we can set one end of the spatial domain to ground, i.e. set . In the one-dimensional case, then the exact solution can be obtained if we enforce . For the nonlinear system we obtain

 Φh=∫x0∫s0ρh(z,t)dzds−x22−CEx,

where , and

 Eh=−(Φh)x=CE+x−∫x0ρh(z,t)dz, (9)

while for the linear system Poisson’s equation gives

 Φh=∫x0∫s0ρh(z,t)dzds−CEx,

where , and

 Eh=−(Φh)x=CE−∫x0ρh(z,t)dz. (10)

From (9) and (10), we see that if (or ), then ; hence, . The above approach uses the classical representation of the solution by Green’s function and will be referred to as the “exact” Poisson solver. It is valid only for the one-dimensional case. For higher dimensions, a suitable elliptic solver needs to be implemented, such as those discussed in [23]. Here we use the exact solver to entirely eliminate discretization error from Poisson’s equation and, thereby, spotlight the performance of the Vlasov solver.

Below we describe positivity-preserving limiters, as summarized in [54]. We only use such a limiter to enforce the positivity of for the nonlinear VP system. For each of the forward Euler steps of the Runge-Kutta time discretization, the following procedures are performed:

• On each cell , we evaluate , where , and denote the Gauss quadrature points, while denote the Gauss-Lobatto quadrature points.

• We compute , where is the cell average of on , and . This limiter has the effect of maintaining the cell average, while “squeezing” the function to be positive at all points in .

• Finally, we use instead of to compute the Euler forward step.

This completes the description of the numerical algorithm.

2.2 Scheme Conservation properties

In the following, we will briefly review and discuss some of the conservation properties of the above RKDG scheme for the nonlinear VP equations without the positivity-preserving limiter. Some of those results have been reported in [22] and [3].

Proposition 1

(charge conservation) For both the and spaces,

 ∑i,jHnonlini,j(fh,Eh,1)=Θ(fh,Eh,1)

which implies

 ∑i,j∫Ii,jfn+1hdxdv = ∑i,j∫Ii,jfnhdxdv+23△t(Θ(f(2)h,E(2)h,1) +14Θ(f(1)h,E(1)h,1)+14Θ(fnh,Enh,1))

for the fully discrete scheme (2.1). Here,

 Θ(fh,Eh,φh)=∑i∫Ji(ˆEhfhφh)x,Nv+12dx−∑i∫Ji(ˆEhfhφh)x,12dx

denotes the contribution from the phase space boundaries located at , and should be negligible if is chosen large enough.

Remark: Charge conservation (or mass conservation or probability normalization as it is sometimes called) states that the total charge will be preserved on the discrete level up to approximation errors associated with the phase space boundaries. The proof is straightforward and, therefore, omitted. The same conclusion can be proven for the linearized system. The positivity preserving limiter does not destroy this property because it keeps the cell averages unchanged.

Proposition 2

(Semi-discrete stability – enstrophy decay) For both the and spaces, . Hence,

 ddt∑i,j∫Ii,jf2hdxdv≤0.

The proof, for an arbitrary field , can be found in [10], Theorem 4, which applies directly here by setting the collisional form in that proof.

For the remainder of this subsection we will assume the DG solution satisfies the velocity boundary conditions . This is a reasonable assumption when is large enough. In particular, this will guarantee exact charge conservation, which implies that is constant in time . Therefore, using the definition of in (9), we can obtain periodicity in , i.e, . Without this assumption the propositions below contain multiple boundary terms and the proof becomes technical.

Proposition 3

(Momentum conservation) Assuming , for both the and spaces when , , which implies

 ∑i,j∫Ii,jfn+1hvdxdv=∑i,j∫Ii,jfnhvdxdv

for the fully discrete scheme.

Proof. Choosing in (2.1), we have

 ∑i,jHnonlini,j(fh,Eh,v) = ∑i,j(∫Ii,jvfh(v)xdxdv−∫Kj(ˆvfhv)i+12,vdv+∫Kj(ˆvfhv)i−12,vdv −∫Ii,jEhfhdxdv+∫Ji(ˆEhfhv)x,j+12dx−∫Ji(ˆEhfhv)x,j−12dx) = −∑i,j∫Ii,jEhfhdxdv=−∑i∫JiEhρhdx,

and using the exact Poisson solver together with the periodicity of and yields the following:

 ∑i∫JiEhρhdx = ∑i∫JiEh(ρh−1)dx+∑i∫JiEhdx = −∑i∫JiEh(Eh)xdx+∑i∫JiEhdx = −(E2h(L)−E2h(0))/2−Φ(L)+Φ(0)=0,

which completes the proof.

Remark: The above proof holds for the linearized system as well. Note, however, it relies on the use of the exact Poisson solver. For a full numerical DG Poisson solver, such as that developed in [23] for the discretization Poisson equation, exact momentum conservation remains true, as was proven in [22] by means of the DFUG method developed there for dealing with the discretized Poisson equation. However, the positivity-preserving limiter we use here will destroy this property because it was not designed to conserve the numerical momentum.

Proposition 4

(Semi-discrete total energy equality) Assuming , for both the and spaces when ,

 ddt(12∑i,j∫Ii,jfhv2dxdv+12∑i∫JiE2h(x)dx)=A(fh,Φh)=A(fh−f,Φh−P(Φh)),

where the operator and is any projection such that and at , for .

Proof. Choosing in (2.1) yields

 ddt∑i,j12∫Ii,jfhv2dxdv+∑i,j∫Ii,jEhfhvdxdv=0

and

 ddt∑i12∫JiE2h(x)dx = ∑i∫JiEh(Eh)tdx=−∑i∫Ji(Φh)x(Eh)tdx = ∑i∫JiΦh(Eh)xtdx=∑i∫JiΦh(1−ρh)tdx = −∑i∫JiΦh(ρh)tdx=−∑i,j∫Ii,jΦh(fh)tdxdv,

where in the second line, we have used the periodicity and continuity of and . Therefore, we have proven that

 ddt(12∑i,j∫Ii,jfhv2dxdv+12∑i∫JiE2h(x)dx)=A(fh,Φh).

On the other hand, upon choosing in (2.1) and using the periodicity and continuity of , we can verify that . The exact solution obviously satisfies from the continuity and periodicity of , and therefore we are done.

The above proof indicates that the variation in the total energy will be related to the error of the solution, , together with the projection error, . In [22, 23], error estimates for DG schemes with NIPG methods for the Poisson equations are provided for multiple dimensions. In [3], optimal accuracy of order for the semi-discrete scheme with spaces has been proven under certain regularity assumptions. We remark that in [3] conservation of the total numerical energy is proven when the Poisson equation is solved by a local DG method with a special flux. Unfortunately, no numerical simulations of linear Landau damping or of the nonlinear VP system, such as those done in [23] or in Section 4 of this present manuscript, have been performed up to this date by the scheme proposed in [3], so a comparison is not possible.

3 On recurrence

In this section we study recurrence, a numerical phenomenon that is known to occur in simulations of Vlasov-like equations. Its study is important because it provides information about the numerical accuracy of a scheme. Recurrence was observed in the semi-Lagrangian scheme of Cheng and Knorr [8], where a simple argument for its occurrence was provided. In this section, we carry out a detailed study of recurrence for the DG method.

We study recurrence of our algorithm applied to the linear advection equation on , since it is tractable and contains the basic recurrence mechanism; results for the full Vlasov system tend to be quite similar. The initial condition we consider is , and the equilibrium distribution is taken to be either the Maxwellian or Lorentzian distribution, viz.

 fM=1√2πe−v2/2orfL=1π1v2+1.

For the Maxwellian equilibrium, , we take , and for the Lorentzian equilibrium, , we take .

The exact solution for the advection equation is . Hence, a simple calculation shows for the Maxwellian distribution; similarly, for the Lorentzian, . Thus, we see how the density for each case should decay to zero. The failure of decay and the occurrence of recurrence noted for the semi-Lagrangian scheme of [8] stems from the finite resolution in the velocity space and, indeed, the recurrence time depends on the mesh size in .

To be specific, we repeat the definition of DG scheme for this equation, which amounts to (2.1) with set to zero: we find (or ) , such that

 ∫Ii,j(fh)tφhdxdv−∫Ii,jvfh(φh)xdxdv+∫Kj(ˆvfhφ−h)i+12,vdv−∫Kj(ˆvfhφ+h)i−12,vdv=0 (11)

holds for any test function (or ). Again is the upwind numerical flux of (5). In the analysis below, we always assume time to be continuous, because recurrence is mainly a phenomenon that comes from the spatial and velocity discretization.

3.1 The case of l=0

For the piecewise constant case, the DG method is equivalent to a simple first order finite volume scheme and we can derive rigorously the behavior for . Suppose we define on cell , and assume uniform grids in both directions. Moreover, we assume to be even for simplicity. With this assumption, the location of the cell center is . Now (11) simply becomes

 dfijdt+vjfij−fi−1,j△x=0ifvj≥0,dfijdt+vjfi+1,j−fij△x=0ifvj<0. (12)

The initial condition chosen is clearly equivalent to . We prove that the scheme above gives

 fij(t)=Re(Aeikxi+sjtfeq(vj)) (13)

where is given in (14) below.

Upon plugging (13) into (12), we have

 sjfij+vj1−e−ik△x△xfij=0ifvj≥0sjfij+vjeik△x−1△xfij=0ifvj<0.

Hence,

 sj=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩−vj1−e−ik△x△x=vjcos(k△x)−1△x−vjsin(k△x)△xiifvj≥0−vjeik△x−1△x=−vjcos(k△x)−1△x−vjsin(k△x)△xi%ifvj<0,

which can be summarized as

 sj=|vj|cos(k△x)−1△x−vjsin(k△x)△xi. (14)

Therefore, the real part of is always negative, this means the magnitude of will always damp, but because of the -dependence it does so at different rates for different cells. Since the density

 ρ(xi)=∑jfij△v=Re(∑jAeikxi+sjtfeq(vj))△v,

the density will damp at a rate between and . Another important fact is that recurring local maxima of the density will have a period that satisfies . If we define , then . When , , and this coincides with the recurrence time obtained in [8].

Next we compare the above theoretical prediction with numerical results. In all of the calculations below, we take , and the mesh size to be . In Figure 1, we display results for numerical runs using piecewise constant polynomials and time discretization using TVD-RK3. (We use the third order method to minimize the time discretization error.) We plot in Figure 1. First, we notice the pattern of has the expected periodic structure with damping for both Maxwellian and Lorentzian equilibria. For the Maxwellian distribution, a simple calculation yields . Similarly, with the formulas above, the smallest damping rate is , while the biggest is . For the Lorentzian distribution,