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

The paper introduces CATAPULT, a CUDA-accelerated timestepper implemented in the f3d Python package that significantly outperforms existing CPU-based methods for simulating alpha particle confinement in stellarators under both equilibrium magnetic fields and Shear Alfven waves.

The Gist

The "CATAPULT" Paper: A Simple Explanation

Imagine you are trying to design a futuristic power plant that runs on the same energy as the sun (nuclear fusion). To make this work, you need to trap super-hot, charged particles (like tiny, angry bees) inside a magnetic cage called a stellarator. If these particles escape, they hit the walls and cool down the reaction, or worse, they damage the machine.

The problem is: How do we know if our magnetic cage is good enough?

1. The Old Way: A Slow, Lonely Hiker

To answer this, scientists use a method called Monte Carlo simulation. Think of it like this:

• You release thousands of "test particles" (our angry bees) into the magnetic cage.
• You watch them fly around to see if they stay trapped or escape.
• To get a reliable answer, you need to simulate millions of these particles.

In the past, scientists did this using powerful computer processors (CPUs). But CPUs are like a team of very smart, but slow, hikers. Even if you have 128 hikers working together, they have to walk one step at a time, checking the map, calculating the next move, and writing it down. It takes a long time to simulate enough particles to get a clear picture.

2. The New Way: The CATAPULT

The authors of this paper built CATAPULT (CUDA-Accelerated Timestepper for Alpha Particles Using Local Tricubics).

• The Hardware (The Super-Team): Instead of using slow hikers, CATAPULT uses GPUs (Graphics Processing Units). Imagine a GPU as a stadium filled with 10,000 tiny, incredibly fast robots. Instead of walking one by one, they all sprint in perfect unison.
• The Software (The Smart Map): The particles move through a complex magnetic field that changes shape. To calculate where a particle goes next, the computer has to look at a "map" of the magnetic field.
  • The old way used a rough, blocky map.
  • CATAPULT uses Local Tricubics. Imagine the map isn't made of flat squares, but of smooth, curved, 3D puzzle pieces (like a high-end video game terrain). This allows the particles to glide smoothly along the magnetic curves without getting stuck or taking wrong turns.

3. How It Works (The Analogy)

Think of the magnetic field as a giant, twisting roller coaster track made of invisible magnetic rails.

• The Particles: These are the roller coaster cars.
• The Simulation: We want to know: "If we launch 100,000 cars, how many will fly off the track?"
• The CPU (Old Way): A single person tries to calculate the path of one car, then stops to calculate the next, then the next. Even with 128 people, they are still calculating one car at a time.
• The GPU (CATAPULT): You have a massive factory where 10,000 robots calculate the path of 10,000 cars simultaneously. They share a single, perfect blueprint of the track (the magnetic field data) so they don't waste time copying maps.

4. The Results: Speed and Smarts

The paper tested CATAPULT on several different magnetic "roller coasters" (stellarator designs).

• Speed: CATAPULT is 5 to 60 times faster than the old CPU methods.
  • Analogy: If the old method took a week to simulate a year of particle movement, CATAPULT does it in a few hours.
• Accuracy: Because it's so fast, scientists can now use a much higher resolution map (more puzzle pieces). This means the simulation is more accurate, catching tiny details that the old, slower methods missed.
• Memory: The old methods needed to copy the map for every single worker (using up all the memory). CATAPULT shares one map among all the robots, allowing them to handle much bigger, more detailed maps without running out of space.

5. Why Does This Matter?

Fusion energy is the "holy grail" of clean power. But building a stellarator is incredibly expensive and difficult.

• Before: Scientists had to guess if a design would work because simulating it took too long.
• Now: With CATAPULT, they can run thousands of simulations quickly. They can tweak the magnetic cage design, run the test, and see if it works in minutes instead of days.

In a nutshell: CATAPULT is a super-fast, super-smart simulator that lets scientists test nuclear fusion designs much faster and more accurately than ever before, bringing us one step closer to unlimited clean energy.

Technical Summary

1. Problem Statement

In magnetic confinement fusion devices, specifically stellarators, high-temperature alpha particles generated by fusion reactions must be confined to maintain plasma temperature and drive subsequent reactions. Evaluating the confinement of these particles is critical for device design and optimization.

Current evaluation methods rely on Monte Carlo simulations where thousands of particle trajectories are traced using Ordinary Differential Equation (ODE) timesteppers until they either escape the plasma or reach a maximum time. While existing parallelized CPU implementations (e.g., SIMSOPT, SIMPLE, ORBIT) are functional, they face limitations in:

• Throughput: The stochastic error in Monte Carlo estimates scales as $O(n^{-1/2})$, requiring a massive number of particle samples ($n$) to achieve high-fidelity results. CPU-based tracing is often too slow to generate sufficient samples for complex physics scenarios (e.g., Shear Alfven Waves).
• Memory Scalability: Parallelized CPU implementations often require duplicating large interpolant grid data for every thread, limiting the achievable grid resolution due to memory constraints.
• Numerical Stability: Issues arise near coordinate singularities (e.g., the magnetic axis in Boozer coordinates) and when particles approach the last closed flux surface (LCFS), where standard interpolation may fail or produce undefined derivatives.

2. Methodology

The authors introduce CATAPULT, a CUDA-accelerated GPU implementation designed to drastically increase particle throughput while maintaining numerical accuracy.

A. Mathematical Model

• Dynamics: Particle motion is modeled using the Littlejohn Lagrangian, which describes drift motion in an equilibrium magnetic field ($B_0$).
• Coordinates: The code supports both Boozer coordinates (flux-aligned, convenient for stellarator geometry) and Cartesian coordinates (necessary for tracing past the LCFS and handling magnetic islands).
• Objective: The goal is to estimate the integral $\phi = \int_\Omega g(y(x))\mu(x)dx$, representing quantities like particle loss fractions or wall loads, using Monte Carlo sampling.

B. Magnetic Field Representation

• Local Tricubics: The magnetic field quantities (e.g., $|B|$, derivatives) are stored on a regular grid. Within each cell, the field is represented using tricubic interpolation on a $4 \times 4 \times 4$ grid of points.
• Continuity: This approach ensures continuity of field quantities across cell boundaries but does not guarantee differentiability at the boundaries (a known trade-off).
• Data Layout: Data is stored in row-major order with coordinates varying from slowest to fastest (e.g., $s, \theta, \zeta$). This layout is optimized for the toroidal dominance of particle motion in well-optimized stellarators, ensuring that particles moving to the next cell access contiguous memory.

C. GPU Implementation & Optimization

• Timestepping: CATAPULT uses the Dormand-Prince 5 (DP5) adaptive timestepper, matching the algorithm used in SIMSOPT for consistency.
• Thread Mapping:
  • A CUDA block of 64 threads (2 warps on an A100) is assigned a group of 8 particles.
  • Each thread handles one particle through the stages of the DP5 integrator (computing stages, interpolants, and derivatives).
  • Shared Memory: Intermediate calculation data for the 8 particles is allocated in shared memory to minimize global memory latency.
• Memory Optimization:
  • Constant Memory: Global constants (DP5 coefficients, simulation parameters) are stored in constant memory.
  • Interpolant Data: Stored in global memory but accessed using const and __restrict__ to aid compiler optimization.
  • Coalesced Access: To maximize memory throughput, $m$ consecutive threads (where $m$ is the number of interpolant values) access data for a single particle in lockstep. Due to the row-major storage, this allows the GPU to issue a single coalesced read for multiple double-precision values, significantly reducing instruction traffic.
• Numerical Fixes:
  • Extrapolation: When particles approach the LCFS ($s > 1$) in Boozer coordinates, the code uses the nearest cell's data as an extrapolation, preventing silent failures common in other codes.
  • Pseudo-Cartesian Error Estimation: To handle the coordinate singularity at the magnetic axis ($s=0$) in Boozer coordinates, error estimates are calculated in a pseudo-Cartesian system $(s \cos \theta, s \sin \theta, \zeta, v_\parallel)$ rather than raw coordinates. This prevents small physical displacements from appearing as large angular errors.
  • Scaling: The absolute tolerance ($atol$) is scaled by the total velocity ($v_{total}$) to account for the large magnitude of parallel velocity compared to spatial coordinates.

3. Key Contributions

1. High-Performance GPU Tracing: A novel CUDA implementation that achieves 5x to 60x speedup over 128-core CPU nodes (Perlmutter) depending on the physics complexity (vacuum vs. Shear Alfven Waves).
2. Memory Efficiency: By sharing a single copy of the interpolant grid across all GPU threads (unlike CPU implementations that duplicate data per thread), CATAPULT can handle much higher grid resolutions without running out of memory.
3. Robust Numerical Schemes: Introduction of pseudo-Cartesian error estimation and specific extrapolation logic to handle coordinate singularities and boundary conditions (LCFS) more robustly than existing codes.
4. Open Source Integration: The code is integrated into the firm3d Python package, making it accessible for the fusion research community.

4. Results

The authors tested CATAPULT on several stellarator configurations (ATEN, HSX, NCSX, QA, QH) scaled to ARIES-CS parameters.

• Accuracy & Convergence:
  • GPU and CPU results showed excellent consistency.
  • Loss fraction estimates converged well with 32,768 particles when grid resolution $\ge 20$ and tracing tolerance $\le 10^{-4}$.
  • Lower fidelity runs (higher tolerances) resulted in higher estimated loss fractions, highlighting the need for high-resolution tracing.
• Performance Scaling:
  • Particle Count: Runtime scales linearly with the number of particles. The GPU becomes saturated around 25,000 particles, at which point speedups reach a factor of 5–10 for vacuum tracing.
  • Fidelity: Runtime is dominated by tracing tolerance (step size) rather than grid resolution, thanks to the cache-locality of the interpolation scheme.
  • Complex Physics: In simulations involving Shear Alfven Waves (SAWs), the speedup increased to 40–60x compared to CPU nodes. This is attributed to the GPU's ability to handle the additional computational load of wave-particle interactions more efficiently than the CPU.
• Coordinate Systems: The code performed effectively in both Boozer and Cartesian coordinates, with Cartesian tracing showing similar speedup factors.

5. Significance

CATAPULT represents a significant leap forward in the computational capability for stellarator design and analysis:

• Enabling High-Fidelity Physics: The massive speedup allows researchers to run simulations with higher grid resolutions and tighter tolerances, leading to more accurate predictions of alpha particle confinement.
• Complex Scenario Analysis: The ability to run 40–60x faster makes it feasible to study complex, time-dependent physics (like Shear Alfven Waves) and perform direct optimization of stellarator configurations, which were previously computationally prohibitive.
• Scalability: The memory-efficient design allows for the use of extremely fine magnetic field grids, which is crucial for capturing small-scale magnetic islands and stochastic regions.
• Future Applications: The framework opens the door to more advanced studies, including full-orbit tracing, collisional effects, and detailed wall-load calculations, ultimately aiding in the design of viable fusion reactors.