Source-linked AI summary
High-Order Structure-Preserving SBP Finite Difference Methods for the Vlasov-Maxwell System on Matrix-Free GPUs
Robin Dymér, Ken Mattsson, Murtazo Nazarov
TL;DR
Fine-scale structures challenge stable Vlasov–Maxwell simulation with central SBP advection operators. The paper develops an upwind SBP method with explicit Runge–Kutta integration and a matrix-free CUDA implementation. It reports high-order accuracy, reproduces standard instabilities, and establishes fully discrete mass conservation with momentum preserved up to truncation error.
Problem
Central SBP advection operators do not damp high-frequency modes in advection-dominated problems such as the Vlasov problem.
Source: §3 · p. 7Method
The method combines high-order upwind SBP finite differences, fourth-order five-stage SSP explicit Runge–Kutta integration, and a matrix-free CUDA implementation.
Source: §1 · p. 2, §7 · p. 23Results
The authors verify high-order accuracy and reproduce diocotron, Weibel, and Kelvin–Helmholtz instabilities; the largest reported run uses 420 million degrees of freedom on one NVIDIA L40.
Source: §1 · p. 2, §7 · p. 23Takeaways & Limitations
At the fully discrete level, mass is conserved exactly and momentum is preserved up to truncation error.
Source: §7 · p. 23Upwind stabilization dissipates total energy, and a uniform L2 decay estimate is not established because dissipation may vanish for constant or nearly constant solutions.
Source: §3.1 · p. 12, §4.1 · p. 14
Abstract
from arXiv · showhide
In this paper, we present a high-order, stable summation-by-parts (SBP) finite difference method for solving the Vlasov-Maxwell system in a 2D2V phase space. Central SBP operators for the advection terms are not stable when the solution becomes non-smooth and fine-scale filamentary structures develop, as is typical in high-dimensional Vlasov-Maxwell simulations. To address this issue, the method is stabilized using high-order upwind SBP operators. High-order explicit Runge-Kutta methods are employed for time integration. We prove that the fully discrete scheme exactly conserves mass and preserves momentum up to truncation error. Furthermore, we present a matrix-free implementation of the method on modern GPU architectures. A range of challenging benchmark problems is solved to demonstrate the accuracy, robustness, and performance of the proposed scheme.
1. Introduction.
The introduction motivates kinetic plasma simulation by the limits of fluid models and the computational challenges of resolving fine-scale Vlasov–Maxwell dynamics. It presents an upwind SBP finite difference method with conservation analysis, a matrix-free implementation, and numerical tests.
1. Introduction.
- Fluid models such as MHD cannot fully capture microscopic particle interactions, including non-equilibrium dynamics, wave-particle interactions, and small-scale turbulence.The Vlasov–Maxwell system provides a kinetic description of collisionless plasmas under self-consistent electromagnetic fields.Source: §1 · p. 1
- Fine-scale phase-space structures make accurate simulation challenging, and failure to capture them may cause numerical instabilities and nonphysical solutions.The Vlasov equation is scalar but is solved over spatial and velocity variables, making the full physical phase space six-dimensional.Source: §1 · p. 1
- Finite difference methods offer efficiency and high-order accuracy, but their application to the fully coupled Vlasov–Maxwell system remains limited.Prior approaches include particle-in-cell and semi-Lagrangian methods, each with stated accuracy, noise, smoothness, or conservation limitations.Source: §1 · p. 2
- The proposed high-order SBP discretization uses tensor-product operators and Lax–Friedrichs-type upwind stabilization; its analysis establishes conservation properties while noting stabilization-induced energy dissipation.The paper also analyzes Forward Euler and discusses extension to explicit Runge–Kutta methods, showing preservation of linear invariants at the fully discrete level.Source: §1 · p. 2
- The matrix-free implementation avoids assembling global operators, and numerical experiments report expected-order accuracy and robustness across convergence studies and standard plasma instabilities.The authors emphasize that the paper's focus is numerical analysis and structure-preserving properties.Source: §1 · p. 2
2. Preliminaries.
The preliminaries specify the single-species 2D2V Vlasov–Maxwell setting and its physical conservation properties. These properties motivate a discretization that carries continuous integration-by-parts structure into the discrete setting.
2.1. Governing equations.
- The analysis adopts nondimensional units, one species with charge-to-mass ratio one, and periodic boundary conditions on bounded phase-space domains.The 2D2V reduction retains the out-of-plane magnetic component Bz and writes the unknown vector as U=(f, Ex, Bz, Ey)^T.Source: §2.1 · p. 3, §2.2 · pp. 3–4
2.2. Vlasov–Maxwell equations in 2D2V.
- In 2D2V, the Vlasov–Maxwell system has continuous conservation properties for mass, momentum, total energy, and the L2 norm of f.The momentum expression includes particle and electromagnetic momentum.Source: §2.2 · pp. 4–5
- The paper targets a numerical discretization that carries these conservation properties forward because SBP operators mimic integration by parts at the discrete level.Source: §2.2 · p. 6
3. Summation-by-parts discretization.
The discretization extends one-dimensional SBP operators to 2D2V phase space and stabilizes central advection operators with upwind dissipation. Its semi-discrete analysis establishes conservation and stability properties while identifying total-energy dissipation.
3. Summation-by-parts discretization.
- Central SBP differences do not damp high-frequency modes in advection-dominated problems, so upwind SBP operators add controlled dissipation while preserving the SBP stability framework.The operators are introduced to suppress spurious modes.Source: §3 · p. 7
- The method extends one-dimensional difference operators to multidimensional phase space using Kronecker products on rectangular spatial and velocity grids.Directional derivatives are defined across the x, y, v, and w coordinates.Source: §3 · pp. 7–8
3.1. Semi-discrete approximation of Vlasov–Maxwell.
- Lax–Friedrichs flux splitting separates advection coefficient matrices into nonnegative- and nonpositive-eigenvalue parts before applying upwind SBP operators.The splitting follows and includes a parameter that can increase numerical dissipation.Source: §3.1 · pp. 8–9
- The semi-discrete upwind SBP scheme conserves mass and momentum, satisfies a discrete L2 stability estimate, propagates an initially satisfied discrete Gauss law, and dissipates total energy.The Gauss-law propagation result is stated for the semi-discrete level.Source: §3.1 · pp. 9, 11, 12
4. Fully discrete approximation.
The fully discrete scheme uses Forward Euler for its conservation analysis and SSP-RK(5,4) for numerical validation. It preserves mass and momentum up to truncation error, while upwind stabilization dissipates total energy.
4. Fully discrete approximation.
- The numerical validations use SSP-RK(5,4), whose convex-combination structure extends the Forward Euler conservation properties to explicit Runge–Kutta time integration.SSP-RK(5,4) is a fourth-order, five-stage method.Source: §4 · p. 12, §4.1 · p. 13
4.1. Forward Euler time discretization.
- Forward Euler exactly preserves discrete mass and preserves discrete momentum up to truncation error.The conservation statement is given for the fully discrete upwind SBP scheme.Source: §4.1 · pp. 13–14
- Although the L2 norm is numerically stable, a rigorous proof is challenging because upwind dissipation may vanish for constant or nearly constant solutions.Source: §4.1 · p. 14
- Upwind stabilization causes total energy to decrease rather than remain conserved.Source: §4.1 · p. 14
- An initially satisfied discrete Gauss law remains satisfied for all time at the semi-discrete level.Source: §4.1 · p. 14
5. GPU implementation.
The GPU implementation applies SBP operators through local stencils without assembling differentiation matrices. It stores the distribution and fields on their respective grids and advances them with separate GPU kernels.
5. GPU implementation.
- The distribution function is flattened onto a one-dimensional array, while electromagnetic fields are stored on the two-dimensional spatial grid, reducing field-variable memory use.Source: §5 · p. 15
- Stencil coefficients are stored in CUDA constant memory, and periodic stencil evaluations implement derivatives in the spatial and velocity directions.Source: §5 · p. 15
- GPU threads evaluate the Vlasov right-hand side at phase-space points, while atomic additions accumulate velocity-grid contributions to spatial current values.Source: §5 · p. 15
- Each Runge–Kutta stage recomputes current and evaluates the Maxwell right-hand side; separate kernels update the distribution and fields.The Maxwell update uses a separate kernel on the two-dimensional spatial grid.Source: §5 · p. 16
- The implementation is fully matrix-free, avoiding explicit differentiation-matrix storage that would be prohibitively expensive in high-dimensional simulations.Tensor-product SBP operators are applied through local stencil operations.Source: §5 · pp. 15–16
6. Numerical experiments.
The experiments verify convergence and examine the solver on diocotron, Weibel-type, and Kelvin–Helmholtz instabilities. They report field- and parameter-dependent behavior, Gauss’ law residuals, and comparisons with prior simulations.
6. Numerical experiments.
- Unless stated otherwise, all simulations use periodic finite-difference operators of order 6, ks = 1 for stabilization, and single precision on one NVIDIA L40 GPU.Source: §6 · p. 16
6.1. Verification of solver.
- Figure 1 and Table 1 show the expected convergence rates for all variables under grid refinement.Table 1 reports L2 errors and observed rates for f, Ex, Ey, and Bz with second-, fourth-, and sixth-order operators.Source: §6.1 · p. 16, Table 1 · p. 18
6.2. Diocotron instability.
- At ε = 1, weak confinement yields no diocotron vortices, while ε = 0.1 confines the plasma and forms six vortices; runtimes rise from 17 minutes to about 17 hours.Both cases use CFL = 0.5; the smaller ε case is stiffer and uses a much smaller time step.Source: §6.2 · p. 18
- The Gauss’ law residual remains small, is largest where the distribution varies rapidly, and decays slightly over time due to stabilization.The residual’s discretization error stems from solving Poisson’s equation with SBP finite differences.Source: §6.2 · pp. 18–19, Figure 4 · p. 21
6.3. Weibel-type instability.
- The Weibel-type benchmark evolves to T = 50 on a 128^4-point phase-space grid and takes about 12 minutes; its distribution is shown in two-dimensional projections and a three-dimensional view.The reported solutions differ slightly from Liu et al. because ions are simulated, Nm is unspecified in that reference, and the field initialization is stochastic.Source: §6.3 · p. 21
6.4. Kelvin–Helmholtz instability.
- The Kelvin–Helmholtz instability grows from kmax and develops similarly to Umeda et al., though the different initialization produces a slightly different growth rate.For computational efficiency, the setup uses mi/me = 16; the simulation takes about 10.8 hours.Source: §6.4 · pp. 22–23
7. Conclusions.
The paper develops a matrix-free, high-order upwind SBP discretization for the 2D2V Vlasov–Maxwell system and demonstrates it on plasma instabilities. Its fully discrete scheme conserves mass exactly and momentum up to truncation error, while upwind stabilization dissipates total energy.
7. Conclusions.
- A tensor-product, high-order upwind SBP discretization with Lax–Friedrichs flux splitting is advanced by a fourth-order, five-stage SSP Runge–Kutta method in a matrix-free CUDA implementation.The implementation avoids assembling differentiation operators and applies them through local stencil evaluations, yielding a low memory footprint.Source: §7 · p. 23
- The method verifies high-order accuracy and reproduces the diocotron, Weibel, and Kelvin–Helmholtz instabilities; its largest run used 420 million degrees of freedom on a single NVIDIA L40.Total energy is not conserved because upwind stabilization dissipates it.Source: §7 · p. 23
- The Kelvin–Helmholtz figure shows f(x, y, t) at mx = 128, my = 160, and mv = mw = 80, with position x and y on the axes.Source: Figure 7 · p. 24
- Extending the framework to 3D3V is identified as the next step, requiring multiple GPUs and likely an implicit treatment of the stiff magnetized regime.Source: §7 · p. 24