Source-linked AI summary
An Energy- and Charge-conserving, Implicit, Electrostatic Particle-in-Cell Algorithm
Guangye Chen, Luis Chacón, Daniel C. Barnes
TL;DR
The paper addresses stability, conservation, accuracy, and efficiency limitations in implicit and explicit PIC plasma simulation. It develops a fully implicit, tightly converged Vlasov-Ampère solver with adaptive orbit integration and particle enslavement. The resulting scheme preserves charge and energy for arbitrary implicit time steps, maintains accuracy with large steps, and demonstrates moderate CPU gains in a 1D ion acoustic wave problem.
Problem
Explicit PIC faces CFL and, in some formulations, Debye-length resolution constraints, while earlier implicit approaches suffered accuracy pitfalls that motivated a more tightly converged formulation.
Method
The paper develops a fully implicit 1D electrostatic Vlasov-Ampère PIC solver using tight nonlinear Newton-Krylov convergence, orbit-averaged particle sub-stepping, and nonlinear particle enslavement.
Results
Large implicit time steps preserve accuracy in multiscale ion acoustic wave simulations, while the unpreconditioned Newton-Krylov implementation demonstrates moderate CPU gains up to 15 in 1D.
Takeaways & Limitations
The approach supports accurate, charge- and energy-conserving kinetic plasma simulations with very large implicit time steps and potential efficiency gains, particularly in multidimensional applications.
Takeaways & Limitations
The study is a proof of principle focused on a 1D electrostatic Vlasov-Ampère model, with future work needed for fully electromagnetic implementations and improved preconditioning strategies.
Abstract
from arXiv · showhide
This paper discusses a novel fully implicit formulation for a 1D electrostatic particle-in-cell (PIC) plasma simulation approach. Unlike earlier implicit electrostatic PIC approaches (which are based on a linearized Vlasov-Poisson formulation), ours is based on a nonlinearly converged Vlasov-Ampère (VA) model. By iterating particles and fields to a tight nonlinear convergence tolerance, the approach features superior stability and accuracy properties, avoiding most of the accuracy pitfalls in earlier implicit PIC implementations. In particular, the formulation is stable against temporal (CFL) and spatial (aliasing) instabilities. It is charge- and energy-conserving to numerical roundoff for arbitrary implicit time steps. While momentum is not exactly conserved, errors are kept small by an adaptive particle sub-stepping orbit integrator, which is instrumental to prevent particle tunneling. The VA model is orbit-averaged along particle orbits to enforce an energy conservation theorem with particle sub-stepping. As a result, very large time steps, constrained only by the dynamical time scale of interest, are possible without accuracy loss. Algorithmically, the approach features a Jacobian-free Newton-Krylov solver. A main development in this study is the nonlinear elimination of the new-time particle variables (positions and velocities). Such nonlinear elimination, which we term particle enslavement, results in a nonlinear formulation with memory requirements comparable to those of a fluid computation, and affords us substantial freedom in regards to the particle orbit integrator. Numerical examples are presented that demonstrate the advertised properties of the scheme. In particular, long-time ion acoustic wave simulations show that numerical accuracy does not degrade even with very large implicit time steps, and that significant CPU gains are possible.
1. Introduction
The paper develops a fully implicit 1D electrostatic PIC solver based on a tightly converged Vlasov-Ampère formulation, targeting the temporal, spatial, accuracy, and efficiency limitations of explicit and earlier implicit methods.
- Motivation: Explicit PIC is constrained by CFL time-step limits, while some momentum-conserving formulations also require resolving the Debye length to avoid finite-grid instabilities.These constraints make explicit PIC computationally demanding for multidimensional applications.
- Approach: The proposed solver iterates field-particle couplings to a tight nonlinear tolerance using a fully implicit Newton-Krylov formulation.The proof-of-principle implementation uses a one-dimensional electrostatic Vlasov-Ampère model.
- Conservation and stability: The formulation is stable against temporal and spatial instabilities and conserves charge and energy to numerical round-off for arbitrary implicit time steps.Momentum is not exactly conserved, but the paper reports small errors supported by adaptive particle sub-stepping.
- Conservation and stability: Orbit-averaging the time-centered Vlasov-Ampère equations enables exact energy conservation with particle sub-stepping and permits very large field-solver time steps without significant accuracy loss.The relevant time steps are constrained by the collective dynamical time scale of interest.
- Algorithmic design: Particle enslavement nonlinearly eliminates new-time particle positions and velocities, yielding fluid-comparable memory requirements and flexibility in the particle orbit integrator.This flexibility enables an adaptive charge-conserving mover.
- Results: Long-time ion acoustic wave simulations show that accurate orbit integration is essential, large implicit steps preserve accuracy, and moderate CPU gains are demonstrated with the unpreconditioned solver.The paper argues that larger gains are possible, especially in multiple dimensions when system size is much larger than the Debye length.
2. Electrostatic Vlasov-Ampère model
The paper formulates the 1D electrostatic plasma problem with the Vlasov-Ampère model, which is continuum-equivalent to Vlasov-Poisson but has different discrete conservation and stability properties.
- Vlasov-Poisson formulation: The electrostatic plasma can be described by Vlasov-Poisson equations involving species distribution functions, charges, masses, electric potential, electric field, charge density, and vacuum permittivity.The distribution functions evolve in phase space and determine the charge density.
- Vlasov-Ampère formulation: An alternate electrostatic formulation is derived from Ampère’s law using the charge continuity equation and the Vlasov-Poisson relations.The derivation introduces the plasma current density and integrates the 1D equation over the periodic domain.
- Vlasov-Ampère formulation: The spatially averaged current provides the solvability condition for the 1D electrostatic Ampère equation and preserves Galilean invariance.The condition appears through the term involving the average current density.
- Discrete properties: Although Vlasov-Poisson and Vlasov-Ampère are equivalent in the continuum, their discrete forms differ: VP is generally momentum- and charge-conserving, whereas VA can conserve energy and charge.The paper selects VA because it generally has better stability properties.
- Multidimensional extension: A multidimensional extension replaces the current correction with its longitudinal component so that the electric field obtained from Ampère’s law remains conservative.An alternative potential evolution requires inverting a Laplace operator to determine the potential.
3. Implicit particle-based discretization of the VA model
The implicit VA discretization combines time-centered particle-field updates with sub-stepped, orbit-averaged particle motion to preserve energy and charge while supporting large field time steps. Adaptive orbit integration and spatial filtering address orbit, tunneling, finite-grid, and aliasing errors, although exact momentum conservation is not enforced.
- 3.1. Exact energy conservation theorem: The Crank-Nicolson VA discretization is implicit, second-order accurate, unconditionally stable, and non-dissipative, enabling an exact energy-conservation theorem.The theorem depends on nonlinear consistency between field and particle updates and identical current-assignment and force-interpolation shape functions.
- 3.1. Exact energy conservation theorem: Exact total energy, including particle kinetic and electric-field energies, is conserved from time step n to n + 1.This conservation holds for the discretized electrostatic system when its coupled equations are updated nonlinearly consistently.
- 3.2. Particle sub-stepping and orbit-averaging: Particle sub-stepping uses small orbit steps satisfying kmaxvth∆τ < 1 while the field solver advances with a larger step ∆t ≥ ∆τ.This separates accurate fast-particle orbit integration from the larger time scale of collective dynamics.
- 3.2. Particle sub-stepping and orbit-averaging: Orbit-averaged current recovers energy conservation with particle sub-stepping, while the particle mover enforces exact charge conservation regardless of nonlinear convergence tolerance.The resulting VA implementation simultaneously conserves energy and charge exactly.
- 3.4. Particle tunneling and adaptive orbit integration: Adaptive orbit integration limits particle tunneling and controls momentum errors, but the energy- and charge-conserving mover does not enforce exact momentum conservation.Tunneling through potential barriers can send particles in the wrong direction and severely affect long-term accuracy.
- 3.5. Space filtering: Energy-conserving schemes can lose momentum through aliasing, so particle-shape and grid smoothing reduce short-wavelength signals while preserving energy and charge conservation.Binomial smoothing preserves the spatial average of current and retains both conservation properties.
- 3.6. Equivalence of VA and VP: Unlike explicit energy-conserving methods, this implicit approach conserves energy exactly for arbitrary ∆t and permits different particle and field time steps.That flexibility supports orbit integration choices for momentum-error control.
4. Formulation of nonlinear residual: nonlinear elimination/particle enslavement
The formulation eliminates new-time particle variables from the nonlinear residual, reducing the solve to a lower-dimensional field residual while preserving the nonlinear solution. This particle enslavement supports flexible particle integration and memory-efficient JFNK solution of the coupled system.
- Motivation: The original residual includes new-time field and particle quantities, giving iterative methods the full dimensionality of the kinetic problem and large storage requirements.Memory is identified as a principal bottleneck for massively parallel computers.
- Nonlinear elimination: Nonlinear elimination rewrites particle variables as X2 = f2(X1), producing a lower-dimensional residual G with the same nonlinear solution as F.The reduced residual requires less storage for Krylov and Newton iterations.
- Particle enslavement: For this PIC formulation, the electric field is X1 and particle positions and velocities are X2, which are solved from the particle equations for a given mesh field.The particle solve is local and uses substepping; Picard iteration with tolerance 10^-10 is sufficient for each sub-timestep.
- Particle enslavement: Particle enslavement removes particle quantities from the nonlinear residual, yielding a field-only residual while retaining particle data as auxiliary variables for residual evaluation.The formulation can be generalized to electromagnetic variables and retains flexibility in particle integration.
- JFNK solution: Krylov preconditioning can improve iteration convergence without changing the converged Jacobian-system solution, although this study demonstrates feasibility with an identity preconditioner.Effective preconditioning is identified as necessary for future multidimensional applications.
5. Numerical examples
The numerical tests show that the fully implicit scheme reproduces standard electrostatic behavior while remaining accurate for underresolved plasma dynamics and long ion-acoustic-wave simulations. Accuracy depends strongly on orbit integration and charge conservation, while large implicit steps can provide efficiency gains over explicit PIC.
- Langmuir wave: The implicit and explicit schemes agree closely for cold-plasma Langmuir waves at Δt = 0.1.The simulation uses L = 2π, Nx = 32, and Np = 2000, resolving the plasma period.
- Langmuir wave: At Δt = 2, the implicit scheme remains stable and accurately reproduces Landau damping despite an underresolved plasma frequency.The explicit scheme is unstable because it violates the CFL condition.
- Two-stream instability: The two-stream instability growth rates agree very well with the analytical value γ = 0.5 for both simple CN and ACC CN movers.Both simulations use Δt = 0.2.
- Two-stream instability: The ACC mover conserves charge exactly and substantially improves momentum conservation, while both CN movers conserve total energy commensurately with ϵt = 10−8.The simple CN mover produces significant charge-conservation errors, especially in the nonlinear regime.
- Ion acoustic wave: For ion acoustic waves, only the ACC mover accurately captures long-term dynamics with a large implicit step, while other movers accumulate substantial errors.Fixed-step CN fails when kmaxvtheΔt > 1; sub-stepping improves results but remains unsatisfactory, and momentum errors can produce nonphysical drifts.
- Ion acoustic wave: The implicit ACC solver matches the explicit solution only when the explicit timestep is reduced to Δt = 0.005, despite the implicit solver using a much larger timestep.This comparison supports the claim that comparable explicit accuracy requires a timestep much smaller than the explicit CFL limit.
- CPU gain and robustness: The ACC mover is the most robust to nonlinear-tolerance changes between 10−8 and 10−6, whereas simple and sub-stepped CN are strongly sensitive.At ϵt = 10−6, total-energy errors remain approximately 10−7 per timestep.
- CPU gain and robustness: Large CPU gains are possible for the implicit scheme, particularly when kλD ≪ 1, although bounded nonlinear iteration counts require suitable preconditioning.Developing such preconditioning strategies is identified as future work.
6. Conclusions
The study demonstrates a fully implicit 1D electrostatic PIC solver with exact energy and charge conservation, flexible particle integration, and fluid-like memory requirements. Multiscale ion acoustic wave experiments show accurate large-step operation and potential efficiency gains, while the study remains proof-of-principle.
- The 1D Vlasov-Ampère formulation exactly conserves energy through orbit averaging and charge through an adaptive charge-conserving sub-stepping mover.The formulation is developed as a proof of principle for electrostatic plasma simulation.
- Particle enslavement removes new-time particle variables from the nonlinear solve, enabling general orbit treatment and memory requirements comparable to fluid simulations.The particle orbit procedure is segregated during nonlinear residual evaluation, while particle states remain auxiliary variables.
- The unpreconditioned JFNK solver and tight nonlinear convergence provide the computational framework for the fully implicit formulation.A tight nonlinear tolerance is required to achieve exact energy conservation.
- Multiscale ion acoustic wave experiments demonstrate that large implicit time steps can be used without accuracy degradation.Comparable accuracy requires explicit approaches to use time steps much smaller than the CFL limit.
- The study argues that implicit methods can achieve large CPU gains, particularly in multidimensional simulations, but reports only moderate 1D gains up to 15.The reported unpreconditioned Newton-Krylov implementation achieved gains up to 15 in 1D.
- The results are explicitly limited to a proof-of-principle study, with future work targeting fully electromagnetic models, preconditioning-oriented moment-PIC formulations, and heterogeneous architectures.
Appendix A. Derivation of exact charge-conserving 1D particle mover
Appendix A derives a charge-conserving 1D particle mover using particle-based current and charge densities, B-spline shape functions, and Crank-Nicolson time discretization. Exact conservation depends on within-cell motion and low-order shape functions, requiring particles to stop at cell boundaries.
- The derivation defines current and charge density from particle charge, velocity, and B-spline shape functions.The current-density shape function is one order lower than the charge-density shape function.
- Crank-Nicolson time discretization is used to derive the discrete charge-conservation relation.
- Exact charge conservation holds when particles remain within a cell and the B-spline order satisfies m ≤2.The second-derivative terms cancel exactly under these conditions.
- Particles must stop at cell boundaries because the B-spline second derivative is only piecewise continuous within a cell.
Appendix B. Local error estimate of particle orbit integrator
Appendix B estimates local orbit-integration errors by comparing Euler prediction with a second-order time-centered correction. The derivation distinguishes first-order Euler accuracy from second-order Heun accuracy and obtains position and velocity error estimates.
- The appendix formulates the particle-orbit initial-value problem in terms of acceleration and initial velocity.
- Forward Euler provides the temporal predictor for the particle orbit discretization.
- Heun’s formula combines the Euler predictor with a second-order time-centered correction.
- Euler is first-order accurate, whereas Heun’s formula is second-order accurate.
- The appendix derives local error estimates separately for particle position and velocity and reports both in the main text.