CATAPULT: A CUDA-Accelerated Timestepper for Alpha Particles Using Local Tricubics

arXiv:2604.07617 · physics.comp-ph, physics.plasm-ph · Submitted 2026-04-08 · Read on arXiv

Listen

Transcript

Introduction to the show: ident: Physics. The day in general physics: instruments, methods and results across the field.

Kai: We're looking at the paper: "CATAPULT: A CUDA-Accelerated Timestepper for Alpha Particles Using Local Tricubics".

Mira: Here's the gist of the paper. CATAPULT introduces a CUDA-Accelerated Timestepper for Alpha Particles Using Local Tricubics (CATAPULT) for use in Monte Carlo calculations of alpha particle confinement in stellarators.

Kai: Core Methodology The paper presents CATAPULT, a GPU implementation designed to be significantly faster than existing parallelized CPU implementations for Monte Carlo calculations of alpha particle confinement in stellarators.

Mira: The method involves estimating the quantity ϕ by sampling initial conditions from a birth distribution and computing an approximation Yi of the last location of interest for a particle born at Xi.

Kai: The Monte Carlo estimate is then calculated as ϕˆ = 1/n Xn i=1 g(Yi).

Title and authors: Kai: So we're looking at this paper, "CATAPULT: A CUDA-Accelerated Timestepper for Alpha Particles Using Local Tricubics." It sounds like they’re building a tool to do Monte Carlo calculations for alpha particle confinement in stellarators.

Mira: Yeah, the title says it uses a CUDA-Accelerated Timestepper and Local Tricubics. That immediately tells us this is about using a GPU for faster simulations of how alpha particles stay confined in these magnetic devices.

Lev: From my side, I'm wondering what kind of speedup they are actually talking about compared to the existing parallelized CPU implementations we use for these kinds of calculations.

Kai: They claim their GPU implementation is significantly faster than those existing parallelized CPU implementations, and it handles both equilibrium magnetic fields and Shear Alfven Waves.

Mira: That’s interesting because handling Shear Alfven Waves adds a layer of complexity to the physics they're trying to model for confinement.

Lev: So what’s the core idea behind this speedup? Is it just running things on the GPU faster, or is there something else about how they're doing the math?

Kai: The core methodology involves estimating a quantity phi by sampling initial conditions from a birth distribution and then computing an approximation Yi of the last location of interest for a particle born at Xi.

Mira: And then they calculate the Monte Carlo estimate as phi hat equals one over n Xn i equals one g(Yi). So, it's about approximating that complex quantity using sampling.

Lev: Sampling is a standard way to do Monte Carlo stuff, but how does this specific approach improve the actual particle tracing compared to what we currently have?

Kai: They use the evolution map f(x, t) where x is the initial point in phase space and f(x, t) is the position in phase space at time t given a particle position of x at time zero <ref:2604.07617#pg2,f(x, t) is the position in phase space at time t>.

Mira: And they describe the drift motion using the Littlejohn Lagrangian L0(R, R˙, v, t), which involves terms like qiA0 plus Miv B0 B0 dot R˙ minus Miv two two minus mu iB0 <ref:2604.07617#pg3,Lagrangian L0(R, R˙ , v||, t>.

Lev: They solve those equations of motion in both Boozer and Cartesian coordinates using the Dormand-Prince adaptive timestepping scheme. That tells me they're dealing with time steps that change based on the error estimates.

Kai: And for magnetic fields, they represent quantities on a regular interpolant grid in Boozer coordinates and use tricubic interpolation on this grid to recover a dense representation of the magnetic field when tracing in Cartesian coordinates.

Mira: So, they're trying to get that high-fidelity picture of the magnetic field geometry needed for accurate particle paths.

The paper's summary: Kai: So, to summarize what CATAPULT is doing, it’s introducing a CUDA-Accelerated Timestepper for Alpha Particles Using Local Tricubics to speed up Monte Carlo calculations of alpha particle confinement in stellarators.

Mira: Exactly. It’s about making those Monte Carlo simulations much faster than the parallelized CPU versions we currently use, while also handling more complex magnetic field structures like Shear Alfven Waves.

Lev: What does that mean for the actual physics simulation? Does this speedup allow us to do something fundamentally different in terms of what we can test?

Kai: It means they've implemented it on a GPU, and the runtime results show a speedup over one hundred twenty-eight CPU cores on a single Perlmutter node is at least five times faster for various vacuum configurations in both Cartesian and Boozer coordinates.

Mira: A five times speedup is substantial, especially when you have to run these simulations repeatedly to get statistically significant results for confinement.

Lev: Five times faster, that’s good overhead reduction on the hardware side. But what about the accuracy trade-off? Are we sacrificing fidelity for that speed?

Kai: They're trying to maintain correctness while getting that speedup, and they note that CATAPULT threads share a single copy of the interpolant grid data, which allows it to scale to higher grid resolutions without running out of memory.

Mira: That shared data structure is key there. It lets them use finer magnetic field representations without hitting memory limits on the GPU side.

Lev: So they’re managing both computational intensity and memory usage simultaneously, which is always a challenge when pushing things onto accelerators like CUDA.

Kai: They tune the adaptive step size of DP5 using error estimates defined by the formula error equals max i ei atolici plus rtoli(xi + dtx˙i).

Mira: That specific error definition is crucial because it shows how they are controlling the time step acceptance criteria for stability and accuracy.

The paper's improvements: Kai: Now let's look at what the authors suggest are the specific improvements in CATAPULT over what came before.

Mira: They focus on two main areas: improving the magnetic field representation and making sure the timestepping is robust, especially near coordinate singularities.

Lev: What exactly is that improvement in magnetic field representation? Is it just using tricubic interpolation, or something more fundamental?

Kai: They use local tricubics interpolation on a regular grid in Boozer coordinates to recover a dense representation of the magnetic field when tracing in Cartesian coordinates, which enforces continuity across cell boundaries.

Mira: Enforcing continuity is important because when you're tracing particles near flux surface transitions, that smooth field representation is what gives you an accurate path.

Lev: And they address the coordinate singularity issue in Boozer coordinates by transforming to a pseudo-Cartesian coordinate system with the mapping s, theta, zeta, v to s cos theta, s sin theta, zeta, v.

Kai: They also account for the relative scaling of v parallel by multiplying atol by vtotal or by rescaling v parallel in the equations of motion when they use that transformation.

Mira: That part about accounting for the velocity scaling is necessary because when you change coordinates, you have to be careful not to lose physical information about how fast a particle is moving along the field lines.

Conclusion: Kai: So, wrapping up CATAPULT, this paper shows they can get a CUDA-Accelerated Timestepper for Alpha Particles Using Local Tricubics that significantly speeds up Monte Carlo calculations and handles both equilibrium magnetic fields and Shear Alfven Waves.

Mira: The main implications are the five times speedup over one hundred twenty-eight CPU cores and their ability to handle both Cartesian and Boozer coordinates efficiently, with memory scaling well.

Lev: For me, the real thing is that this shows how deterministic particle dynamics can be run on real hardware at a high throughput, which opens the door for things like full orbit tracing or wall load studies.

Kai: They leave room to tune to other GPUs or porting to other GPU libraries for future hardware optimization.

Mira: It leaves ample avenues for exploring charged particle motion and understanding stellarator design by allowing more complex physics than was previously feasible with these simulation methods.

Department of Statistics and Data Science, Cornell University · Department of Applied Physics and Applied Mathematics, Columbia University

physics.comp-ph, physics.plasm-ph

Submitted: 2026-04-08

Updated: 2026-10-29

License: http://arxiv.org/licenses/nonexclusive-distrib/1.0/

Importance score: 87/100

The gist: CATAPULT introduces a CUDA-Accelerated Timestepper for Alpha Particles Using Local Tricubics (CATAPULT) for use in Monte Carlo calculations of alpha particle confinement in stellarators.

Key concepts

CUDA-Accelerated Timestepper: This is a program designed to run on a Graphics Processing Unit (GPU) using NVIDIA's CUDA platform. It speeds up the time-stepping process in simulations, making it much faster than traditional CPU calculations for tracking particle movement over time.
Local Tricubics Interpolation: This technique is used to represent magnetic field data on a grid. By using tricubic interpolation, the system creates a dense and continuous representation of the magnetic field, which is essential for accurately modeling how particles move near boundaries or transitions in the stellarator's magnetic structure.
Dormand-Prince Adaptive Timestepping: This is a numerical scheme used to solve equations of motion. It automatically adjusts the size of each time step based on an estimated error. This adaptive nature ensures high accuracy by taking smaller steps when the particle dynamics change rapidly and larger steps when they are steady.

Terminology

Fidelity — 16× in this paper · explained in 2 episodes
Core — 10× in this paper · explained in 2 episodes
Accuracy — 3× in this paper · explained in 1 episode
National Energy Research Scientific Computing Center — 2× in this paper
Trace — 2× in this paper · explained in 2 episodes
Consistency — 1× in this paper · explained in 3 episodes

Summary

CATAPULT introduces a CUDA-Accelerated Timestepper for Alpha Particles Using Local Tricubics (CATAPULT) for use in Monte Carlo calculations of alpha particle confinement in stellarators.

Core Methodology

The paper presents CATAPULT, a GPU implementation designed to be significantly faster than existing parallelized CPU implementations for Monte Carlo calculations of alpha particle confinement in stellarators. The method involves estimating the quantity ϕ by sampling initial conditions from a birth distribution and computing an approximation Yi of the last location of interest for a particle born at Xi. The Monte Carlo estimate is then calculated as ϕˆ = 1/n Xn i=1 g(Yi).

Particle Orbit Model

The deterministic particle dynamics are described by the evolution map f(x, t) where x is the initial point in phase space and f(x, t) is the position in phase space at time t given a particle position of x at time 0. The drift motion of a particle with mass Mi and charge qi in an equilibrium magnetic field B0 is described by the Littlejohn Lagrangian L0(R, R˙, v, t). CATAPULT solves these equations of motion in both Boozer and Cartesian coordinates using the Dormand-Prince adaptive timestepping scheme.

Magnetic Field Representation

The magnetic field quantities in Boozer coordinates are represented on a regular interpolant grid, and CATAPULT uses tricubic interpolation on this grid to recover a dense representation of the magnetic field. When tracing in Cartesian coordinates, the authors use a regular grid in the cylindrical coordinates (R, ϕ, z).

Performance and Results

The GPU implementation of CATAPULT is significantly faster than existing parallelized CPU implementations and handles both equilibrium magnetic fields and Shear Alfven Waves. Runtime results show that the speedup over 128 CPU cores on a single Perlmutter node is at least 5 times faster for various vacuum configurations in both Cartesian and Boozer coordinates. The authors note that CATAPULT threads share a single copy of the interpolant grid data, which allows it to scale to higher grid resolutions without running out of memory.

Numerical Details

The adaptive step size of DP5 is tuned using error estimates defined by the formula error = max i ei atolici + rtoli(xi + dtx˙i). To address coordinate singularity issues in Boozer coordinates, CATAPULT transforms to a pseudo-Cartesian coordinate system with the mapping (s, θ, ζ, v∥) → (s cos θ, s sin θ, ζ, v∥). The relative scaling of v∥ is accounted for by multiplying atol by vtotal or by rescaling v∥ in the equations of motion.

Future Work

The higher throughput of particle tracing allows the consideration of more complicated physics, such as full orbit tracing or collisions, and opens the door to new work that requires more samples, such as studying wall loads for alpha particles or directly optimizing confinement. There is also the possibility of tuning to other GPUs other than A100s, porting to other GPU libraries, or optimizing for new generations for hardware. The implementation leaves ample avenues for further exploration of charged particle motion and understanding of stellarator design.

Acknowledgements

The research acknowledges funding through the U.S. Department of Energy under contracts DE-SC0024630, DE-SC0024548, DE-AC02-09CH11466, and DE-AC52 07NA27344. Funding was also provided through the Simons Foundation collaboration “Hidden Symmetries and Fusion Energy,” Grant No. 601958. This research used resources of the National Energy Research Scientific Computing Center (NERSC) under award ERCAP0031926. The paper also cites several related works, including Bonofiglo et al. [1], Bindel et al. [2], and Landreman & Buller [14].

References

[1] P. Bonofiglo, D. Dudt, C. Swanson, Fast ion confinement in quasiaxisymmetric stellarator equilibria, Nuclear Fusion 65 (2) (2025) 026050.

[2] D. Bindel, M. Landreman, M. Padidar, Direct optimization of fast-ion confinement in stellarators, Plasma Physics and Controlled Fusion 65 (6) (2023) 065012.

[3] M. Landreman, B. Medasani, F. Wechsung, A. Giuliani, R. Jorge, C. Zhu, Simsopt: a flexible framework for stellarator optimization, Journal of Open Source Software 6 (65) (2021) 3525.

[4] C. G. Albert, S. V. Kasilov, W. Kernbichler, Accelerated methods for direct computation of fusion alpha particle losses within, stellarator optimization, Journal of Plasma Physics 86 (2) (2020) 815860201.

[5] R. B. White, M. S. Chance, Hamiltonian guiding center drift orbit calculation for plasmas of arbitrary cross section, The Physics of Fluids 27 (10) (1984) 2455–2467.

[6] E. Hirvijoki, O. Asunta, T. Koskela, T. Kurki-Suonio, J. Miettunen, S. Sipilä, A. Snicker, S. Äkäslompolo, ASCOT: Solving the kinetic equation of minority particle species in tokamak plasmas, Computer Physics Communications 185 (4) (2014) 1310–1321.

[7] D. Zarzoso, D. del Castillo-Negrete, R. Lacroix, P.-E. Bernard, S. Touzet, Transport and losses of fusion-born alpha particles in the presence of tearing modes using the new toroidal accelerated particle simulator (TAPaS), Plasma Physics and Controlled Fusion 64 (4) (2022) 044003.

[8] R. G. Littlejohn, Variational principles of guiding centre motion, Journal of Plasma Physics 29 (1) (1983) 111–125.

[9] L.-M. Imbert-Gérard, E. J. Paul, A. M. Wright, An introduction to stellarators: from magnetic fields to symmetries and optimization, SIAM, 2024.

[10] J. R. Dormand, P. J. Prince, A family of embedded Runge-Kutta formulae, Journal of computational and applied mathematics 6 (1) (1980) 19–26.

[11] Boost.Numeric.Odeint, Boost.Odeint: Solving ordinary differential equations in C++, https://www.boost.org/libs/numeric/odeint (2025).

[12] B. Medasani, F. Wechsung, M.

Improvements for AI systems

GPU acceleration for particle tracing enables higher fidelity simulations by significantly reducing stochastic error through CUDA implementation, allowing for more complex physics like full orbit tracing or wall load studies that require more samples than currently feasible.
Implementation of Local Tricubics interpolation on a grid allows the system to handle magnetic field quantities in both Boozer and Cartesian coordinates with continuity enforced across cell boundaries, which is crucial for accurately tracing particles near flux surface transitions.
The adaptive timestepper, CATAPULT using Dormand-Prince 5, is tuned by transforming error estimates to pseudo-Cartesian coordinates when in Boozer space to mitigate coordinate singularities and improve the accuracy of time step acceptance criteria.

Sources

On shear Alfv'en wave-induced energetic ion transport in optimized stellarators