FIRM3D: GPU-Accelerated 3D Particle Dynamics

• FIRM3D is a GPU-accelerated open-source suite for modeling collisionless energetic particle dynamics in 3D magnetic confinement.
• It integrates guiding-center equations using both symplectic and Runge–Kutta methods, providing diagnostics such as Poincaré maps and orbit classification.
• The framework features minimal dependencies, interoperability with stellarator optimization tools, and scalable CPU/GPU parallelization for advanced fusion research.

Searching arXiv for the primary FIRM3D paper and directly related software/method papers. FIRM3D is an open-source Python/C++/CUDA software suite for modeling energetic particle dynamics in 3D magnetic fields, with a specific focus on fast-ion reduced modeling in fully three-dimensional magnetic confinement geometries. It targets the collisionless energetic particle regime relevant to alpha particles and fast ions born by fusion reactions or auxiliary heating, and it operates at the interface between stellarator equilibrium modeling and MHD wave stability/structure computations. In practice, it enables Monte Carlo sampling and rapid guiding-center integration of energetic particles in static equilibrium fields superposed with time-dependent MHD waves, while providing diagnostics such as Poincaré maps, orbit classification, and weighted Birkhoff averaging. The framework is packaged as a standalone, lightweight, GPU-accelerated code with minimal dependencies and tight interoperability with BOOZ_XFORM, AE3D, FAR3D, and SIMSOPT (Paul et al., 16 May 2026).

1. Definition, purpose, and scope

FIRM3D was developed to efficiently predict energetic-particle orbit confinement, loss mechanisms, and transport in optimized 3D equilibria, including quasisymmetric stellarators, while also incorporating the effects of Alfvénic activity that can degrade energetic-particle confinement. Its motivating use case is the regime in which energetic particles are largely collisionless over the timescales of interest, so that orbit integration in prescribed electromagnetic fields becomes a tractable reduced model for confinement and transport analysis. This need persists even in equilibria engineered for energetic-particle confinement via quasisymmetry, because Alfvén eigenmodes excited by energetic-particle free energy can resonantly transport those same particles (Paul et al., 16 May 2026).

The code’s novelty lies in its emphasis on a lightweight, GPU-accelerated, open-source framework with minimal dependencies, strong integration with the stellarator optimization ecosystem, and support for both symplectic and high-order Runge–Kutta integrators. Relative to the guiding-center functionality inherited from SIMSOPT, FIRM3D adds energetic-particle-specific physics and diagnostics not typically required in stellarator optimization, notably orbit classification, Poincaré maps, and weighted Birkhoff averaging. Its intended users include plasma physicists and fusion engineers studying energetic-particle-driven wave stability and transport in stellarators and advanced tokamaks, as well as private fusion companies pursuing stellarator concepts (Paul et al., 16 May 2026).

2. Physical model and coordinate representation

FIRM3D integrates the standard nonrelativistic guiding-center equations in prescribed fields. The state is evolved in reduced phase space $(X, v_\parallel, \mu)$ under the guiding-center approximation, which assumes gyroradius much smaller than equilibrium scale lengths and treats the magnetic moment as approximately invariant when fields vary slowly in time. Coulomb collisions are neglected, consistent with the intended collisionless energetic-particle regime. The model therefore describes reduced energetic-particle dynamics rather than full-orbit Lorentz dynamics or self-consistent wave-particle evolution (Paul et al., 16 May 2026).

The formulation uses the standard definitions

$$\mathbf{b} = \frac{\mathbf{B}}{|\mathbf{B}|}, \qquad \mu = \frac{m v_\perp^2}{2B},$$

with Hamiltonian

$H = \frac{1}{2} m v_\parallel^2 + \mu B + q\Phi,$

where $\Phi$ is the electrostatic potential of the perturbation. The effective fields are

$$\mathbf{B}^* = \mathbf{B} + \frac{mc}{q} v_\parallel \nabla \times \mathbf{b}, \qquad B_\parallel^* = \mathbf{b}\cdot \mathbf{B}^*,$$

and

$$\mathbf{E} = -\nabla \Phi - \frac{1}{c}\frac{\partial \mathbf{A}}{\partial t}, \qquad E_\parallel = \mathbf{b}\cdot \mathbf{E}.$$

The parallel-velocity evolution is written as

$$\dot{v}_\parallel = -\frac{\mu}{m}\,\mathbf{b}\cdot \nabla B - \frac{q}{m} E_\parallel.$$

In time-independent fields, the Hamiltonian is conserved; in perfectly quasisymmetric fields, the corresponding canonical momentum is conserved (Paul et al., 16 May 2026).

A central design choice is the use of Boozer coordinates $(\psi,\theta,\zeta)$ for equilibria and wave fields. In these coordinates, the static MHD equilibrium field is supplied through BOOZ_XFORM and represented as a Fourier series on a uniform radial grid,

$$B(\psi,\theta,\zeta) = \sum_{m,n} B_{mn}(\psi)\cos(m\theta - n\zeta + \phi_{mn}),$$

with radial evaluation by Lagrange interpolation. Time-dependent perturbations from AE3D or FAR3D are represented in the same coordinate system and superposed on the equilibrium field at specified mode frequency $\omega$. This coordinate representation is central to the code’s interoperability with stellarator equilibrium and wave-stability workflows (Paul et al., 16 May 2026).

Energetic-particle distributions are introduced by Monte Carlo sampling from user-specified distribution functions, including fusion birth spectra. Helper functions allow sampling that preserves invariants such as energy $$\mathbf{b} = \frac{\mathbf{B}}{|\mathbf{B}|}, \qquad \mu = \frac{m v_\perp^2}{2B},$$0 or canonical toroidal momentum $$\mathbf{b} = \frac{\mathbf{B}}{|\mathbf{B}|}, \qquad \mu = \frac{m v_\perp^2}{2B},$$1 when desired. This suggests a deliberate emphasis on phase-space-consistent initialization for transport and resonance studies (Paul et al., 16 May 2026).

3. Software architecture and numerical integration

The core numerics are implemented in C++ for performance, exposed through a Python front end via pybind11, and accelerated on GPUs through CUDA kernels. The code is designed to remain standalone with minimal dependencies while maintaining direct interfaces to BOOZ_XFORM for equilibria and AE3D/FAR3D for wave fields. Equilibrium and mode data are supplied as Fourier harmonics in Boozer coordinates on a uniform radial grid, and FIRM3D performs Lagrange interpolation in radius while evaluating angular harmonics on the fly (Paul et al., 16 May 2026).

Two main integration strategies are provided. The first is a symplectic explicit–implicit Euler scheme in noncanonical form, chosen for long-time conservation of invariants and exact phase-space volume preservation. The second is Dormand–Prince 5 Runge–Kutta integration, available both through Boost.odeint with adaptive step control and dense output and through a custom implementation with minimum step control to avoid excessively small steps in stiff regions. The practical guidance embedded in the framework is explicit: symplectic integration is preferred for long-time transport, chaos studies, and conservation, whereas Dormand–Prince is preferred when adaptive control, dense output, or event timing is required (Paul et al., 16 May 2026).

Parallelization is implemented on both CPUs and GPUs. On CPUs, MPI is used over Fourier harmonics for parallel field setup and OpenMP over interpolation nodes, while particle-level Monte Carlo parallelism is trivial. On GPUs, CUDA kernels implement Boozer harmonic evaluation, Lagrange interpolation, and guiding-center integration. The computational cost scales as $$\mathbf{b} = \frac{\mathbf{B}}{|\mathbf{B}|}, \qquad \mu = \frac{m v_\perp^2}{2B},$$2 for $$\mathbf{b} = \frac{\mathbf{B}}{|\mathbf{B}|}, \qquad \mu = \frac{m v_\perp^2}{2B},$$3 particles and $$\mathbf{b} = \frac{\mathbf{B}}{|\mathbf{B}|}, \qquad \mu = \frac{m v_\perp^2}{2B},$$4 time steps, with field-evaluation cost proportional to the number of Fourier harmonics and interpolation nodes. For sufficiently large ensembles, the GPU path becomes the dominant execution mode (Paul et al., 16 May 2026).

The CUDA implementation described in the related CATAPULT work is distributed within the firm3d Python package. CATAPULT implements Dormand–Prince 5 on local tricubic grids, supports equilibrium magnetic fields and optional Shear Alfvén waves, and is designed for high-throughput Monte Carlo confinement studies. That work places the GPU execution model on a more explicit footing, including local tricubic interpolation, double-precision arithmetic, and performance-oriented CUDA memory layout (Czekanski et al., 8 Apr 2026).

4. Diagnostics, resonance analysis, and transport attribution

FIRM3D includes an unusually rich set of diagnostics for a reduced energetic-particle code. Kinetic Poincaré sections in Boozer angle space are used to identify resonance islands, chaotic layers, and banana-drift diffusion, and to classify orbits. Orbit classification itself is topology-based, distinguishing banana, ripple, barely trapped, and passing trajectories, while also detecting class transitions. These diagnostics are then combined with loss detection to separate prompt losses from delayed losses and to distinguish convective from diffusive tendencies (Paul et al., 16 May 2026).

A key diagnostic for adiabaticity and integrability is the parallel invariant

$$\mathbf{b} = \frac{\mathbf{B}}{|\mathbf{B}|}, \qquad \mu = \frac{m v_\perp^2}{2B},$$5

together with its normalized variation

$$\mathbf{b} = \frac{\mathbf{B}}{|\mathbf{B}|}, \qquad \mu = \frac{m v_\perp^2}{2B},$$6

In the code’s transport interpretation, this normalized variation functions as an integrability metric. The paper associates banana losses in a quasiaxisymmetric example with exceeding a 1% nonintegrability threshold in this metric, which suggests that loss attribution is treated as a dynamical-systems problem rather than purely a confinement statistic (Paul et al., 16 May 2026).

For chaos detection, FIRM3D implements weighted Birkhoff averaging. For an observable $$\mathbf{b} = \frac{\mathbf{B}}{|\mathbf{B}|}, \qquad \mu = \frac{m v_\perp^2}{2B},$$7 along a trajectory $$\mathbf{b} = \frac{\mathbf{B}}{|\mathbf{B}|}, \qquad \mu = \frac{m v_\perp^2}{2B},$$8 over $$\mathbf{b} = \frac{\mathbf{B}}{|\mathbf{B}|}, \qquad \mu = \frac{m v_\perp^2}{2B},$$9 samples,

$H = \frac{1}{2} m v_\parallel^2 + \mu B + q\Phi,$0

with a smooth bump weight

$H = \frac{1}{2} m v_\parallel^2 + \mu B + q\Phi,$1

This suppresses endpoint effects and yields rapid, often spectral, convergence for quasiperiodic trajectories, while chaotic trajectories exhibit poor convergence. In the quasiaxisymmetric example, digit accuracy of the weighted Birkhoff average applied to $H = \frac{1}{2} m v_\parallel^2 + \mu B + q\Phi,$2 is used as a practical chaos indicator, with fewer than three digits signaling chaos (Paul et al., 16 May 2026).

The framework does not enforce a single built-in resonance formula, but its diagnostics support classical resonance analysis. The paper explicitly notes the relevance of conditions such as

$H = \frac{1}{2} m v_\parallel^2 + \mu B + q\Phi,$3

for Alfvénic Landau or cyclotron resonances, as well as bounce and drift resonances in 3D geometry computed from orbit frequencies. An example use case is a single-harmonic Alfvén eigenmode perturbation with $H = \frac{1}{2} m v_\parallel^2 + \mu B + q\Phi,$4, $H = \frac{1}{2} m v_\parallel^2 + \mu B + q\Phi,$5, and $H = \frac{1}{2} m v_\parallel^2 + \mu B + q\Phi,$6, used to reveal chaos through weighted Birkhoff averaging of canonical momenta (Paul et al., 16 May 2026).

5. Validation, benchmarking, and example studies

Validation in FIRM3D is organized around conservation properties, cross-code comparison, and parallel scaling. In time-independent fields, the code conserves energy. The symplectic integrator exhibits a stable moving time-average of energy over $H = \frac{1}{2} m v_\parallel^2 + \mu B + q\Phi,$7, whereas Dormand–Prince shows net energy drift. In a perfectly quasisymmetric field obtained by removing non-QS harmonics, toroidal canonical momentum $H = \frac{1}{2} m v_\parallel^2 + \mu B + q\Phi,$8 is conserved; with Dormand–Prince tolerance $H = \frac{1}{2} m v_\parallel^2 + \mu B + q\Phi,$9 and 64 radial Lagrange nodes, the relative error converges to approximately $\Phi$0 (Paul et al., 16 May 2026).

Cross-code benchmarking is reported against the symplectic SIMPLE code. In the precise quasisymmetric-helical equilibrium, a trapped $\Phi$1 alpha shows a relative difference in the $\Phi$2 coordinate of $\Phi$3 at $\Phi$4. For 5000 alphas sampled from a fusion birth distribution in the precise quasiaxisymmetric equilibrium and traced to $\Phi$5, FIRM3D and SIMPLE report identical loss fractions. GPU scaling is measured on NERSC Perlmutter: GPU launch latency favors CPUs below roughly $\Phi$6 particles, whereas above roughly $\Phi$7 particles the GPU is about an order of magnitude faster than 128 CPU threads (Paul et al., 16 May 2026).

Benchmark	Reported result	Context
Energy conservation	Stable moving time-average over $\Phi$8	Symplectic integrator
Canonical momentum conservation	Relative error $\Phi$9	QS field, tolerance $$\mathbf{B}^* = \mathbf{B} + \frac{mc}{q} v_\parallel \nabla \times \mathbf{b}, \qquad B_\parallel^* = \mathbf{b}\cdot \mathbf{B}^*,$$0, 64 Lagrange nodes
SIMPLE comparison	$$\mathbf{B}^* = \mathbf{B} + \frac{mc}{q} v_\parallel \nabla \times \mathbf{b}, \qquad B_\parallel^* = \mathbf{b}\cdot \mathbf{B}^*,$$1 relative difference in $$\mathbf{B}^* = \mathbf{B} + \frac{mc}{q} v_\parallel \nabla \times \mathbf{b}, \qquad B_\parallel^* = \mathbf{b}\cdot \mathbf{B}^*,$$2 at $$\mathbf{B}^* = \mathbf{B} + \frac{mc}{q} v_\parallel \nabla \times \mathbf{b}, \qquad B_\parallel^* = \mathbf{b}\cdot \mathbf{B}^*,$$3	Trapped $$\mathbf{B}^* = \mathbf{B} + \frac{mc}{q} v_\parallel \nabla \times \mathbf{b}, \qquad B_\parallel^* = \mathbf{b}\cdot \mathbf{B}^*,$$4 alpha in precise QH
Ensemble comparison	Identical loss fractions	5000 alphas in precise QA traced to $$\mathbf{B}^* = \mathbf{B} + \frac{mc}{q} v_\parallel \nabla \times \mathbf{b}, \qquad B_\parallel^* = \mathbf{b}\cdot \mathbf{B}^*,$$5
GPU scaling	$$\mathbf{B}^* = \mathbf{B} + \frac{mc}{q} v_\parallel \nabla \times \mathbf{b}, \qquad B_\parallel^* = \mathbf{b}\cdot \mathbf{B}^*,$$6 faster than 128 CPU threads	Perlmutter, above $$\mathbf{B}^* = \mathbf{B} + \frac{mc}{q} v_\parallel \nabla \times \mathbf{b}, \qquad B_\parallel^* = \mathbf{b}\cdot \mathbf{B}^*,$$7 particles

An application highlighted in the paper traces $$\mathbf{B}^* = \mathbf{B} + \frac{mc}{q} v_\parallel \nabla \times \mathbf{b}, \qquad B_\parallel^* = \mathbf{b}\cdot \mathbf{B}^*,$$8 alphas for $$\mathbf{B}^* = \mathbf{B} + \frac{mc}{q} v_\parallel \nabla \times \mathbf{b}, \qquad B_\parallel^* = \mathbf{b}\cdot \mathbf{B}^*,$$9 in a quasiaxisymmetric equilibrium at $$\mathbf{E} = -\nabla \Phi - \frac{1}{c}\frac{\partial \mathbf{A}}{\partial t}, \qquad E_\parallel = \mathbf{b}\cdot \mathbf{E}.$$0. The reported loss is $$\mathbf{E} = -\nabla \Phi - \frac{1}{c}\frac{\partial \mathbf{A}}{\partial t}, \qquad E_\parallel = \mathbf{b}\cdot \mathbf{E}.$$1. Among lost particles, banana-class accounts for $$\mathbf{E} = -\nabla \Phi - \frac{1}{c}\frac{\partial \mathbf{A}}{\partial t}, \qquad E_\parallel = \mathbf{b}\cdot \mathbf{E}.$$2, ripple $$\mathbf{E} = -\nabla \Phi - \frac{1}{c}\frac{\partial \mathbf{A}}{\partial t}, \qquad E_\parallel = \mathbf{b}\cdot \mathbf{E}.$$3, barely trapped $$\mathbf{E} = -\nabla \Phi - \frac{1}{c}\frac{\partial \mathbf{A}}{\partial t}, \qquad E_\parallel = \mathbf{b}\cdot \mathbf{E}.$$4, transition classes $$\mathbf{E} = -\nabla \Phi - \frac{1}{c}\frac{\partial \mathbf{A}}{\partial t}, \qquad E_\parallel = \mathbf{b}\cdot \mathbf{E}.$$5, and $$\mathbf{E} = -\nabla \Phi - \frac{1}{c}\frac{\partial \mathbf{A}}{\partial t}, \qquad E_\parallel = \mathbf{b}\cdot \mathbf{E}.$$6 prompt exit. The paper further states that banana losses correlate with exceeding the 1% nonintegrability threshold via $$\mathbf{E} = -\nabla \Phi - \frac{1}{c}\frac{\partial \mathbf{A}}{\partial t}, \qquad E_\parallel = \mathbf{b}\cdot \mathbf{E}.$$7 metrics and with low $$\mathbf{E} = -\nabla \Phi - \frac{1}{c}\frac{\partial \mathbf{A}}{\partial t}, \qquad E_\parallel = \mathbf{b}\cdot \mathbf{E}.$$8, implicating orbital chaos rather than convection (Paul et al., 16 May 2026).

6. Workflow, ecosystem position, and limitations

A typical FIRM3D workflow begins with a VMEC equilibrium transformed to Boozer harmonics via BOOZ_XFORM. The resulting equilibrium is loaded as a Boozer magnetic field and paired with a choice of radial Lagrange interpolation nodes. If wave effects are required, an Alfvén eigenstructure is computed with AE3D or FAR3D, exported as Boozer harmonics at frequency $$\mathbf{E} = -\nabla \Phi - \frac{1}{c}\frac{\partial \mathbf{A}}{\partial t}, \qquad E_\parallel = \mathbf{b}\cdot \mathbf{E}.$$9, and superposed in FIRM3D with explicit time dependence. An energetic-particle ensemble is then sampled from a chosen distribution, after which either the symplectic or Dormand–Prince integrator is selected according to the study objective. CPU execution is recommended for modest ensemble sizes below $$\dot{v}_\parallel = -\frac{\mu}{m}\,\mathbf{b}\cdot \nabla B - \frac{q}{m} E_\parallel.$$0, while CUDA execution is recommended for large ensembles at or above $$\dot{v}_\parallel = -\frac{\mu}{m}\,\mathbf{b}\cdot \nabla B - \frac{q}{m} E_\parallel.$$1. Diagnostics such as Poincaré sections, orbit classification, $$\dot{v}_\parallel = -\frac{\mu}{m}\,\mathbf{b}\cdot \nabla B - \frac{q}{m} E_\parallel.$$2 variation, weighted Birkhoff averages of $$\dot{v}_\parallel = -\frac{\mu}{m}\,\mathbf{b}\cdot \nabla B - \frac{q}{m} E_\parallel.$$3, and loss detection are then configured to identify resonances, chaotic layers, and loss channels (Paul et al., 16 May 2026).

Within the broader ecosystem, FIRM3D is positioned as a focused energetic-particle code rather than a general plasma simulation platform. BOOZ_XFORM supplies Boozer-harmonic equilibria; AE3D and FAR3D supply perturbation eigenstructures and frequencies; SIMSOPT provides the numerical heritage for guiding-center routines; and ASCOT5 is described as a complementary high-fidelity tokamak-centric Monte Carlo tool. FIRM3D’s distinguishing emphasis is on stellarators, Boozer-coordinate interoperability, symplectic and Runge–Kutta integration options, CUDA-native acceleration, and energetic-particle diagnostics in a lightweight package (Paul et al., 16 May 2026).

The framework’s limitations are explicit. It models collisionless guiding-center dynamics and does not include Fokker–Planck collisions. Fields are prescribed, so there is no self-consistent energetic-particle feedback on waves. The guiding-center approximation assumes small gyroradius relative to equilibrium scales. Accuracy also depends on Boozer harmonic truncation and the choice of radial interpolation nodes. These constraints do not diminish the code’s intended use; rather, they define its role as a reduced model for rapid confinement, loss, and resonance studies in 3D magnetic configurations (Paul et al., 16 May 2026).

FIRM3D is distributed at the ColumbiaStellaratorTheory repository, documented through a dedicated documentation site, released under the MIT license, and accompanied by unit and regression tests on CPUs and GPUs via GitHub Actions. This software profile reinforces the code’s stated objective: a modern, modular, and performant platform for studying collisionless energetic-particle dynamics and Alfvén-eigenmode-induced transport in stellarators and related devices (Paul et al., 16 May 2026).