Conservative discontinuous Galerkin/Hermite Spectral Method for the Vlasov-Poisson System

04/06/2020 ∙ by Francis Filbet, et al. ∙ 0

We propose a class of conservative discontinuous Galerkin methods for the Vlasov-Poisson system written as a hyperbolic system using Hermite polynomials in the velocity variable. These schemes are designed to be systematically as accurate as one wants with provable conservation of mass and possibly total energy. Such properties in general are hard to achieve within other numerical method frameworks for simulating the Vlasov-Poisson system. The proposed scheme employs discontinuous Galerkin discretization for both the Vlasov and the Poisson equations, resulting in a consistent description of the distribution function and electric field. Numerical simulations are performed to verify the order of accuracy and conservation properties.



There are no comments yet.


page 16

page 17

page 20

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

One of the simplest model that is currently applied in plasma physics simulations is the Vlasov-Poisson system. Such a system describes the evolution of charged particles (electrons and ions) under the effects of a self-consistent electric field. For each species labelled , the unknown , depending on the time , the position , and the velocity , represents the distribution function of particles in phase space. This model can be used for the study of beam propagation or a collisionless plasma


where represents the single charge and represents the mass of one particle , whereas the electric field satisfies the Poisson equation


where is the vacuum permittivity. On the one hand, for a smooth and nonnegative initial data , the solution to (1.1) remains smooth and nonnegative for all . On the other hand, for any function , we have

which leads to the conservation of mass, norms, for and kinetic entropy,

We also get the conservation of momentum

and total energy

The numerical resolution of the Vlasov-Poisson system (1.1)-(1.2) is often performed by particle methods (PIC) which consist in approximating the plasma by a finite number of macro-particles. The trajectories of these particles are computed from the characteristic curves given by the Vlasov equation, whereas self-consistent fields are computed on a mesh of the physical space. This method allows to obtain satisfying results with a few number of particles (we refer to Birdsall-Langdon for more details [4]). However, it suffers from poor accuracy since the numerical noise only decreases in when the number of particles

is increased, which is not satisfying for some specific problems. Therefore, different approaches, discretizing the Vlasov equation on a mesh of phase space, have been proposed. Among them, the Fourier-Fourier transform is based on a Fast Fourier Transform of the distribution function in phase space, but this method is only valid for periodic boundary conditions


. Consequently, for non periodic boundary conditions, Gibbs oscillations form at the boundary of the grid and become source of spurious oscillations which propagate into the distribution function. Later semi-Lagrangian methods, which consist in computing the distribution function at each grid point by following the characteristic curves backward, were also used. To compute the origin of the characteristics, a high order interpolation method is needed. E. Sonnendrücker

proposed the cubic spline reconstruction which gives very good results [31], but the use of spline interpolation destroys the local character of the reconstruction. Nakamura and Yabe also presented the Cubic Interpolated Propagation (CIP) method based on the approximation of the gradients of the distribution function in order to use a Hermite interpolation [32]. This method is very expensive in memory computation since it needs the storage of , , and . Another scheme for the Vlasov equation is the Flux Corrected Transport (FCT) [5] or more recently finite volume methods [16, 17, 13, 30]: the basic idea is to compute the average of the Vlasov equation solution in each cell of the phase space grid by a conservative method. In the sixties, a variant has been proposed: rather to discretize the function in velocity space, Galerkin methods with a “small” finite set of orthogonal polynomials are used [1, 24]. More recently, in [33, 22] and [29], the authors have shown the merit to use rescaled orthogonal basis like the so-called scaled Hermite basis. In [29], the authors show that, according to some symmetry properties of the Hermite weights, the numerical method conserves certain integral quantities like the number of particles, the momentum, the energy or the

-norm. These properties lead to algorithms that can provide long time numerical stability while using a small set of basis functions. The aim of the present work is to apply these techniques in order to handle with the velocity space and design conservative discretization with a reduced degree of freedom.

For the space discretization, we adopt the point of view of discontinuous Galerkin methods. These methods are similar to finite element methods but use discontinuous polynomials and are particularly well suited to handle complicated boundaries which may arise in many practical applications. Moreover, their local construction endows the methods with good local conservation properties without sacrificing the order of accuracy. Furthermore, they are extremely flexible in handling -adaptivity, the boundary conditions are imposed weakly and the mass matrices are block-diagonal which results in very efficient time-stepping algorithms in the context of time-dependent problems, as it is the case here. Among the computational works, we mention the use of discontinuous Galerkin schemes for the Vlasov-Poisson system in [11, 20]; and Vlasov-Maxwell [10]. Several theoretical works have been performed to analyze a family of semi-discrete discontinuous Galerkin schemes for the Vlasov-Poisson system with periodic boundary conditions, for the one [3] and multi-dimensional cases [2] , but also for the Vlasov-Maxwell system [35]

. The authors show optimal error estimates for both the distribution function and the electrostatic field, and they study the conservation properties of the proposed schemes. Let us emphasize that these former works apply a discontinuous Galerkin method using a phase space mesh, whereas here we adopt this approach only in physical space. Furthermore, we propose to modify the fluxes to ensure conservation of mass, momentum and total energy. We prove these conservations for the semi-discrete and fully discrete cases (first and second order time discretizations).

This paper is organized as follows: in the first part (Section 2), we briefly reformulate the Vlasov equation using the Hermite basis in velocity, recalling some properties of the solution. Then, in Section 3, we present a variant of the discontinuous Galerkin method for the space discretization such that mass, momentum and energy are preserved for semi-discrete and fully discrete systems. In the last section (Section 4), we present numerical results for Landau damping, two stream instability and bump on the tail problems to illustrate accuracy and conservations of our discretization technique.

2. Hermite spectral form of the Vlasov equation

For simplicity, we now set all the physical constants to one and consider the one dimensional Vlasov-Poisson system for a single species with periodic boundary conditions in space, where the Vlasov equation reads


with , position and velocity . Moreover, the self-consistent electric field is determined by the Poisson equation


where the density is given by

and is a constant ensuring the quasi-neutrality condition of the plasma


2.1. Hermite basis

To discretize the velocity variable, our approach is based on spectral methods, where we expand the velocity part of the distribution function in basis functions (typically Fourier or Hermite), leading to a truncated set of moment equations for the expansion coefficients. Actually, spectral methods are commonly used to approximate the solution to partial differential equations

[6, 21] and in particular the Vlasov-Poisson system [28, 25, 14, 22, 29].

For instance, in their seminal work [28], M. Shoucri and G. Knorr applied Chebyshev polynomials, whereas in [25, 14] the authors used Fourier basis but these methods does not conserve neither momentum, nor total energy. Later, J. P. Holloway and J. W. Schumer [22, 29] applied Hermite functions. Indeed, the product of Hermite polynomials and a Gaussian, seems to be a natural choice for Maxwellian-type velocity profiles [19]. Moreover, compared with other classical polynomials, they possess a fairly simple expression for their derivatives and allow to get conservation of mass, momentum and total energy. Recently, these methods generate a new interest leading to new techniques to improve their efficiency [26, 7, 9, 27, 8, 12] .

Let us summarize the main characteristics of this approach. The solution is approximated by a finite sum which corresponds to a truncation of a series


where is the number of modes. The issue is then to determine a well-suited class of basis functions and to find the expansion coefficients . Here, we have chosen the following basis of normalized scaled asymmetrically weighted Hermite functions:


where corresponds to the scaling velocity and we have set , and for , has the following recursive relation


The Hermite basis (2.5) has the following properties


and is the Kronecker delta function. With those properties, one can substitute the expression (2.4) for into the Vlasov equation (2.1) and using the orthogonality property (2.7), it yields an evolution equation for each mode ,


Meanwhile, we first observe that the density satisfies

then the Poisson equation becomes


Observe that if we take in the expression (2.4), we get an infinite system (2.8)-(2.9) of equations for and , which is formally equivalent to the Vlasov-Poisson system (2.1)-(2.2).

2.2. Conservation properties

There are different kinds of Hermite approximations based on the choice of Hermite functions and weights. In [22, 29], J. W. Schumer and J. P. Holloway have discussed precisely the different choices of orthogonal Hermite basis functions depending on the form of their weight functions. If the basis is asymmetrically weighted, hence mass, momentum and total energy are preserved but not the norm. In the case of symmetrically weighted basis functions, particles number, total energy (for odd) or momentum (for even), as well as -norm are conserved and the numerical stability is shown to be better. Moreover, it has been shown in [33] and later in [26] that it is required to introduce a velocity scaling factor to make the method more accurate and stable. Here, the basis is asymmetrically weighted, hence we prove the following result for the truncated system (2.8)-(2.9).

Proposition 2.1.

For any , consider the distribution function given by the truncated series

where is a solution to the Vlasov-Poisson system written as (2.8)-(2.9). Then mass, momentum and total energy are conserved, that is,

and for the total energy,


We consider the first three equations of (2.8) given by


which will be related to the conservation of mass, momentum and energy. Thus, let us take a look at the conservation properties from (2.10), by considering periodic or compact support boundary conditions.

First the total mass is defined as

hence the conservation of mass directly comes from (2.10) by integrating the first equation with respect to .

Then we define the momentum as

due to the second equation (2.10) and using the Poisson equation (2.9), we get

so that the conservation of momentum is obtained.

Finally to prove the conservation of total energy , we compute the variation of the kinetic energy defined as

Thus, combining the first and third equations in (2.10) and integrating over , we get


On the other hand, multiplying the first equation in (2.10) by and integrating over , it yields

Using the Poisson equation (2.9), we have , hence

Combining this latter equality with (2.11), it yields the conservation of total energy

Remark 2.2.

The natural space associated to our asymmetrically weighted basis is

Unfortunately there is no estimate of the associated norm for the solution to the Vlasov-Poisson system (2.1)-(2.2), hence there is no warranty to get such estimate for the system (2.8)-(2.9). It is worth to mention that the symmetrically weighted basis with the norm would be a good choice, but it would affect the efficiency of the conservation algorithm which will proposed in the Section 3 (see Remark 3.2).

2.3. Filtering technique

Filtering is a common procedure to reduce the effects of the Gibbs phenomenon inherent to spectral methods [21]. The filter will consist in multiplying some spectral coefficients in (2.4) by a scaling factor in order to reduce the amplitude of high frequencies, for any ,

Here, we simply apply a filter, called Hou-Li’s filter, already proposed in [23] for Fourier spectral method, which reads

where the coefficient is chosen as . Actually this filter achieved great success in Fourier [23] and Hermite [12] spectral methods. We again refer to [12] for a detailled discussion on the application to the Vlasov-Poisson system using a Hermite spectral method.

Remark 2.3.

Observe that the filter is applied only when , hence the filtering process does not modify the coefficients , that is, conservations of mass, momentum and total energy are not affected.

3. Discontinuous Galerkin scheme

This section is devoted to the space and time discretizations of the Vlasov-Poisson system written in the form (2.8)-(2.9). We apply a local discontinuous Galerkin method for space and introduce a slight modification of the fluxes to preserve mass, momentum and local energy for this semi-discrete system. Then, we focus on the time discretization and prove the conservation for the fully-discrete case.

3.1. Semi-discrete discontinuous Galerkin scheme

In this section, we will define the discontinuous Galerkin (DG) scheme for the Vlasov-Poisson system with Hermite spectral basis in velocity (2.8)-(2.9).

Let us first introduce some notations and start with , a partition of . Here , . Each element is denoted as with its length , and .

Given any non-negative integer , we define a finite dimensional discrete piece-wise polynomial space

where the local space consists of polynomials of degree at most on the interval . We further denote the jump and the average of at defined as,

where . We also denote

From these notations, we apply a semi-discrete discontinuous Galerkin method for (2.8) as follows. We look for an approximation , such that for any , we have


where is defined as


The numerical flux in (3.2) could be taken as any consistent numerical fluxes. As an example, the global Lax-Friedrichs flux is used, which is defined as


with the numerical viscosity coefficient .

For the Poisson equation (2.9), we introduce the potential function , such that

Hence we get the one dimensional Poisson equation

for which we also define a local discontinuous Galerkin scheme: we look for a couple , such that for any and , we have


where the numerical fluxes and in (3.4) here are taken as in [3, 2],


with either being a positive constant or a constant multiplying by and is the mesh size.

3.2. Conservation properties for the semi-discrete scheme

The local discontinuous Galerkin method presented in the previous section naturally preserves the mass since the equation on only contains a convective term. Unfortunately, the momentum and total energy are not conserved due to the contribution of the source terms in (3.1). Therefore, we propose a slight modification of the fluxes for the unknowns in order to ensure the right conservations.

In this subsection, without any confusion, we will drop the subindex for simplicity. The aim is to modify the scheme (3.1)-(3.5) in order to ensure the conservation of mass, momentum and energy without deteriorating the order of accuracy. Thus, we follow the proof of Proposition 2.1 and consider the scheme (3.1) for combined with (3.4), by adding two residual terms to ensure the conservation of momentum and total energy


and the residual terms are given by


Apparently, these residues are given as products of two cell interface jump terms, which both are of high order, hence it should not deteriorate the order of accuracy. Therefore, we prove that this choice will ensure the discrete conservation properties of the scheme.

Proposition 3.1.

For any , we consider the solution to the system (3.1)-(3.5) combined with the residual terms (3.6)-(3.7). Then, discrete mass and momentum are conserved

and the discrete total energy, defined as

is also preserved with respect to time.


On the one hand, the conservation of mass easily follows by choosing in (3.1) for and summing over , hence since there is no source term we get

On the other hand, for the conservation of momentum, taking in the modified scheme (3.6) and summing over , we have


where the contribution of the residual term, using periodic boundary conditions, gives


Then we compute the contribution of the first source term on the right hand side of (3.8). Choosing in the first equation of (3.4), we get that

whereas taking in the second equation of (3.4), we have

hence adding the latter both equalities, it yields

Summing on , we obtain


Finally, using the definition of the flux in (3.5), we have

hence gathering (3.9) and (3.10) in (3.8), we get the conservation of momentum

Now let us show the conservation of total energy. Choosing in the modified scheme (3.6) and summing over , we get


where the contribution of the residual term, using periodic boundary conditions, gives


Let us evaluate the first term on the right hand side of (3.11). We consider the first equation of (3.4) and take ,

whereas in (3.1), we choose ,

Adding the latter two equalities, it yields

Summing on and using periodic boundary conditions, we get

To evaluate the right hand side of the latter term, we use that

and from the definition of the flux for the Poisson equation (3.5), we get


Now, it remains to evaluate the second term on the left hand side in (3.13), hence we compute the time derivative of the second equation in (3.4) and choose the test function , it gives

Then in the first equation of (3.4), we take and obtain