PaperlayerView exact revision on arXiv

Source-linked AI summary

Discontinuous Galerkin algorithms for fully kinetic plasmas

J. Juno, A. Hakim, J. TenBarge, E. Shi, W. Dorland

arXiv:1705.05407v3physics.plasm-phastro-ph.SRcs.NAmath.NAphysics.comp-phphysics.space-ph

TL;DR

Full kinetic plasma simulations need a Vlasov-Maxwell solver, but particle-in-cell methods can introduce counting noise. The paper develops an explicit discontinuous Galerkin solver with a reduced basis and reports benchmarks spanning simple waves and shocks through five-dimensional turbulence.

  • Problem

    Full Vlasov-Maxwell simulations are important for low-signal kinetic plasma problems, while particle-in-cell counting noise can complicate analysis and kinetic reductions are not valid in every context.

    Source: §1 · p. 2

  • Method

    The paper discretizes the multi-species Vlasov-Maxwell system with a high-order discontinuous Galerkin method, explicit SSP-RK3 time stepping, and a reduced high-order basis to reduce cost.

    Source: §1 · p. 2, §3.4 · p. 14, §4 · pp. 14–15

  • Results

    The solver is verified against electrostatic and fully electromagnetic benchmarks ranging from one configuration and velocity dimension to five-dimensional turbulence, and the Orszag-Tang distribution shows greater than eight orders of magnitude of signal.

    Source: §6 · pp. 43–44

  • Takeaways & Limitations

    The benchmarks demonstrate a robust DG Vlasov-Maxwell solver and a noise-free distribution function suited to diagnostics such as field-particle correlations.

    Source: §6 · pp. 43–44

    Velocity-space boundaries remain unphysical, and particle acceleration or heating can bring the plasma into contact with the velocity-space boundary.

    Source: §6 · p. 44

Abstract

from arXiv · show

We present a new algorithm for the discretization of the Vlasov-Maxwell system of equations for the study of plasmas in the kinetic regime. Using the discontinuous Galerkin finite element method for the spatial discretization, we obtain a high order accurate solution for the plasma's distribution function. Time stepping for the distribution function is done explicitly with a third order strong-stability preserving Runge-Kutta method. Since the Vlasov equation in the Vlasov-Maxwell system is a high dimensional transport equation, up to six dimensions plus time, we take special care to note various features we have implemented to reduce the cost while maintaining the integrity of the solution, including the use of a reduced high-order basis set. A series of benchmarks, from simple wave and shock calculations, to a five dimensional turbulence simulation, are presented to verify the efficacy of our set of numerical methods, as well as demonstrate the power of the implemented features.

1. Introduction

Kinetic plasmas can require a full Vlasov-Maxwell description, but reduced models have limited scope and particle-in-cell simulations introduce counting noise. This paper adds a discontinuous Galerkin Vlasov-Maxwell solver to Gkeyll and outlines its numerical approach and benchmarks.

1. Introduction

  • Fluid models may need sub-grid closures, while kinetic reductions such as gyrokinetics are not valid in all contexts, including some large-amplitude waves above the cyclotron frequency.The introduction identifies PIC as useful where a tractable, physically reasonable asymptotic reduction does not exist.Source: §1 · p. 2
  • The paper adds a multi-species Vlasov-Maxwell solver to Gkeyll, using discontinuous Galerkin spatial discretization and an explicit strong-stability-preserving Runge-Kutta time update.The introduction highlights high-order accuracy and data locality as advantages for accuracy and parallelization.Source: §1 · p. 2
  • For low-signal problems requiring full kinetic physics, directly solving Vlasov can provide a noise-free solution, unlike particle-in-cell methods with counting noise that decreases as 1/N.Delta-f PIC noise mitigation can break down when the distribution function deviates substantially from its initial value.Source: §1 · p. 2
  • The paper proceeds from kinetic equations and conservation properties to implementation details, then evaluates the algorithm with scaling studies and high-dimensional benchmarks.The introduction names a five-dimensional turbulence simulation among the benchmark studies.Source: §1 · p. 3

2. The Vlasov-Maxwell system

The paper formulates the non-relativistic, collisionless Vlasov-Maxwell system and reviews conservation properties of its continuous equations. The distribution evolves in phase space while its moments determine charge and current that couple it self-consistently to electromagnetic fields.

2.1. Basic equations

  • The Vlasov equation evolves each species' distribution in configuration and velocity space under the Lorentz force, with the non-relativistic and collisionless limits assumed.The Lorentz force depends on the electric and magnetic fields; the paper notes that collisions or relativistic effects require extensions to the stated system.Source: §2.1 · p. 3
  • Charge and current densities are moments of the distribution function, and couple the Vlasov equation to Maxwell's equations so the fields evolve self-consistently.The equations may also include external electric or magnetic fields in the Lorentz force.Source: §2.1 · p. 4

2.2. Conservation Properties of the Vlasov-Maxwell system

  • Under the stated boundary assumptions, the continuous system conserves particle number, the distribution's L2 norm and entropy, and total particle-plus-field momentum and energy.The conservation review assumes the distribution vanishes sufficiently rapidly at velocity-space infinity and either periodic configuration space or vanishing distribution at its boundaries.Source: §2.2 · pp. 4–6

3. The Discrete Vlasov-Maxwell system

This section formulates discontinuous Galerkin discretizations for the Vlasov and Maxwell systems and examines conservation and stability in the semi-discrete scheme. It also describes the explicit SSP-RK3 update and its light-speed and velocity-space acceleration CFL constraints.

3.1. The semi-discrete Vlasov equation

  • The Vlasov discretization uses a Serendipity polynomial space and local Lax-Friedrichs flux, which can be reduced to central or upwind fluxes.Surface and volume integrals are evaluated using Gaussian quadrature of an appropriate order.Source: §3.1 · p. 7

3.2. The semi-discrete Maxwell system

  • The Maxwell system is discretized with DG interface fields obtained from central averages or upwind fluxes that solve a local face Riemann problem.The discretization omits Maxwell’s divergence constraints, which may therefore be violated during simulation.Source: §3.2 · pp. 8–9

3.3. Conservation properties of the semi-discrete system

  • The scheme conserves total particle number, while the distribution function’s L2 norm decays monotonically.The L2 result follows from the discrete scheme’s phase-space incompressibility.Source: §3.3 · pp. 9–10
  • Discrete phase-space incompressibility holds because the Vlasov-Maxwell phase-space flow is continuous across the interfaces considered, a property not guaranteed for other systems.The velocity is continuous across configuration-space faces, while the velocity-space flow is continuous across velocity-space faces.Source: §3.3 · pp. 9–10
  • If the discrete distribution remains positive, the scheme’s discrete entropy grows monotonically, reflecting regularization across phase space.This behavior follows from particle conservation and monotonic decay of the distribution function’s L2 norm.Source: §3.3 · pp. 10–11
  • Electromagnetic energy is conserved exactly with central Maxwell fluxes and bounded with upwind fluxes; upwinding causes monotonic decay of total energy.The total-energy result requires |v|^2 to belong to the approximation space.Source: §3.3 · pp. 11, 13
  • With central Maxwell fluxes and |v|^2 in the approximation space, the semi-discrete scheme conserves total particle-plus-field energy exactly.The condition on |v|^2 can be met using at least piecewise quadratic polynomials in velocity space.Source: §3.3 · p. 13
  • A projected energy remains conserved with piecewise linear basis functions, while using piecewise quadratic polynomials can represent |v|^2 directly.The conserved quantity in the linear-basis case is the projection of |v|^2.Source: §3.3 · p. 13
  • The discrete system does not conserve momentum: field discontinuities prevent global integration by parts, and increasing configuration-space resolution can strengthen momentum conservation.Momentum conservation depends only weakly on velocity-space resolution.Source: §3.3 · pp. 13–14

3.4. The time discretization of the semi-discrete Vlasov-Maxwell system

  • The fully explicit update uses SSP-RK3 to balance memory footprint and accuracy; higher-order time discretization may use single-step ADER methods such as Balsara et al..Source: §3.4 · p. 14
  • The explicit step is limited by both the speed-of-light configuration-space CFL condition and a velocity-space Lorentz-acceleration condition, which can be more restrictive.The v × B force near the velocity-grid edge can impose a particularly restrictive step limit for moderate magnetic fields.Source: §3.4 · p. 14

4. Efficient Implementation of Kinetic Update

The implementation controls high-dimensional cost through a reduced Serendipity basis, tailored quadrature, and precomputed update operations. These choices preserve conservation-relevant integration and reduce work, with structured-grid assumptions and memory trade-offs constraining their scope.

4. Efficient Implementation of Kinetic Update

  • Serendipity removes super-linear monomials above order p, reducing high-dimensional degrees of freedom while retaining tensor-space convergence on structured grids.Arbitrary unstructured-grid refinements can destroy that convergence order.Source: §4 · pp. 15–16
  • Recursive nodal layouts make each higher-dimensional element face a unisolvent expansion, reducing 5D surface-integral degrees of freedom by a factor of two to three.Only nodal information local to each face is needed for those surface integrals.Source: §4 · pp. 15–16, Figure 1 · p. 16
  • In five dimensions at p=2, exact volume integration naively requires 1024 quadrature points; under-integration risks poor particle and energy conservation and numerical instability.Source: §4 · p. 17
  • Velocity-space advection requires degree 2p+1 rather than 3p, enabling anisotropic quadrature; gains are absent at p=1 and moderate for other tested cases.Source: §4 · p. 17, Table 4 · p. 17
  • Precomputing configuration-space streaming integrals reduces update work from O(NqNp) to O(Np^2), with gains matching the Serendipity basis's quadrature-scaling reduction.Source: §4 · p. 18
  • Precomputing nonlinear force updates saves at most a factor of 2 in operations for the considered high-dimensional problems, but its additional memory was not worthwhile.Source: §4 · p. 21

5. Applications

The authors describe the benchmark set as non-exhaustive and point to separate Vlasov-Maxwell studies of the plasma sheath and current filamentation instability.

5. Applications

  • Separate studies examine the plasma sheath and current filamentation instability; the latter reports long-wavelength growth agreeing with linear theory and nonlinear saturation amplitudes comparable to two-fluid simulations.Source: §5 · p. 21

5.1. Scaling studies

Scaling studies assess solver performance using strong and weak scaling setups with the Serendipity basis.

5.1. Scaling studies

  • With a 1X3V setup, order 4, and 136 Serendipity degrees of freedom per cell, strong scaling demonstrates operation on thousands of processors; weak scaling grows from 64 × 83 cells on 16 processors to 1024 × 163 on 2048.The reported runs span hundreds of time steps; results are presented in Figure 2.Source: §5.1 · p. 21

5.2. Electrostatic applications

These electrostatic benchmarks solve Maxwell’s equations in full, though their tests resemble the Vlasov-Ampere system; the distribution functions are initialized as drifting Maxwellians. The tests examine conservation, Landau damping, and collisionless shock behavior.

5.2. Electrostatic applications

  • The electrostatic tests solve Maxwell’s equations in their entirety, while drawing an analogy to Vlasov-Ampere; their benchmark distributions are initialized as drifting Maxwellians.The initial conditions for the conservation tests use strong asymmetric flows.Source: §5.2 · p. 21, §5.2.1 · p. 22

5.2.1. Momentum and energy conservation tests

  • Energy-error convergence orders are 2.12, 2.70, and 2.75 with successive CFL halvings, while momentum errors show super-linear convergence under grid refinement.The energy orders are comparable to the expected third-order scaling of the strong-stability-preserving Runge-Kutta method; momentum orders for polynomial orders 2 and 3 are 1.46, 2.27, 3.76 and 1.99, 3.08, 6.04, respectively.Source: §5.2.1 · p. 23, Figure 3 · p. 23
  • Positivity is not guaranteed by the DG scheme, so the authors use a modified entropy-like quantity to assess entropy behavior despite small positivity violations.The proof of discrete-entropy growth assumes a positive distribution function; violations persist at high resolution because round-off errors arise where the distribution is very small.Source: §5.2.1 · p. 24
  • In the 1D conservation test, electric-field divergence errors are small relative to the largest field variation, sub-dominant to the dynamics, and decrease with resolution.The charge density used in the divergence relation comes from distribution-function moments and is not explicitly evolved in the equation system.Source: §5.2.1 · p. 24

5.2.2. Landau damping of Langmuir waves

  • Langmuir-wave damping transfers electromagnetic energy to electron thermal energy, and the reported simulations use resolutions ranging from about one cell per Debye length to 40 Debye lengths per cell.The authors report good energy conservation and decreasing divergence errors with increasing resolution.Source: §5.2.2 · pp. 26–27, Figure 6 · p. 27

5.2.3. Electrostatic shocks

  • At Mach 5, the collisionless Vlasov run forms no shock, while the fluid simulation still shocks; at Mach 1.5, density and momentum compare well but particle energy does not.The proton particle energy agrees especially poorly, consistent with protons being far from thermal equilibrium after significant acceleration.Source: §5.2.3 · pp. 27–28

5.3. Electromagnetic applications

The electromagnetic benchmarks compare kinetic Vlasov-Maxwell simulations with fluid models, showing differences in wave-particle energy exchange and turbulent current structure. The Orszag-Tang calculation also examines distribution-function structure, energy spectra, and correlation diagnostics, while noting positivity and magnetic-divergence errors.

5.3.1. Vlasov-Maxwell Riemann problem

  • In the Riemann problem, wave-particle resonances such as Landau damping shift the electron particle-energy profile upward, while total energy is conserved within 1 percent.Ripples in the electron distribution function provide evidence of phase-mixing associated with damping.Source: §5.3.1 · pp. 31–32, Figure 11 · p. 32

5.3.2. The Orszag-Tang vortex

  • Vlasov-Maxwell currents are about 20 percent stronger in developing Orszag-Tang current layers than in the 10-moment fluid model, and the difference continues to grow by 40 Ω−1.The large-scale current structures otherwise show qualitative agreement.Source: §5.3.2 · pp. 32, 35
  • Orszag-Tang spectra steepen across proton and electron inertial scales, while the collisionless Vlasov spectrum builds up at small scales because the system has no means of dissipating energy.The polynomial representation resolves beyond the grid scale.Source: §5.3.2 · p. 37, Figure 17 · p. 37
  • The k-filtered Orszag-Tang distributions retain substantial structure across scales, with signal spanning eight orders of magnitude from the full distribution to the smallest-scale filter.Small positivity violations occur in the full distribution; negative values after k-filtering are not alarming because the filter is a proxy for δf.Source: §5.3.2 · p. 34, Figure 18 · p. 38
  • Field-particle correlations provide proxies for Landau, transit-time, and cyclotron damping, and their signal is striking despite velocity-space gradients that can worsen particle-counting noise.A future publication is to explore the significance of the presented correlations.Source: §5.3.2 · pp. 34, 39
  • Small positivity violations are unphysical, and sufficiently large negative distribution values complicate analysis, although the reported fluid moments remain positive and the simulation physics appears little affected.The authors defer a systematic study of kinetic limiters; increasing velocity resolution becomes more intractable at higher dimensions.Source: §5.3.2 · pp. 39, 41
  • At the coarse resolution explored, magnetic-divergence errors are about 2 percent but appear small relative to the dynamics of interest and likely do not affect the simulated physics.The authors discuss divergence cleaning and other approaches for suppressing or avoiding such errors.Source: §5.3.2 · pp. 41–42

6. Summary and Future Work

The paper reports a conservative semi-discrete DG Vlasov-Maxwell solver with an efficient nodal implementation, validated across benchmarks through five-dimensional turbulence. It also identifies timestep, physics-scope, and velocity-space discretization limits alongside possible extensions.

6. Summary and Future Work

  • The semi-discrete DG Vlasov-Maxwell discretization conserves particles exactly and conserves energy exactly with central Maxwell fluxes, while its nodal basis implementation reduces multiplications.The particle-conservation result is independent of the Vlasov numerical flux function.Source: §6 · p. 42
  • In the five-dimensional turbulence run, the acceleration CFL is up to two times more restrictive than the speed-of-light CFL; an implicit velocity-space update is proposed.The authors propose splitting configuration- and velocity-space updates to address this timestep constraint.Source: §6 · p. 43
  • Preliminary modal-basis tests support reduced operation-count scaling, including a 3D-to-4D turnover where modal surface integrals become cheaper.The proposed structured-grid modal formulation could reduce configuration-space update scaling from O(N^2p) to O(Np).Source: §6 · p. 43
  • The implemented model is non-relativistic and collisionless; proposed extensions face aliasing and stiff nonlinear collision-term challenges, with suggested workarounds not general for highly relativistic, highly collisional plasmas.The authors identify local Taylor expansion of Lorentz factors and sufficiently low collisionality as possible routes for some applications.Source: §6 · pp. 43–44
  • Velocity-space boundaries can become problematic under particle acceleration or heating; the authors also identify possible negative distributions and accumulated Maxwell divergence errors.Proposed improvements include velocity-space grid flexibility, a positivity-preserving limiter, and a curl-conforming DG method for Maxwell’s equations.Source: §6 · p. 44
  • Greater than eight orders of magnitude of signal appears in the Orszag-Tang distribution function, within benchmarks spanning electrostatic and fully electromagnetic cases through five-dimensional turbulence.The authors describe the solver as robust and emphasize the value of its noise-free distribution function.Source: §6 · pp. 43–44
Loading 1705.05407v3…