GANDALF: A hardware-agnostic spectral solver for kinetic reduced MHD turbulence

TL;DR

GANDALF addresses the accessibility barrier to kinetic reduced MHD turbulence research by delivering a hardware-agnostic, spectral solver implemented in JAX. It combines Fourier spatial discretization with a Hermite expansion in velocity space and an integrating factor time-stepping method that exactly propagates linear Alfvén waves, enabling efficient simulations on laptops and GPUs. The paper validates GANDALF across linear Alfvén waves, nonlinear Orszag-Tang dynamics, and fully developed turbulence, showing $E(k_\perp) \propto k_\perp^{-5/3}$ in the inertial range and a Hermite spectrum $C_m^{+} \propto m^{-1/2}$ for phase mixing, while achieving machine-precision wave dispersion in several timesteps. By lowering infrastructure requirements and preserving spectral accuracy, GANDALF broadens participation in plasma turbulence research, supports rapid prototyping and education, and complements existing production codes rather than replacing them, with a roadmap toward adaptive timestepping and differentiable physics.

Abstract

We present GANDALF, a JAX-based spectral solver for Kinetic Reduced MHD (KRMHD) turbulence designed to lower infrastructure barriers to plasma turbulence research. Existing production codes require specialized HPC infrastructure and compilation expertise, limiting participation to well-resourced institutions. GANDALF addresses this barrier by leveraging JAX's hardware abstraction to run transparently on laptops, desktop GPUs, and Apple Silicon without modification, enabling single-command installation via pip. We employ Fourier spectral methods for spatial discretization and Hermite spectral basis for velocity space, combined with an exponential integrating factor method that exactly propagates linear Alfvén waves, eliminating associated numerical stiffness. Verification demonstrates research-grade accuracy: linear Alfvén waves achieve machine precision (~10^{-15} relative error), the Orszag-Tang vortex conserves energy to 10^{-6} over two Alfvén times, and driven turbulence reproduces the expected k_perp^{-5/3} cascade spectrum. GANDALF enables rapid prototyping, parameter surveys, and educational applications on commodity hardware. The code complements rather than replaces established solvers like AstroGK and Viriato, prioritizing accessibility for researchers without HPC resources. By removing infrastructure barriers while maintaining spectral accuracy, GANDALF broadens participation in fundamental plasma turbulence research, particularly benefiting students, small research groups, and institutions in developing regions.

Figures (6)

• Figure 1: Alfvén wave dispersion benchmark convergence studies. (a) Spatial convergence: relative frequency error versus grid resolution $N$ (for $N^3$ grids) with fixed timestep $\Delta t = 0.01$. The flat profile indicates temporal error dominates - spatial error has already converged at $N=32^3$. (b) Temporal convergence: relative error versus timestep $\Delta t$ with fixed resolution $N = 64^3$. Squares show measured non-zero errors ($\sim 10^{-7}$, likely from numerical roundoff); triangles mark three timesteps achieving machine precision (identically zero error), demonstrating that the exponential integrating factor treats linear Alfvén wave propagation analytically exactly. Gray dotted line shows $O(\Delta t^2)$ scaling expected from standard RK2 for comparison.
• Figure 2: Orszag-Tang vortex energy evolution. Top: Total, kinetic, and magnetic energy versus time. Selective decay causes $E_{\mathrm{mag}}/E_{\mathrm{kin}}$ to increase from 1.0 to 1.59 over 2 Alfvén crossing times, characteristic of 2D MHD inverse cascade physics. Bottom: Energy conservation error $|\Delta E/E_0|$ remains at $\sim 10^{-6}$ level despite strong nonlinear interactions, validating the energy-conserving discretization. Resolution: $N = 128^2 \times 2$, timestep $\Delta t = 0.01$.
• Figure 3: Orszag-Tang vortex 2D field structures at $t = 2\tau_A$ showing (a) vorticity $\omega = \nabla^2\Phi$, (b) current density $J_\parallel = \nabla^2\Psi$, (c) stream function $\Phi$, and (d) vector potential $\Psi$. Current sheets (intense localized $J_\parallel$) form at flow stagnation points, demonstrating nonlinear cascade to small scales. Spectral methods resolve steep gradients without oscillations. Resolution: $N = 128^2 \times 2$.
• Figure 4: Orszag-Tang convergence studies. (a) Spatial: Energy conservation error versus resolution shows exponential convergence ($\alpha = 0.076$) typical of spectral methods for nonlinear problems, reaching $10^{-4}$ levels even at modest $N$. (b) Temporal: Only smallest timestep ($\Delta t = 0.0125$) remained stable; larger timesteps developed instabilities despite passing linear CFL criterion, highlighting nonlinear stiffness beyond integrating factor's linear wave treatment.
• Figure 5: Time-averaged turbulent energy spectrum $E(k_\perp)$ from production GANDALF simulation at $N = 64^3$ resolution. Left: Kinetic energy spectrum. Right: Magnetic energy spectrum. Both components exhibit clean $k_\perp^{-5/3}$ Kolmogorov scaling (dashed reference line) across the inertial range $k_\perp \in [2, 12]$. Balanced Elsasser forcing at $k_\perp \in [1, 2]$ with $|k_z| \leq 1$ restriction provides stable energy injection respecting RMHD ordering. Time-averaged over $t \in [180, 200]\tau_A$ after initial transients decay. Total energy: $1.73 \times 10^4$ (normalized units); magnetic fraction: 0.46 (approaching equipartition). Data from checkpoint at $t = 200\tau_A$.