Gkeyll Plasma Simulation Framework

• Gkeyll is a continuum plasma simulation framework that leverages high-order discontinuous Galerkin methods across kinetic, gyrokinetic, and fluid models.
• It supports a broad model hierarchy, ranging from full kinetic Vlasov solvers to extended-fluid and general relativistic finite-volume formulations.
• It incorporates conservative discretization, exact conservation, and efficient basis strategies to accurately simulate turbulence, reconnection, and boundary interactions.

Gkeyll is a continuum plasma simulation framework built around high-order discontinuous Galerkin methods for kinetic, gyrokinetic, fluid, and extended-fluid equations. In the recent literature it appears as a full-$f$ long-wavelength gyrokinetic solver for edge and scrape-off-layer turbulence, a fully kinetic continuum Vlasov–Maxwell code for reconnection and weakly collisional plasma dissipation, a five-moment and ten-moment extended-MHD platform, a finite-volume framework for general relativistic electromagnetism and hydrodynamics, and, when wrapped externally, a high-fidelity C/CUDA plasma simulation suite usable inside differentiable multiphysics workflows (Hakim et al., 2016, Pezzi et al., 2021, TenBarge et al., 2019, Gorard et al., 2024, Coughlin et al., 17 Nov 2025).

1. Computational scope and model hierarchy

Gkeyll is described as a “general plasma simulation framework,” a “highly extensible code framework,” and a “multi-physics, discontinuous-Galerkin code” that includes solvers for full Vlasov, gyrokinetics, and fluid models (Pezzi et al., 2021, Mandell et al., 2021, Hoffmann et al., 13 Oct 2025). In different studies it is used as a full-$f$, continuum, electrostatic gyrokinetic code; a full-$f$ electromagnetic gyrokinetic code in the long-wavelength or drift-kinetic limit; a fully kinetic continuum Vlasov–Maxwell solver; a ten-moment and five-moment two-fluid extended-MHD solver; and a finite-volume code for general relativistic electromagnetism and hydrodynamics (Bernard et al., 2020, Mandell et al., 2021, TenBarge et al., 2023, TenBarge et al., 2019, Gorard et al., 2024, Hanebring et al., 15 Jan 2026).

Within gyrokinetics, Gkeyll is repeatedly used in full-$f$ mode, meaning that it evolves the total distribution function $f_s$ rather than a perturbation $\delta f_s$ around a fixed background. The cited edge and SOL studies emphasize this choice because fluctuations can be order-unity, profiles evolve on comparable timescales to turbulence, and source and sink terms are intrinsic to the dynamics (Hoffmann et al., 13 Oct 2025, Shi et al., 2017). In reconnection and dissipation studies, by contrast, the same framework is used as a continuum Vlasov–Maxwell solver retaining full phase-space information for electrons and ions in 2.5D–3V or 2D–3V settings (Pezzi et al., 2021, TenBarge et al., 2023).

The model hierarchy spans reduced and unreduced descriptions. At one end are full kinetic continuum formulations of Vlasov–Poisson, Vlasov–Maxwell, and VPFP systems; at the other are five-moment, ten-moment, gyrofluid, and general relativistic finite-volume formulations (Cagas et al., 2016, TenBarge et al., 2019, Gorard et al., 2024). This breadth is central to Gkeyll’s role in the literature: it is not restricted to a single plasma regime or closure family. A plausible implication is that Gkeyll is best understood as a numerical framework whose identity is defined less by one equation set than by a common discretization philosophy and by interoperability across model fidelities.

2. Governing equations and physical formulations

A recurring formulation in Gkeyll is the Hamiltonian kinetic form. In one generic presentation, the core evolution equation is written as

$f_t - \{H,f\} = 0,$

with the conservative form

$$(\mathcal{J} f)_t + \nabla\cdot(\mathcal{J}\boldsymbol{\alpha}f) = 0,$$

where $\mathcal{J}$ is the phase-space Jacobian and $\boldsymbol{\alpha}$ the characteristic phase-space velocity (Hakim et al., 2016). In full-$f$0 long-wavelength gyrokinetics, the distribution $f$1 evolves according to

$f$2

with Hamiltonian

$f$3

in the electrostatic case, or with explicit $f$4 dependence in the electromagnetic symplectic formulation (Shi et al., 2017, Bernard et al., 2020, Mandell et al., 2021, Hoffmann et al., 13 Oct 2025).

The corresponding long-wavelength gyrokinetic Poisson or quasineutrality equation appears in several variants. One representative form is

$f$5

with $f$6 built from reference densities and $f$7 factors, often in a Boussinesq or linearized-polarization approximation (Mandell et al., 2021, Mandell et al., 2021). In electromagnetic gyrokinetics, Gkeyll also solves parallel Ampère’s law,

$f$8

or an equivalent time-differentiated form for $f$9 designed to avoid the Ampère cancellation problem (Mandell et al., 2021).

In fully kinetic mode, Gkeyll discretizes the Boltzmann–Maxwell or Vlasov–Maxwell system. A representative species equation is

$f$0

coupled to Maxwell’s equations for $f$1 and $f$2 (Pezzi et al., 2021). In the VPFP sheath application wrapped into a differentiable workflow, Gkeyll serves specifically as a continuum kinetic solver for a Vlasov–Poisson–Fokker–Planck sheath model whose steady-state current closes a macroscopic Z-pinch circuit ODE (Coughlin et al., 17 Nov 2025).

Gkeyll’s extended-fluid branch evolves moment hierarchies rather than full distributions. In the ten-moment formulation, for each species,

$f$3

$f$4

$f$5

A closure for $f$6 is then required, and in several studies this is provided by a local isotropization or Landau-fluid-inspired model (TenBarge et al., 2019, Hanebring et al., 15 Jan 2026).

In general relativistic finite-volume applications, Gkeyll treats generic conservation laws on curved spacetime,

$f$7

and solves them by transforming states to a locally flat tetrad basis at cell faces, applying a special relativistic Riemann solver, and then performing a geometric flux correction (Gorard et al., 2024). This places the relativistic finite-volume branch in the same framework as the kinetic and extended-fluid branches, but with a distinct discretization backend.

3. Numerical architecture and conservative discretization

The most characteristic numerical ingredient of Gkeyll is its discontinuous Galerkin discretization. The 2016 edge/SOL overview describes an energy-conserving mixed DG/CG method applicable to a general class of Hamiltonian equations, with DG for the distribution function and CG for the fields or Hamiltonian (Hakim et al., 2016). Later gyrokinetic studies use DG in all phase-space dimensions with piecewise-linear or piecewise-quadratic basis functions, explicit strong-stability-preserving Runge–Kutta time stepping, and quadrature rules chosen to preserve the algebraic conservation structure (Shi et al., 2017, Francisquez et al., 2020).

A major numerical theme is exact discrete conservation. The paper on the nonlinear gyro-averaged Dougherty operator states that the operator analytically has the advective-diffusive form of Fokker–Planck operators, has a non-decreasing entropy functional, and conserves particles, momentum and energy; discretely, these conservative properties are maintained exactly as well, independent of numerical resolution (Francisquez et al., 2020). That implementation uses weak equality and recovery-based DG constructions to define primitive moments and second derivatives consistently. In the same paper, validation includes relaxation tests, collisional Landau-damping benchmarks, and 5D gyrokinetic turbulence on helical open field lines (Francisquez et al., 2020).

Positivity control is another persistent concern. In the early full-$f$8 open-field-line gyrokinetic work, negative values introduced by collisions are repaired by setting negative nodes to zero, rescaling to restore density, and then applying energy-correction drag operators in $f$9 and $f$0; the correction is designed to vanish with refinement (Shi et al., 2017). The 2016 overview likewise emphasizes positivity as a principal reason for preferring DG over standard high-order advection schemes in full-$f$1 edge applications (Hakim et al., 2016).

The framework also supports specialized basis and collision designs for efficiency. A Maxwellian-weighted DG basis was reported to achieve the same error with 4 times less computational cost in 1v, or 16 times lower cost in the 2 velocity dimensions of gyrokinetics, assuming memory bandwidth is the limiting factor (Hakim et al., 2016). In a later ASDEX-Upgrade SOL study, a conservative implicit BGK collision operator was added to the full-$f$2 gyrokinetic solver. In that collisional regime, the implicit BGK operator enabled 2D axisymmetric simulations to run 56 times faster to completion than simulations with the more expensive Lenard–Bernstein–Dougherty operator, while agreeing well with the latter in the more collisional limit (Liu et al., 30 Jul 2025).

Not all Gkeyll solvers are DG-based. The five-moment, ten-moment, and general relativistic modules use finite-volume and wave-propagation or Riemann-solver-based algorithms. The 2024 general relativity paper extends an existing curvilinear finite-volume infrastructure by adding a spacetime tetrad and a post-hoc geometric flux correction so that special relativistic Riemann solvers can be reused on curved spacetimes (Gorard et al., 2024). This coexistence of DG kinetic solvers and finite-volume fluid or relativistic solvers is a defining architectural feature rather than an inconsistency.

4. Geometry, boundaries, and wall interaction models

Gkeyll is used in geometries ranging from straight open-field-line slabs to helical simple-magnetized-torus models, Miller-parameterized tokamak edge equilibria, axisymmetric divertor/SOL domains, and black-hole spacetimes (Shi et al., 2017, Bernard et al., 2020, Hoffmann et al., 13 Oct 2025, Shukla et al., 11 Apr 2025, Gorard et al., 2024). A common strategy in magnetized plasma applications is field-aligned coordinates. In straight open-field-line turbulence studies the coordinates are Cartesian with $f$3; in helical SOL studies a non-orthogonal field-aligned coordinate system is constructed so that magnetic shear enters through the metric; in predictive tokamak edge/SOL simulations the coordinates are

$f$4

with $f$5 derived from EFIT equilibrium data (Shi et al., 2017, Mandell et al., 2021, Hoffmann et al., 13 Oct 2025).

Boundary conditions at material surfaces are especially important in Gkeyll’s edge and sheath literature. For open-field-line gyrokinetics, the framework uses conducting-sheath boundary conditions at the ends of field lines. The sheath potential at the boundary is obtained from the gyrokinetic Poisson equation, low-energy electrons are reflected, and ions and sufficiently energetic electrons are lost to the wall (Shi et al., 2017, Bernard et al., 2020, Mandell et al., 2021). In the 2017 straight-field-line turbulence paper, the sheath model explicitly allows currents to flow through the walls, in contrast to logical-sheath formulations that impose local zero net current (Shi et al., 2017). In the 2016 classical-sheath comparison, the kinetic solver uses absorbing-wall particle boundary conditions with $f$6 at both walls, while the five-moment fluid solver employs “novel boundary conditions” based on vacuum ghost states and a Riemann solver so that the sheath forms self-consistently without prescribed wall fluxes (Cagas et al., 2016).

Open and closed field lines can coexist in a single simulation. The predictive tokamak edge/SOL work uses twist-and-shift boundary conditions on closed field lines and conducting-sheath conditions on open field lines, with a limiter connection handled by enforcing $f$7 at the sheath/twist-and-shift corner for continuity (Hoffmann et al., 13 Oct 2025). Axisymmetric STEP-like divertor studies add mirror-force physics through the kinetic characteristics and show that this substantially changes parallel acceleration and electrostatic potential drops in long-leg geometries (Shukla et al., 11 Apr 2025). In general relativistic applications, boundary treatment includes excision at black-hole horizons with zero-flux interior conditions in the finite-volume solver (Gorard et al., 2024).

The literature also makes clear that boundary models are often the dominant approximation layer. Several papers explicitly note missing physics such as magnetic presheaths, oblique-angle incidence, secondary electron emission, sputtering, detailed neutral recycling, or geometrically exact limiter and divertor corners (Cagas et al., 2016, Shi et al., 2017, Hoffmann et al., 13 Oct 2025). The framework is therefore best characterized as providing a hierarchy of boundary models whose sophistication is problem-dependent.

5. Scientific applications across plasma and relativistic regimes

In magnetic reconnection and kinetic dissipation studies, Gkeyll functions as a fully kinetic continuum Vlasov–Maxwell code. In the dissipation-measures comparison, Gkeyll is the continuum counterpart to VPIC and HVM, using discontinuous Galerkin phase-space discretization, piecewise quadratic Serendipity basis functions, and a conservative DG implementation of the nonlinear Dougherty collision operator (Pezzi et al., 2021). For the weakly collisional reconnection problem studied there, the Gkeyll results display good agreement with VPIC, with overall structure in both energy-based and distribution-function-based diagnostics agreeing better for electrons than for protons (Pezzi et al., 2021). In the later Eulerian-versus-Lagrangian energization analysis, Gkeyll’s Eulerian Vlasov description is used to decompose $f$8 into parallel, diamagnetic, curvature, agyrotropic, and polarization contributions, leading to the conclusion that away from the X-point, diamagnetic and agyrotropic drifts dominate bulk electron energization, whereas direct acceleration dominates at the X-point (TenBarge et al., 2023).

In fusion edge and SOL physics, Gkeyll has been used from early proof-of-feasibility studies to predictive full-$f$9 workflows. The 2016 and 2017 open-field-line gyrokinetic papers demonstrated turbulence in LAPD-like and SOL-like geometries with sources, sheath losses, and cross-field transport, including intermittency, broadband fluctuation spectra, and sheath-regulated parallel losses (Hakim et al., 2016, Shi et al., 2017). In the Texas Helimak limiter-biasing study, Gkeyll is used as a full-$f_s$0, continuum, electrostatic gyrokinetic code with conducting-sheath boundary conditions and biased wall potentials, leading to the conclusion that turbulence is mostly driven by the interchange instability and that shear rates are mostly less than local linear growth rates; the main effect of biasing is on bulk transport and equilibrium density gradients rather than direct shear stabilization (Bernard et al., 2020). Electromagnetic open-field-line simulations later extended this program to high-$f_s$1 helical SOL models, where a 10% increase in cross-field transport near the midplane yielded broader electron heat-flux widths and a 25% reduction of the peak electron heat flux to the endplates (Mandell et al., 2021). A related magnetic-shear study found that stronger shear reduces perpendicular transport and steepens profiles, while electromagnetic effects slightly increase transport in strong-shear cases (Mandell et al., 2021).

A further step is represented by predictive global edge/SOL gyrokinetics. In the TCV study, Gkeyll uses only magnetic geometry, heating power, and particle inventory as inputs, with an adaptive sourcing algorithm that controls energy injection and emulates particle sourcing due to neutral recycling (Hoffmann et al., 13 Oct 2025). For discharge #65125, the simulations compare reasonably well with Thomson scattering and Langmuir probe data and reproduce blob transport and self-organized electric fields (Hoffmann et al., 13 Oct 2025). Applying the same framework to TCV discharges #65125 and #65130 suggests a mechanism for the improved confinement of negative triangularity: the simulations indicate that negative triangularity increases the $f_s$2 flow shear by about 20% in those cases, correlating with reduced turbulent losses and a modest redistribution of power exhaust to the vessel wall (Hoffmann et al., 13 Oct 2025).

Gkeyll is also used as a kinetic complement to fluid SOL/divertor tools. In the STEP-like axisymmetric comparison with SOLPS-ITER, kinetic Gkeyll simulations show significant upstream ion mirror trapping, a larger electrostatic potential drop, failure of the assumption of equal main-ion and impurity temperatures, and much stronger impurity confinement to the divertor region. For cold argon cases, the reduction in upstream impurity density relative to the fluid calculation can be by orders of magnitude (Shukla et al., 11 Apr 2025). In the more collisional ASDEX-Upgrade axisymmetric SOL study, the new implicit BGK operator yields macroscopic profiles and divertor heat fluxes in good agreement with converged LBD runs while dramatically reducing cost (Liu et al., 30 Jul 2025).

Outside fusion edge plasma physics, Gkeyll’s extended-fluid branch has been applied to space and astrophysical plasmas. In the MMS Burch-event study, the ten-moment model reproduces fast reconnection, pressure agyrotropy at the X-line, and electron temperature anisotropy magnitudes comparable to MMS, while the five-moment model over-amplifies lower-hybrid drift turbulence and fails to maintain a coherent diffusion region (TenBarge et al., 2019). In the 2026 collisionless-dynamo study, the 10-moment solver with $f_s$3 captures Weibel seed-field generation, transition to a collisionless turbulent dynamo, and saturation near equipartition, with the electron heat-flux closure acting as an effective regulator of magnetic Reynolds number (Hanebring et al., 15 Jan 2026).

Finally, Gkeyll’s finite-volume branch has been extended to general relativistic electromagnetism and hydrodynamics. The “tetrad-first” algorithm transforms primitive and conserved variables into a locally flat spacetime basis at each face, applies a special relativistic Riemann solver, and then applies a geometric correction. In the reported tests, this improves convergence and stability relative to direct curved-space Riemann solvers, especially near rapidly spinning black holes (Gorard et al., 2024).

6. Differentiable workflows, computational cost, and open limitations

A recent development situates Gkeyll inside differentiable multiphysics pipelines. In the Z-pinch case study, Gkeyll is explicitly described as a non-differentiable C/CUDA plasma simulation suite used as the high-fidelity plasma kinetic engine in a JAX-based workflow (Coughlin et al., 17 Nov 2025). The Tesseract abstraction layer and tesseract_jax adapter treat Gkeyll, neural surrogates, and analytic approximations as interchangeable implementations of the same $f_s$4 map, while custom finite-difference VJP/JVP rules provide gradients for Newton solves and outer optimization loops (Coughlin et al., 17 Nov 2025). The paper’s practical implication is direct: Gkeyll does not need to be rewritten in JAX in order to be embedded in an end-to-end differentiable multiphysics stack (Coughlin et al., 17 Nov 2025).

Computational cost remains a defining constraint. In the differentiable Z-pinch study, a single Gkeyll run for one circuit-ODE trajectory with kinetic closure takes about 18 hours on a single NVIDIA A100 GPU (Coughlin et al., 17 Nov 2025). In predictive tokamak edge/SOL simulations, milliseconds of medium-sized turbulence with kinetic electrons and ions are reported to be feasible in a few hundred GPU hours, with the baseline TCV case costing about 750 GPU-hours per millisecond on 8 A100 GPUs and the fine case about 2000 GPU-hours per millisecond on 16 GPUs (Hoffmann et al., 13 Oct 2025). Axisymmetric reactor-SOL kinetic runs are expensive but tractable, typically taking 1–3 days on 8 multi-GPU nodes in the STEP-like impurity study (Shukla et al., 11 Apr 2025). These numbers indicate that the framework is already being used in parameter studies, but they also delimit the accessible dimensionality, collisionality, and resolution.

The literature repeatedly identifies limitations rather than presenting a single universal model. Common approximations include electrostatic or long-wavelength gyrokinetics, Boussinesq or linearized polarization, simplified source models, static or reduced neutral treatments, model collision operators in place of full Landau dynamics, reduced mass ratios in some reconnection and dynamo studies, and simplified geometry relative to true diverted tokamak equilibria or fully three-dimensional wall configurations (Mandell et al., 2021, Hoffmann et al., 13 Oct 2025, Shukla et al., 11 Apr 2025, TenBarge et al., 2019). In differentiable use cases, finite-difference gradient estimates are approximate and can be noisy, especially when coupled to expensive stiff kinetic solves (Coughlin et al., 17 Nov 2025). In continuum-versus-PIC comparisons, Gkeyll’s noise-free phase-space representation is advantageous for entropy and agyrotropy diagnostics, but spatial resolution can be coarser and some diagnostics can appear broader as a result (Pezzi et al., 2021).

Future directions stated across the cited work are correspondingly diverse: more realistic geometry, neutral and impurity models, improved collision operators and implicit schemes, stronger electromagnetic and FLR fidelity, data-driven or higher-fidelity surrogate closures, adjoint or tangent-linear differentiability support beyond finite differences, and extensions to dynamic spacetimes or relativistic multi-fluid systems (Liu et al., 30 Jul 2025, Hoffmann et al., 13 Oct 2025, Coughlin et al., 17 Nov 2025, Gorard et al., 2024, Hanebring et al., 15 Jan 2026). Taken together, these studies suggest that Gkeyll is less a fixed code for a single plasma model than a numerically conservative continuum framework whose research program is to push high-fidelity kinetic and extended-fluid simulation into regimes where turbulence, boundaries, geometry, and multiscale coupling all matter simultaneously.