An efficient, conservative, time-implicit solver for the fully kinetic arbitrary-species 1D-2V Vlasov-Ampère system

https://doi.org/10.1016/j.jcp.2020.109686Get rights and content

Highlights

  • Algorithm for fully kinetic arbitrary-species (including electrons) 1D-2V Vlasov-Ampére system.

  • Fully implicit time integration through a HOLO-accelerated nonlinear iterative approach.

  • Nonlinear constraint functions ensure conservation of mass, momentum, and energy (and the discrete Gauss law).

  • Velocity-space adaptivity through analytic transformation of governing equations.

Abstract

We consider the solution of the fully kinetic (including electrons) Vlasov-Ampère system in a one-dimensional physical space and two-dimensional velocity space (1D-2V) for an arbitrary number of species with a time-implicit Eulerian algorithm. The problem of velocity-space meshing for disparate thermal and bulk velocities is dealt with by an adaptive coordinate transformation of the Vlasov equation for each species, which is then discretized, including the resulting inertial terms. Mass, momentum, and energy are conserved, and Gauss's law is enforced to within the nonlinear convergence tolerance of the iterative solver through a set of nonlinear constraint functions while permitting significant flexibility in choosing discretizations in time, configuration, and velocity space. We mitigate the temporal stiffness introduced by, e.g., the plasma frequency through the use of high-order/low-order (HOLO) acceleration of the iterative implicit solver. We present several numerical results for canonical problems of varying degrees of complexity, including the multiscale ion-acoustic shock wave problem, which demonstrate the efficacy, accuracy, and efficiency of the scheme.

Introduction

In recent years, it has become apparent that kinetic effects (i.e., particle long mean-free-path) can play a significant role in the evolution of high-energy-density (HED) plasma systems, such as inertial confinement fusion (ICF) capsule implosions [1], [2], [3], [4], [5]. To study these systems, radiation-hydrodynamic models are typically used; however, to resolve the long mean-free-path effects, it is necessary to employ a kinetic approach. Vlasov-Fokker-Planck codes, such as iFP [6] and FPion [7], have been developed with the goal of resolving ion kinetic effects in weakly collisional regimes with arbitrary Knudsen numbers. However, they continue to treat the electrons as a quasineutral, ambipolar fluid, including only an electron temperature equation. The fluid electron assumption neglects important kinetic plasma effects such as nonlocal electron heat transport, which may be necessary to correctly describe HED plasma phenomena such as shocks and ablation fronts in ICF implosions. This study proposes an efficient and accurate algorithmic solution for simulating the fully kinetic 1D-2V ion-electron system.

There have been attempts to account for nonlocal electron effects within fluid models. A common approach is limiting the electron heat flux to some fraction of the free-streaming flux [8]; another strategy is to spatially convolve the electron heat flux [9]. However, in order to describe electron kinetic effects accurately, it is necessary to solve a fully kinetic model. To this end, the Vlasov-Fokker-Planck(VFP)-Maxwell system of equations may be taken as a first-principles representation of a weakly-coupled plasma. In a one-dimensional spatial system, and assuming an electrostatic field response, this may – without loss of generality – be reduced to a 2D velocity space described by longitudinal (parallel) and perpendicular velocities. This leads to the 1D-2V VFP-Ampère (or VFP-Poisson) system. The 1D Ampère equation describes the evolution of the longitudinal electric field based on the moments of all species' velocity distribution functions, while the 1D-2V VFP equation describes the evolution of these distribution functions. To simplify the current development and presentation, we explore in this work only the collisionless aspects of the algorithm, i.e., the 1D-2V Vlasov-Ampère system. This in no way compromises our goal of developing a fully kinetic simulation tool for plasmas of arbitrary collisionality, as the numerical details of extending the approach described herein to include the Fokker-Planck (FP) collision operator may be considered independently. Interested readers are referred to e.g., Refs. [10], [11], [12] for a compatible fully implicit finite-difference approach.

To solve the Vlasov-Ampère system, there are a variety of possible approaches including temporally implicit or explicit applications of particle-in-cell (PIC) methods [13], [14], [15], [16], semi-Lagrangian grid-based methods [17], [18], [19], or Eulerian grid-based approaches [20], [21], [22]. In addition, it is possible to utilize semi-implicit strategies [23], which aim at treating only the stiff physics implicitly, as well as various asymptotic-preserving (AP) schemes [24], [25], which propose discrete formulations capable of capturing the correct asymptotic limit when stiff physics are stepped over. In this work, we apply a temporally fully implicit grid-based Eulerian approach. An implicit approach has significant advantages over explicit schemes, particularly for grid-based approaches where the advective Courant-Friedrichs-Lewy (CFL) time-step limit is determined by the fastest speed on the velocity grid. With a fully-implicit nonlinear iterative solver, highly temporally multiscale problems – in which the system dynamics are driven on time-scales much longer than the fastest supported time-scales – may become much more tractable. A classic collisionless example is the ion-acoustic shock wave, where the dynamic time-scale is roughly 100 times longer than the inverse electron plasma frequency, and may be 1000 times (or more) longer than the explicit time-scale based on the maximum grid velocity [20]. In addition, a fully implicit iterative scheme exhibits advantages over semi-implicit and AP schemes primarily in the straightforward application of existing strategies to enforce discrete conservation. In the current work, the high-order(HO)/low-order(LO) scheme (HOLO) is used so that the LO (moment) system of equations accelerates convergence of the HO (kinetic) system. The LO system consists of the moments of the plasma species' Vlasov equations coupled to Ampère's equation, while the HO system consists of the Vlasov equations. HOLO approaches have been used to solve a variety of systems [26], from neutron [27] and thermal radiation transport [28] to BGK gas-kinetics [29]. More importantly, the HOLO approach has been applied to the solution of collisionless [13], [20] and collisional [30] plasma systems. In an earlier study by Taitano and Chacón [20], the HOLO approach was used to accelerate Vlasov-Ampère convergence by using the LO system to efficiently evaluate the electric field with the higher-order moment closure provided by the HO system. In this work, we generalize this study both by applying a non-centered time integration scheme and by considering the adaptive velocity-space strategy proposed in Ref. [11].

When solving the Vlasov-Ampère system, a static velocity mesh may become inefficient in problems where the species temperatures and bulk velocities exhibit significant temporal and spatial variations. Specifically, the mesh must be large enough to capture both the shift in the bulk velocity and temperature evolution at the hottest location in space and time (the largest thermal speed), while maintaining a sufficient resolution for the coldest location (the smallest thermal speed) for all species. In contrast, a mesh that dynamically expands/contracts in space and time while shifting the center to track changes in their bulk velocities may efficiently resolve the hottest/coldest regions of each species. The present work applies an analytic transformation to the Vlasov equation for each species α, scaling it by a normalizing speed vα (which is a function of the thermal speed vth,α) and shifting the velocity space by an offset velocity u,α (which is a function of the bulk velocity u,α). This is similar to the approach described in Refs. [6], [7], [25], [31], [32].

To preserve the numerical accuracy of long-time simulations, we desire a discretization scheme for which the continuum symmetries of the governing equations (leading to mass, momentum, and energy conservation) are preserved in the discrete. Without a discrete conservation principle, long-term simulations may produce significant violations of the conservation properties due to accumulated discretization errors, which can manifest as numerical plasma heating or cooling [30], or a departure of the solution from the asymptotic hydrodynamic manifold. Indeed, as we shall demonstrate later, the failure to ensure discrete charge conservation (i.e., enforcing the discrete Gauss's law) leads to catastrophic failure in simulations, with significant departure from the correct solution. Further, in the case of the Vlasov-Fokker-Planck equation, Taitano et al. [6] showed that even neglecting to ensure discrete momentum and energy conservation relationships only in the Vlasov equation (while enforcing it in the Fokker-Planck collision operator) leads to extremely large numerical errors [O(1)]. Thus, discrete conservation is key for achieving high fidelity and accuracy.

Broadly, we will distinguish between two different strategies for achieving discrete conservation. The first, which we term a “passive” approach, relies on specifically chosen discretizations that “passively” preserve the structure of the governing equations. This is a general catch-all for the symplectic and Hamiltonian-preserving techniques for plasma physics systems described by Morrison [33]. Such techniques have also been called “structure-preserving” [33], [34], and are in general very effective at preserving invariants – for the Vlasov-Maxwell system, Shiroto et al. [34] demonstrated conservation errors on the order of machine precision. However, from our perspective, there are two significant shortcomings of this approach. The first is that they are generally only possible through central differencing schemes, which are not monotonic, positivity-preserving, or non-oscillatory (all of which are desirable properties). The second is that, in the case of the velocity-space transformed Vlasov equation used in this work, it is not readily apparent whether such “structure-preserving” discretizations are possible. The second strategy, which we term an “active” approach, is the strategy we employ in this work to achieve discrete conservation. To implement this approach, we introduce Lagrange-multiplier-like constraint functions into the discretized governing equations. These “nonlinear constraint functions” are defined so as to actively enforce certain continuum symmetries of the governing equations, ensuring conservation of e.g., mass, momentum, and energy in the discretized system. The primary benefit of the “active” strategy is that it permits a choice of arbitrary discretization schemes in both time (e.g., backward Euler, or BDF2) and phase space (e.g., SMART [35] or WENO [36]).

The present approach is similar to the strategies used in Refs. [20] and [31]; however, there are some important differences. In the previous Vlasov-Ampère implementation, the approach relied on a time-centered Crank-Nicolson integrator to achieve energy conservation. The current work uses a BDF2 temporal integration scheme (which is more appropriate for an eventual application to the collisional system), and can in principle be applied to an arbitrary temporal integration scheme. Further, the approach for enforcing discrete conservation with the velocity-space adaptivity (as in Ref. [31]) must be modified because of interaction with the additional constraint functions. Thus, the current work delivers an implicit algorithm for the fully kinetic, arbitrary-species 1D-2V Vlasov-Ampère system, which conserves mass, momentum, and energy to within nonlinear convergence tolerance. The algorithm is adaptive in the velocity space to ease meshing requirements due to temporal and spatial variations in the local bulk velocity and thermal speed of each species, while the nonlinear constraint functions that ensure conservation also allow substantial freedom of choice for temporal and advective discretizations.

The rest of this paper is organized as follows. Section 2 gives an overview of the governing equations for the Vlasov-Ampère/Poisson system in 1D-2V, its transformation in the velocity space, and the continuum-conservation symmetries of the system. Section 3 describes the discretization of the Vlasov system. Section 4 provides details of our strategy for ensuring the continuum conservation symmetries in the discretized system. In Sec. 5 we present our nonlinear iterative strategy for solving the discretized system implicitly in time using a HOLO acceleration scheme. We present numerical results highlighting the accuracy and performance of the algorithm for several canonical problems of varying difficulty in Sec. 6, and provide concluding remarks in Sec. 7.

Section snippets

Vlasov-Ampère system of equations

The Vlasov-Ampère/Poisson system may be regarded as a first-principles representation for a fully ionized electrostatic collisionless plasma. The governing equations are the Vlasov equations for each species α,tfα+x(vfα)+qαmαEv(fα)=0, which describe the evolution in phase space of distribution functions, fα, and Ampère's equation,ϵ0tE+αqαΓα=0,which describes the evolution of the electric field, E. In Eqs. (1) and (2), v is the particle velocity, and qα and mα are the particle charge and

Discretization of the transformed Vlasov equation

The Vlasov equation, Eq. (10), is discretized using finite differences in the transformed phase space as follows. The discrete cylindrical cell volume in the velocity space for a uniform velocity mesh isΩ˜j,k2πc,kΔcΔc, while the total discrete volume including the configuration space on a uniform mesh isΔV˜i,j,k=ΔxΩ˜j,k. The quantities Δc, Δc, and Δx are the mesh spacings for the parallel velocity, perpendicular velocity, and configuration space, respectively. The domains are defined to be

Discrete conservation strategy for charge, momentum, and energy

As we saw in Sec. 2.2, there are certain symmetries of the continuum equations that must be satisfied in order to conserve charge, momentum, and energy. For an arbitrary discretization, these symmetries will pose conflicting constraints. As a result, it will not generally be possible to satisfy all of them simultaneously unless we design our discretization such that it includes elements that ensure these properties. In the following development, we will present the discrete definitions of the

Solving the discretized Vlasov-Ampère system

To solve the discretized Vlasov-Ampère system, we use the high-order/low-order (HOLO) nonlinear acceleration iterative strategy [26]. HOLO accelerates the nonlinear convergence of the temporally implicit high-order (HO) Vlasov-Ampère system through a low-order (LO) representation, which efficiently exposes the stiff physics. The LO system is obtained from the velocity-space moments of the HO system. This approach has been successfully employed to solve the Vlasov–Ampère and

Numerical results

In this section, we demonstrate the accuracy, convergence, and conservation properties of the proposed numerical scheme. We do so using several canonical collisionless problems of increasing complexity, ranging from the linear Landau damping to an ion-acoustic shock wave. Unless specified otherwise, for all the problems presented we normalize the particle mass and charge to the electron mass, me, and proton charge, qp, while normalizing the temperature, density, velocity, and time to the

Conclusions

We have presented a fully conserving, adaptive algorithm for numerically integrating the 1D-2V multi-species Vlasov-Ampère system. The algorithm is applicable for the fully kinetic Vlasov system with an arbitrary number of species of arbitrary mass ratio. The velocity-space adaptivity scheme allows each species' velocity-space mesh to evolve according to variations in their bulk velocity and temperature resulting in very efficient velocity-space meshing. Conservation of the total mass,

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work was supported by the Thermonuclear Burn Initiative of the Advanced Simulation and Computing Program, used resources provided by the Institutional Computing Program at Los Alamos National Laboratory, and was performed under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory, managed by Triad National Security, LLC under contract 89233218CNA000001.

References (50)

  • William T. Taitano et al.

    Charge-and-energy conserving moment-based accelerator for a multi-species Vlasov-Fokker-Planck-Ampère system, part I: collisionless aspects

    J. Comput. Phys.

    (2015)
  • Richard B. Horne et al.

    A new code for electrostatic simulation by numerical integration of the Vlasov and Ampère equations using MacCormack's method

    J. Comput. Phys.

    (2001)
  • Yingda Cheng et al.

    Energy-conserving discontinuous Galerkin methods for the Vlasov-Maxwell system

    J. Comput. Phys.

    (2014)
  • Pierre Degond et al.

    Asymptotic-preserving methods and multiscale models for plasma physics

    J. Comput. Phys.

    (May 2017)
  • L. Chacón et al.

    Multiscale high-order/low-order (HOLO) algorithms and applications

    J. Comput. Phys.

    (2017)
  • William T. Taitano et al.

    Charge-and-energy conserving moment-based accelerator for a multi-species Vlasov-Fokker-Planck-Ampère system, part II: collisional aspects

    J. Comput. Phys.

    (2015)
  • Francis Filbet et al.

    A rescaling velocity method for dissipative kinetic equations: applications to granular media

    J. Comput. Phys.

    (Sep. 2013)
  • Takashi Shiroto et al.

    Quadratic conservative scheme for relativistic Vlasov-Maxwell system

    J. Comput. Phys.

    (2019)
  • Guang-Shan Jiang et al.

    Efficient implementation of weighted ENO schemes

    J. Comput. Phys.

    (1996)
  • P.J. Mardahl et al.

    Charge conservation in electromagnetic PIC codes; spectral comparison of Boris/DADI and Langdon-Marder methods

    Comput. Phys. Commun.

    (1997)
  • Barry Marder

    A method for incorporating Gauss' law into electromagnetic PIC codes

    J. Comput. Phys.

    (1987)
  • John Villasenor et al.

    Rigorous charge conservation for local electromagnetic field solvers

    Comput. Phys. Commun.

    (1992)
  • Yuxi Chen et al.

    Gauss's law satisfying energy-conserving semi-implicit particle-in-cell method

    J. Comput. Phys.

    (2019)
  • B.P. Leonard

    A stable and accurate convective modelling procedure based on quadratic upstream interpolation

    Comput. Methods Appl. Mech. Eng.

    (1979)
  • Jeffrey Willert et al.

    Leveraging Anderson acceleration for improved convergence of iterative solutions to transport systems

    J. Comput. Phys.

    (2014)
  • Cited by (9)

    • A charge-momentum-energy-conserving 1D3V hybrid Lagrangian–Eulerian method for Vlasov–Maxwell system

      2022, Journal of Computational Physics
      Citation Excerpt :

      Therefore, the balance of energy exchange is exactly maintained even if the numerical solution is polluted by truncation errors. This is the reason why the development of fully-conservative Eulerian scheme is relatively easy [12–17]. In this paper, we propose a hybrid Lagrangian–Eulerian (HLE) method for the Vlasov–Maxwell system; the time and configuration space are discretized by the Eulerian, whereas the velocity space is described by the Lagrangian in this method.

    View all citing articles on Scopus
    View full text