On the velocity space discretization for the Vlasov–Poisson system: Comparison between implicit Hermite spectral and Particle-in-Cell methods

https://doi.org/10.1016/j.cpc.2015.09.002Get rights and content

Abstract

We describe a spectral method for the numerical solution of the Vlasov–Poisson system where the velocity space is decomposed by means of an Hermite basis, and the configuration space is discretized via a Fourier decomposition. The novelty of our approach is an implicit time discretization that allows exact conservation of charge, momentum and energy. The computational efficiency and the cost-effectiveness of this method are compared to the fully-implicit PIC method recently introduced by Markidis and Lapenta (2011) and Chen et al. (2011). The following examples are discussed: Langmuir wave, Landau damping, ion-acoustic wave, two-stream instability. The Fourier–Hermite spectral method can achieve solutions that are several orders of magnitude more accurate at a fraction of the cost with respect to PIC.

Introduction

The Vlasov equation represents the cornerstone for the kinetic modeling of collisionless plasmas. It describes the time evolution of the distribution function of a population of charged particles that respond self-consistently to the effect of self- and externally induced electromagnetic fields. The numerical solution of the Vlasov equation for collisionless plasmas represents a very active area of research. The main challenge resides in the fact that the distribution function lives in a six dimensional phase space. Unarguably, the most widely used technique to solve the Vlasov equation is the Particle-In-Cell (PIC) method  [1], [2]. The main idea, originally developed in the 1950s, is to sample the distribution function in velocity space by means of a discrete number of super-particles. The electromagnetic field is represented on a grid in the computational domain, and the super-particles move through the grid according to the Lorentz force that is calculated by interpolation from the grid to the particles position. The particles interactions are therefore mediated by the grid, and in this way the number of operations per time step is reduced from Np2 (if the total force on a particle is calculated by summing the pair interaction with every other particle) to Np (with Np the number of super-particles). A comprehensive review of the PIC methodology can be found in  [3].

The primary shortcoming of PIC simulations is related to the numerical noise: even starting from an equilibrium configuration, the discrete nature of the particles generates instantaneous (i.e. within the first time step) unphysical perturbations that produce a ‘noise ground’ level in the fields, below which any physical signal is lost. The most obvious way to reduce the noise in PIC is by increasing the number of super-particles, i.e. refining the discretization in velocity space. The main problem is that while the simulation time scales roughly linearly with Np, the noise ground level decreases only as Np1/2, as typical of Monte Carlo methods, implying that problems that require a high signal-to-noise ratio require significant computational resources. Some examples that illustrate the impact of the particle noise on the efficiency of PIC codes are the recent studies in space plasmas on the acceleration of particles in the solar wind and the coupling between large scale turbulence and small scale (i.e. kinetic) dissipative effects  [4], [5], [6]. PIC codes have demonstrated very poor performances for these problems. In two-dimensional (2D) simulations,  [7] have shown that at least 10,000 particles per cell (equivalent to 106 in 3D) are needed to achieve a sufficiently large signal-to-noise ratio. For comparison, one of the biggest and most advanced PIC simulations in the world (of a magnetic reconnection problem) was done with 250 particles per cell and required the use of the highest performance, petascale computer available at that time  [8], [9]. We note that methods based on δf or re-mapping of the distribution function have been developed in the PIC community to reduce particle noise. This has led to some very interesting work (see for instance  [10], [11], [12]), but has not yet resulted in a commonly accepted denoising technique adopted by the PIC community.

The non-optimal scaling of the noise level with the number of particles in the PIC method has long been recognized in the computational plasma physics community  [13], [14], [15], [16], and it has perhaps been the main motivation for developing alternative ways of solving the Vlasov equation numerically. Nevertheless, despite its intrinsic noise, the PIC method still represents the preferred approach in the community, probably due to its simplicity and robustness, the relative ease to parallelize it, and the recent impressive advances in available computing power. However, one should always keep in mind the resource-intensive nature of PIC and the related cost. As a figure of merit on the cost of large scale PIC simulations, a single run using the global gyrokinetic PIC code developed within the SciDAC Gyrokinetic Particle Simulation Center (GPSC)  [17] for investigating the long-time evolution of turbulent transport in nuclear fusion devices was estimated, in 2008, to take approximately 25 million CPU hours, on 100,000+ cores  [18]. The monetary cost of such a simulation can be estimated at about $1,000,000  [19].

An alternative to the PIC method is represented by solving the Vlasov equation directly in phase-space (namely with a sixth dimensional computational grid in space and velocity), by means of a so-called Eulerian Vlasov code. This has been done, usually for problems with reduced dimensionality, with a variety of techniques such as high-order finite differences  [20], [21], finite elements  [22], [23], or finite volumes  [24], [25]. The reader is referred to  [26], [27] for a recent review of different methods. Notably, successful algorithms include the semi-Lagrangian schemes  [28], [29], [30], [22], [31], [32], [33], [34], [35], [36], [37], [38] and positive flux conservative methods  [39].

Yet another class of techniques is represented by spectral methods where the velocity space is projected onto a complete orthogonal basis  [40], [41], [42], [43], [44], [45].

Clearly, all these methods have advantages and shortcomings. Historically, the PIC method has been favored for large-scale, multidimensional problems which are still not doable with Vlasov solvers (mainly due to memory limitations)  [46]. On the other hand, Vlasov codes allow the study of fine scale details of the distribution functions that are typically inaccessible in PIC codes, due to the above mentioned noise problem. Vlasov codes, however, can suffer from serious numerical problems that can lead to a lack of discrete conservation properties, the occurrence of unphysical oscillations, and the generation of non-positive values in the distribution functions  [24].

In this paper, we focus on a spectral method where the distribution function is projected onto an Hermite basis in velocity space. This gives rise to a set of nonlinear, time- and space-dependent PDEs for the coefficients of the expansion. The use of Hermite functions to discretize the velocity space in the Vlasov equation dates back to the works of  [47], [40], [48], and  [41]. More recently this approach has been investigated by  [49], [50], [45] and, in the context of linearized equations, by  [51] and  [52].

The expansion of the distribution function using an Hermite basis in velocity space can be appealing for several reasons. First, the Hermite functions are a complete orthogonal basis with respect to a Gaussian weight function. As such, they are the optimal basis to represent Maxwellian-like distributions. In fact, an exact Maxwellian in velocity can be represented by only one term of the Hermite basis. Second, as already noted by  [47], the low order expansion coefficients are directly related to the low order moments of the distribution functions (i.e. density, mean velocity, energy, heat flux, etc.) that describe the plasma macroscopically and are usually the physical quantities of interest. As a consequence, the use of the Hermite basis allows a smooth transition from a fluid to a kinetic description by simply increasing the number of coefficients retained  [53]. This is an important and crucial feature of this method in approaching multi-scale problems and in assessing the importance of kinetic effects, which is not available in PIC or Vlasov methods that are forced to treat the full distribution function.

This paper builds largely upon the work of  [50]. In Ref.  [50] two different Hermite bases (so called ‘asymmetrically’ and ‘symmetrically’ weighted) and their properties were discussed, and a qualitative comparison between the Fourier–Hermite (FH) method and the PIC method was presented for simulations of Landau damping and bump-on-tail instability. Here we expand that work presenting an implicit time discretization, and quantitatively comparing the new scheme with an implicit PIC. Furthermore, we will show that with this discretization and periodic boundary conditions it is possible to conserve the total charge, momentum, and energy exactly in a finite time step (in contrast, as we discuss below, only charge and energy is conserved in the implicit PIC of  [54]).

Explicit PIC has often been criticized due to its stringent requirements on the choice of the time step and grid size, for numerical stability reasons. In fact, an explicit PIC code requires the resolution of the smallest time scale and the shortest spatial scale, even when the physics of interest only involves larger time/spatial scales. Of course, this generally translates into the requirement for large computational resources. Historically, the search for more efficient PIC schemes based on implicit time discretization dates back to the 1980s, when the implicit moment  [55], [56], [57] and the direct implicit  [58], [59] methods were introduced. Both methods rely on a linearization of the equations for the electromagnetic fields, and, thus, they should be more properly regarded as semi-implicit methods. Using these techniques, the numerical stability is greatly improved (with respect to explicit PIC), but energy is still not conserved exactly and, at each time step, there is a small inconsistency between the charge and current densities calculated from the particles and the one that is used for advancing the fields. Some authors have suggested that such limitations are responsible for the accumulation of numerical errors that preclude semi-implicit PIC simulations to run for long time intervals  [54], [60]. Recently, however, [61] and  [54] have formulated and successfully implemented a fully-implicit, one-dimensional PIC code (see also  [62] for an electromagnetic extension). The main feature that characterizes fully-implicit PIC is that particles and fields are advanced simultaneously through a Jacobian-Free Newton–Krylov (JFNK) solver, and converged nonlinearly within a certain tolerance. Moreover, by using the so-called particle enslavement, the nonlinear solver converges on a residual that does not contain particle positions and velocities. In the fully-implicit PIC method, energy is conserved within the level imposed by the nonlinear tolerance (i.e. almost exactly).  [54] have shown that by implementing a sub-stepping procedure that makes particles stop every time they cross a cell boundary, the charge is also conserved exactly. On the other hand, momentum is not conserved and must be monitored throughout the simulation (in contrast to semi-implicit PIC that exactly conserves the momentum, but not the energy).

When comparing different schemes, the key information is represented by the CPU time needed to obtain a solution with a certain accuracy (this metric was not considered in  [50]), in order to be able to assess which method must be preferred (and for which conditions). Hence, our comparisons between FH and PIC are presented in terms of computational efficiency and efficacy. The former is essentially represented by the CPU time required to achieve a certain accuracy in the solution, while the latter is a measure of the cost-effectiveness of the algorithm, i.e. how much an increased resolution and consequently a longer CPU time is paid off in terms of better accuracy  [63].

For the examples presented in this paper involving mostly near-equilibrium Maxwellian plasmas in one dimension, the FH method is several orders of magnitude more accurate than PIC for equal CPU time or, conversely, results with the same level of accuracy can be obtained in a fraction of CPU time.

The paper is organized as follows. In Section  2, we introduce the Fourier–Hermite expansion of the Vlasov–Poisson system in one dimension (1D). We use a fully-implicit discretization of the equation based on the Crank–Nicolson scheme, which leads to a non-linear system of equations that is solved numerically by means of a JFNK solver. We also show that the fully-implicit discretization ensures exact charge, momentum and energy conservation. Section  3 presents the comparison between FH and PIC for four standard cases: Langmuir wave, linear Landau damping, two-stream instability, and ion acoustic wave. We present conclusions in Section  4.

Section snippets

Numerical method: Fourier–Hermite basis

We study the Vlasov–Poisson system in the 2-dimensional phase space (x,v), where x denotes position and v velocity. The phase space is assumed to be periodic in x. In order to describe the method we specialize to the case of a plasma consisting of an electron and a singly charged ion species. The quantities are normalized as follows. Time is normalized on the electron plasma frequency ωpe, velocities on the electron thermal velocity vte=kTe/me, lengths on the electron Debye length λD, the

Filamentation and artificial collisionality

The development of increasingly smaller phase space structures in a collisionless plasma is very well known in plasma physics and typically referred to as the filamentation process. A classical and well-studied example where filamentation occurs is linear Landau damping, i.e., the damping in time of an initial electric field perturbation due to wave–particle resonances  [67]. In this case filamentation is due to the presence of the factor exp[ikvt] (k is the wavevector of the perturbation) in

Results

In this section we compare the performance of the Fourier–Hermite (FH) method with the PIC method, both implemented with fully-implicit time discretization. For the fully-implicit PIC method, we follow the approach of  [54], that nonlinearly solves the Ampere equation discretized in time with a Crank–Nicolson scheme. The current density is self-consistently calculated from the particles. The only difference between the fully-implicit PIC used in this paper and the one employed in  [54] is that

Conclusions

We have described a spectral method to numerically solve the Vlasov–Poisson equations that describe the evolution of a collisionless plasma. The velocity and configuration spaces are discretized by means of an Hermite and Fourier basis, respectively. In this paper, we have introduced an implicit Crank–Nicolson discretization in time, which is charge, momentum, and energy conserving. The simultaneous conservation of these three quantities is a very important property of this method, especially

Acknowledgments

The authors would like to thank L. Chacon, S. Markidis, V. Roytershteyn, X. Tang for useful conversations. This work was funded by the Laboratory Directed Research and Development program (LDRD) under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy by Los Alamos National Laboratory, operated by Los Alamos National Security LLC under contract DE-AC52-06NA25396.

References (80)

  • W. Deng et al.

    Comput. Phys. Comm.

    (2014)
  • R. Sydora

    J. Comput. Appl. Math.

    (1999)
  • J. Whealton et al.

    J. Comput. Phys.

    (1986)
  • N. Besse et al.

    J. Comput. Phys.

    (2003)
  • R. Heath et al.

    J. Comput. Phys.

    (2012)
  • M. Shoucri

    Commun. Nonlinear Sci. Numer. Simul.

    (2008)
  • M. Shoucri et al.

    J. Comput. Phys.

    (1974)
  • C.Z. Cheng et al.

    J. Comput. Phys.

    (1976)
  • E. Sonnendrücker et al.

    J. Comput. Phys.

    (1999)
  • K. Imadera et al.

    J. Comput. Phys.

    (2009)
  • N. Crouseilles et al.

    J. Comput. Phys.

    (2010)
  • J.-M. Qiu et al.

    J. Comput. Phys.

    (2010)
  • J.-M. Qiu et al.

    J. Comput. Phys.

    (2011)
  • J.A. Rossmanith et al.

    J. Comput. Phys.

    (2011)
  • F. Filbet et al.

    J. Comput. Phys.

    (2001)
  • A.J. Klimas

    J. Comput. Phys.

    (1983)
  • B. Eliasson

    J. Comput. Phys.

    (2003)
  • S. Le~Bourdiec et al.

    Comput. Phys. Comm.

    (2006)
  • A. Mangeney et al.

    J. Comput. Phys.

    (2002)
  • J. Schumer et al.

    J. Comput. Phys.

    (1998)
  • J. Vencels et al.

    Procedia Comput. Sci.

    (2015)
  • G. Chen et al.

    J. Comput. Phys.

    (2011)
  • R. Mason

    J. Comput. Phys.

    (1981)
  • J. Brackbill et al.

    J. Comput. Phys.

    (1982)
  • B. Cohen et al.

    J. Comput. Phys.

    (1982)
  • A.B. Langdon et al.

    J. Comput. Phys.

    (1983)
  • S. Markidis et al.

    J. Comput. Phys.

    (2011)
  • G. Chen et al.

    Comput. Phys. Comm.

    (2014)
  • G. Lapenta et al.

    J. Comput. Phys.

    (2006)
  • J. Canosa et al.

    J. Comput. Phys.

    (1974)
  • A.J. Klimas

    J. Comput. Phys.

    (1987)
  • A.J. Klimas

    J. Comput. Phys.

    (1994)
  • H.R. Lewis

    J. Comput. Phys.

    (1970)
  • E. Evstatiev et al.

    J. Comput. Phys.

    (2013)
  • C. Birdsall et al.

    Plasma Physics Via Computer Simulaition

    (2005)
  • R. Hockney et al.

    Computer Simulation Using Particles

    (2010)
  • J.P. Verboncoeur

    Plasma Phys. Control. Fusion

    (2005)
  • F. Valentini et al.

    Phys. Rev. Lett.

    (2008)
  • L. Matteini et al.

    Space Sci. Rev.

    (2012)
  • C.T. Haynes et al.

    Astrophys. J.

    (2014)
  • Cited by (0)

    View full text