plasma-fusion-sim

How it is built

Architecture overview

A fusion plasma loses heat through small-scale turbulence, and how fast it loses it decides how large and expensive a reactor must be. A turbulence simulator evolves the plasma in a magnetic geometry over many turbulent eddy times and reports the heat and particle fluxes. Two model families are used. Gyrokinetic codes (GENE, CGYRO) follow the full distribution of particles over space and velocity, with the fast gyration around the field line averaged out, which leaves a five-dimensional field (three space dimensions, two velocity dimensions). Fluid codes (BOUT++) evolve a few moments (density, momentum, pressure, vorticity) of the plasma in three dimensions and are cheaper and more flexible about geometry. In both, the one idea that makes it work is a split of the problem by direction: a spectral method in the directions where the physics is local in wavenumber and periodic, finite differences or implicit solvers along the field line and in velocity, and a data transpose between the layouts that each step needs, which turns the whole system into a stencil-plus-FFT-plus-all-to-all machine.

Parts:

• Discretisation — which basis each dimension gets — Spectral methods, Finite differences, Fast Fourier transform, discrete-cosine-transform (no page yet)
• Time advance — explicit for the fast linear and nonlinear terms, implicit for stiff ones — Runge-Kutta 4, Adaptive Runge-Kutta, Stiff solvers (implicit BDF), Explicit Euler, Newton root finding
• Linear solves — field equations, Laplacian inversion, Newton-Krylov steps — gmres (no page yet), Incomplete factorization preconditioners, Multigrid, thomas-tridiagonal-solver (no page yet)
• Decomposition and communication — splitting the 5D array over ranks and re-splitting it — Halo exchange, Collective communication, tensor-layout-transforms (no page yet)
• Node kernels — the per-GPU or per-core loops — CUDA and GPU languages, kernel-autotuning (no page yet), gpu-memory-coalescing-banking (no page yet), SIMD programming, CPU caches and memory hierarchy
• Particle methods — marker-based (particle-in-cell) alternative to the continuum grid — Monte Carlo, rejection-sampling (no page yet)
• Output and restart — diagnostics, profiles, checkpoints — Chunked array storage, parallel-file-system-striping (no page yet), batch-checkpoint-restart (no page yet)

 equilibrium (shape, profiles, gradients) ─► setup ─► initial perturbation
                                                         │
        ┌───────────────── time loop (per step) ─────────▼──────────────────────┐
        │ collisionless step (explicit RK):                                     │
        │   streaming + drifts (finite diff) ─► field solve (allreduce)         │
        │   nonlinear bracket: FFT (kx,ky) ◄─ transpose ─► real space ─► FFT    │
        │ collision step (implicit): transpose ─► matrix-vector per (kx,θ)      │
        └───────────────────────────────┬───────────────────────────────────────┘
                                        ▼ fluxes, spectra, snapshots, checkpoints
        transport solver: profiles ◄── turbulent fluxes (reduced model or full code)

One run answers one question (the flux for given gradients at one radius, in a local "flux tube"). A profile solver wraps many such runs, or a fast surrogate model fitted to them, in an outer loop to predict how the whole plasma profile settles, which is why GACODE ships both the turbulence codes and a profile-evolution tool (TGYRO) and a reduced transport model (TGLF) (source: GACODE README).

How the parts fit together

Discretisation

CGYRO, the clearest documented case, treats each dimension with the method that suits it (source: GACODE CGYRO computational-performance page):

• Perpendicular plane: Fourier in the radial and binormal wavenumbers (kx, ky), so the gyroaverage becomes a multiplication by Bessel functions of k⊥ρ and the gyroradius operator is diagonal. This is Spectral methods with a plane wave per mode.
• Along the field line (θ): finite differences (Finite differences) with upwind dissipation for the streaming term; the equilibrium is not periodic here.
• Velocity space: pseudospectral expansions: Legendre in pitch angle, generalised Laguerre in energy (parameters N_XI, N_ENERGY; defaults 16 and 8 in its docs). These make the collision operator a dense matrix of modest rank.
• Nonlinearity: the E×B convolution (a Poisson bracket) is evaluated in real space by FFT with 3/2 zero-padding to prevent aliasing, using FFTW on CPU and cuFFT or rocFFT on GPU (Fast Fourier transform). Padding by 3/2 is the standard de-aliasing rule for a quadratic nonlinearity; omit it and energy piles up in the highest wavenumbers.

BOUT++ is the finite-difference counterpart: curvilinear coordinates, derivative operators on staggered or unstaggered grids, and a "shifted metric" or flux-coordinate-independent (FCI) parallel transform so the grid can follow sheared field lines (source: BOUT++ manual, "Parallel transforms"). Its design goal is that equations appear in readable form: ddt(rho) = -V_dot_Grad(v, rho) - rho*Div(v); (source: BOUT++ README), and the framework handles the grid, the operators, MPI and the time integrator.

Time advance

CGYRO uses operator splitting (source: same page). The collisionless part (parallel streaming, magnetic drifts, the instability drive, the nonlinear bracket) is advanced explicitly with RK4 or adaptive RK5(4)/RK7(6) (Runge-Kutta 4, Adaptive Runge-Kutta); its time step is limited by the fastest wave or flow crossing a grid cell (a CFL limit; its docs say the explicit scheme needs a smaller time step than a long-wavelength code). The collisional part is stiff, so it is advanced implicitly with a generalisation of Crank-Nicolson that reduces to a single matrix-vector multiply per (kx, θ) point; the matrix rank is the number of pitch-angle points times energy points times species. This is the plasma version of Stiff solvers (implicit BDF) done by a precomputed dense solve because the stiff coupling is local in space and small.

BOUT++ takes the other route: it hands the whole right-hand side to an ODE library. The manual lists solvers from a built-in RK4 to PVODE and the SUNDIALS family (CVODE, IDA, ARKODE), selected with the option solver:type, and calls CVODE the most commonly used one (source: BOUT++ manual, "Time integration"). CVODE is an implicit BDF/Adams solver with a Newton-Krylov inner loop, so its diagnostics print Newton iterations per step and linear iterations per Newton iteration, and it benefits from a preconditioner that the physics model supplies (BOUT++ has a preconditioning chapter). That is Newton root finding with a Krylov linear solver, i.e. gmres (no page yet) accelerated by something like Incomplete factorization preconditioners. Explicit Euler shows up only as a baseline or inside small subsystems.

Linear solves

Two solves occur in every step. The field solve (quasi-neutrality and Ampère's law in a gyrokinetic code) is, for each (kx, ky, θ), a tiny linear system whose right-hand side is an integral over velocity space; in CGYRO it is computed with a reduction over the velocity-space ranks (its kernel table lists an MPI_ALLREDUCE for this). In a fluid code the equivalent is the Laplacian inversion: from vorticity to the electric potential. BOUT++ ships several: a parallel cyclic solver, PETSc-based solvers, a geometric multigrid, an iterative solver for non-Boussinesq terms, and a serial tridiagonal solver (thomas-tridiagonal-solver, no page yet) (source: BOUT++ manual, "Laplacian inversion"). Multigrid is chosen when the operator is a full 3D elliptic one and the grid is large; the tridiagonal solve is O(n) and wins when the problem separates by Fourier mode.

Decomposition and communication

The distribution function is a large multi-index array, and each kernel wants a different index to be local. CGYRO defines three layouts, one for the linear terms (streaming, field), one for the collisions (velocity indices local, space distributed) and one for the nonlinear FFT (kx, ky local); it moves between them with MPI all-to-all transposes, uses all-reduces for the field solve, and its documentation names the communication cost of each kernel next to its compute cost (source: CGYRO performance page). The cost model is therefore: compute per kernel plus a transpose per layout change. This is the central trade of the system and is the reason tensor-layout-transforms (no page yet) matter: a layout change costs a global exchange (Collective communication), so the choice of process grid decides how much of that exchange stays inside a node.

BOUT++ is a simpler case: a 3D domain decomposition in x and y (z is kept local for the spectral and shifted-metric steps), and a Halo exchange of guard cells between neighbours before each derivative. Graph partitioning enters only for unstructured grids and is not the usual case here.

Node kernels

Inside a rank, the work is arithmetic on large arrays with little reuse, so the kernels are bandwidth-bound: the Roofline model is the standard check. CGYRO is written in Fortran with GPU offload by OpenACC or OpenMP, GPU-aware MPI, and has been ported from NVIDIA to AMD systems (source: CGYRO performance page). The knobs are the usual ones for CUDA and GPU languages: coalesced loads along the fastest-varying index (gpu-memory-coalescing-banking, no page yet), kernel launch shapes tuned per machine (kernel-autotuning, no page yet), SIMD programming on CPU, and layouts that fit CPU caches and memory hierarchy. Precision is a lever too: CGYRO has flags to run the upwind communication and the nonlinear term in 32-bit precision (UPWIND_SINGLE_FLAG, NL_SINGLE_FLAG), trading a little accuracy for speed and message size (source: CGYRO input parameter docs).

Particle methods

The alternative to a grid in velocity is to follow markers: a sample of particles carrying weights, advanced through the fields and deposited back on a mesh (particle-in-cell). It scales to high dimension because cost does not grow as a grid product, and handles complex geometry well, but the result carries sampling noise that falls only as 1/√N (Monte Carlo), so noise control (weight schemes, conservation fixes, sorting for cache efficiency) is its main engineering problem. Marker loading draws from the equilibrium distribution, which is where rejection-sampling (no page yet) comes in. The Eulerian (grid) codes above have no noise and exact conservation properties for the discretisation, at the price of a velocity grid.

Output and restart

A run is long (many eddy times, on thousands of devices), so it must survive failures: periodic checkpoints of the distribution function, written in parallel (Checkpoint and restart), and diagnostics (fluxes as time series, spectra, snapshots) written far more often but much smaller. Large arrays go into chunked, compressed containers (Chunked array storage, Zstandard) on a striped parallel file system (parallel-file-system-striping, no page yet), and the allocation logic that resubmits a job from its last checkpoint is batch-checkpoint-restart (no page yet).

Data and control flow

One time step of a gyrokinetic code (CGYRO-style):

1. Evaluate the field solve: velocity-space integrals of the distribution on each rank, an all-reduce, the three field potentials per (kx, ky, θ).
2. Evaluate the streaming kernel (linear terms) with finite differences along θ and add the drive from the gradients.
3. Evaluate the nonlinear term: transpose to the (kx, ky)-local layout, apply the 3/2-padded 2D FFTs, multiply in real space, transform back, transpose back.
4. Advance the explicit scheme (RK stages repeat 1-3, with the field solve each stage).
5. Collision step: transpose so velocity indices are local, apply the precomputed matrix, transpose back.
6. Every N steps: record fluxes and spectra; at a longer interval, write a checkpoint. Time goes to the nonlinear kernel (FFT plus its transposes), the collision matrix-vector product and the all-to-all communication; the CGYRO performance page tabulates the dominant operation of each kernel and its communication.

One physics answer end to end:

1. Choose a flux surface, geometry and the local gradients from an experimental or predicted profile.
2. Run a linear scan first (cheap): which wavenumbers are unstable and how fast they grow. This sets the resolution and time step for the expensive nonlinear run.
3. Run the nonlinear simulation until the fluxes reach a statistically steady state, discard the transient, average.
4. Feed the flux into the profile solver or fit a reduced model to many such runs; iterate until the profiles stop changing.

Scaling

In the order they bite:

1. Resolution in five dimensions. Cost grows with the product of the grid sizes and the number of steps, and steps grow as the grid refines. Levers: flux-tube (local) geometry instead of the full torus, fewer toroidal modes, velocity spectral methods with few points, adaptive time step.
2. Communication. The transposes are global exchanges and stop scaling when the message size per pair gets small. Levers: a decomposition that keeps the big all-to-alls inside a node or a fast sub-network, 32-bit communication, overlap with compute.
3. Per-device speed. Levers: GPU offload, fewer passes over the array, fused kernels, layouts matched to the fastest index.
4. Statistics. Turbulence needs long averages; a short run gives a noisy flux. Levers: longer runs, ensembles, and the reduced-model ladder (a fast model screens, the full code verifies).
5. IO and reliability. Checkpoints are the size of the state; at scale they dominate if too frequent. Levers: interval from expected failure rate, striping, compression of diagnostics. For measured scaling, read the CGYRO papers on petascale and GPU-portable performance listed on its documentation page rather than quoting numbers here.

Failure modes

Real implementations

The alternative: a particle-in-cell gyrokinetic code trades the velocity grid for markers: it handles geometry and high dimension more easily and has a noise floor instead of a resolution floor. The continuum approach suits cases where quiet, conservative results and clean linear-to-nonlinear convergence studies matter; the particle approach suits large-scale global runs with edge geometry. Use a fluid model (BOUT++ style) when the physics is collisional enough that moments are a good description or when edge geometry dominates and a kinetic code is too costly.

Sources

• Candy, Belli and Bravenec, A high-accuracy Eulerian gyrokinetic solver for collisional plasmas, Journal of Computational Physics (2016). doi
• Candy and Waltz, An Eulerian gyrokinetic-Maxwell solver, Journal of Computational Physics (2003). doi
• Dudson et al., BOUT++: A framework for parallel plasma fluid simulations, Computer Physics Communications (2009). doi
• Görler et al., The global version of the gyrokinetic turbulence code GENE, Journal of Computational Physics (2011). doi
• Jenko, Dorland, Kotschenreuther et al., Electron temperature gradient driven turbulence, Physics of Plasmas (2000). doi
• GACODE documentation: CGYRO and computational performance.
• BOUT++ manual: bout-dev.readthedocs.io (time integration, Laplacian inversion, parallel transforms, preconditioning); GENE.