CGYRO: Gyrokinetic-Maxwell Fusion Solver

• CGYRO is a high-fidelity simulation tool that solves the five-dimensional gyrokinetic Vlasov-Maxwell system to predict plasma microturbulence and transport.
• It employs spectral, pseudospectral, and high-order finite-difference methods with optimized MPI data decomposition to achieve efficient performance on GPU systems.
• CGYRO supports both nonlinear and linear analyses to deliver insights in transport prediction, pedestal physics, and validation of surrogate models in fusion applications.

CGYRO is a state-of-the-art, Eulerian, gyrokinetic-Maxwell solver that has become a community workhorse for first-principles simulation of plasma microturbulence in magnetically confined fusion devices. Across the literature it is described as an Eulerian gyrokinetic turbulence solver designed and optimized for collisional, electromagnetic, multiscale simulation, and as a local or flux-tube solver for the gyrokinetic Vlasov-Maxwell system with multiple kinetic species. Its use spans nonlinear transport prediction, linear stability analysis, pedestal studies, integrated profile evolution, high-performance computing optimization, and data generation for reduced models and machine-learning surrogates (Sfiligoi et al., 2023, Belli et al., 2022, Howard et al., 2024).

1. Scientific remit and governing equations

CGYRO addresses the problem of multiscale plasma turbulence, in which slow, large-scale motions of the thermal ions are coupled to rapid, small-scale electron dynamics. In one formulation, it solves the five-dimensional gyrokinetic Vlasov-Maxwell system for the perturbed gyrocenter distribution $h_s(\mathbf{R},v_\parallel,\mu,t)$, with fluctuating fields obtained from quasineutrality and parallel Ampère’s law. A representative local electromagnetic $\delta f$ form used in reactor and SPARC studies is

$$\partial_t\,\delta f_s + \dot{\mathbf R}_s\!\cdot\nabla \delta f_s + \dot v_{\parallel s}\,\partial_{v_\parallel}\delta f_s = C_{\rm Sugama}[\delta f_s] + \frac{Z_s e}{T_s}\,F_{0s}\Bigl(\partial_t\langle\chi\rangle + \dot{\mathbf R}_s\!\cdot\nabla\langle\chi\rangle\Bigr),$$

where $$\chi=\delta\phi-v_\parallel\delta A_\parallel+(\mu/B)\delta B_\parallel$$, $F_{0s}$ is the Maxwellian background, and $\dot{\mathbf R}_s$ and $\dot v_{\parallel s}$ contain the usual $\mathbf E\times\mathbf B$, magnetic-drift, and parallel-streaming terms (Rodriguez-Fernandez et al., 2024). Other descriptions emphasize the non-adiabatic distribution $h_s$, gyroaveraged potentials, and linearized collision operators, but the common structure is the same: a kinetic evolution equation coupled to gyrokinetic field equations for $\phi$, $\delta f$0, and, in fully electromagnetic configurations, $\delta f$1 (Sfiligoi et al., 2023, Hall et al., 14 Jun 2026).

The code is used in both nonlinear and linear regimes. Nonlinear applications include first-principles simulation of turbulence saturation, multichannel heat and particle transport, impurity transport, and alpha-particle effects. Linear applications include eigenmode identification in pedestal plasmas, mode classification by $\delta f$2, $\delta f$3, and parity, and generation of reference databases for reduced models (Yang et al., 16 Sep 2025, Hu et al., 29 Apr 2026). This breadth reflects a consistent modeling core rather than a single fixed operating mode.

The underlying physical scope varies by study. Some runs are fully electromagnetic and multispecies, retaining $\delta f$4, $\delta f$5, and $\delta f$6 with kinetic electrons, ions, impurities, and fast alpha particles, as in SPARC, ITER, ARC-class, and electron-scale MAST-U studies (Rodriguez-Fernandez et al., 2024, Howard et al., 2024, Hall et al., 14 Jun 2026, Speirs et al., 2024). Other databases use electrostatic, ion-scale limits, such as the impurity-transport study in DIII-D or the SAT3-NN training set (Odstrcil et al., 2020, Sar et al., 1 Apr 2026). This suggests that CGYRO is best understood as a general-purpose gyrokinetic platform whose active physics set is selected to match the target regime.

2. Numerical methods and internal structure

CGYRO adopts an Eulerian approach in phase space. In one high-level description, four of the five dimensions—two spatial directions perpendicular to the magnetic field and the two velocity-space coordinates—are handled spectrally or pseudo-spectrally using Fourier methods, while the remaining poloidal angle coordinate $\delta f$7 is discretized with high-order finite differences. A complementary algorithmic description states that the solver uses spectral discretization in $\delta f$8, pseudospectral treatment of $\delta f$9, and a conservative 5th-order upwind finite-difference scheme in $$\partial_t\,\delta f_s + \dot{\mathbf R}_s\!\cdot\nabla \delta f_s + \dot v_{\parallel s}\,\partial_{v_\parallel}\delta f_s = C_{\rm Sugama}[\delta f_s] + \frac{Z_s e}{T_s}\,F_{0s}\Bigl(\partial_t\langle\chi\rangle + \dot{\mathbf R}_s\!\cdot\nabla\langle\chi\rangle\Bigr),$$0 (Sfiligoi et al., 2023, Belli et al., 2022). Time integration is consistently described as implicit-explicit or semi-implicit, with stiff parallel streaming and collision terms treated implicitly and nonlinear advection and drifts treated explicitly (Sfiligoi et al., 2023, Rodriguez-Fernandez et al., 2024).

The code is organized around a small set of major modules, often called kernels (Sfiligoi et al., 2023):

Kernel	Role
field	computes the right-hand sides of the moment equations and invokes Maxwell solvers
stream	handles parallel streaming and $$\partial_t\,\delta f_s + \dot{\mathbf R}_s\!\cdot\nabla \delta f_s + \dot v_{\parallel s}\,\partial_{v_\parallel}\delta f_s = C_{\rm Sugama}[\delta f_s] + \frac{Z_s e}{T_s}\,F_{0s}\Bigl(\partial_t\langle\chi\rangle + \dot{\mathbf R}_s\!\cdot\nabla\langle\chi\rangle\Bigr),$$1 advection
shear	applies background shear effects
nonlinear	evaluates the pseudo-spectral convolutions via batched two-dimensional FFTs
collision	applies the linearized Fokker-Planck operator as dense matrix-vector products

Data decomposition is carried out over MPI ranks in various 2D layouts so that each kernel sees data locally in memory in an order optimized for its own computational pattern; global transposes or all-reduce calls reshuffle data between kernels (Sfiligoi et al., 2023). In the nl03 benchmark description, parallel decomposition uses two orthogonal MPI communicators: a toroidal communicator of size $$\partial_t\,\delta f_s + \dot{\mathbf R}_s\!\cdot\nabla \delta f_s + \dot v_{\parallel s}\,\partial_{v_\parallel}\delta f_s = C_{\rm Sugama}[\delta f_s] + \frac{Z_s e}{T_s}\,F_{0s}\Bigl(\partial_t\langle\chi\rangle + \dot{\mathbf R}_s\!\cdot\nabla\langle\chi\rangle\Bigr),$$2 and a velocity-species communicator spanning $$\partial_t\,\delta f_s + \dot{\mathbf R}_s\!\cdot\nabla \delta f_s + \dot v_{\parallel s}\,\partial_{v_\parallel}\delta f_s = C_{\rm Sugama}[\delta f_s] + \frac{Z_s e}{T_s}\,F_{0s}\Bigl(\partial_t\langle\chi\rangle + \dot{\mathbf R}_s\!\cdot\nabla\langle\chi\rangle\Bigr),$$3, with each time step interleaving MPI_AllToAll and MPI_AllReduce (Belli et al., 2022). The repeated appearance of communicator-sensitive layouts in performance studies indicates that communication topology is a first-order part of the algorithm, not a secondary implementation detail.

Published setups show that numerical resolution is strongly regime dependent. SPARC first-campaign nonlinear runs used, for a typical mid-radius case, $$\partial_t\,\delta f_s + \dot{\mathbf R}_s\!\cdot\nabla \delta f_s + \dot v_{\parallel s}\,\partial_{v_\parallel}\delta f_s = C_{\rm Sugama}[\delta f_s] + \frac{Z_s e}{T_s}\,F_{0s}\Bigl(\partial_t\langle\chi\rangle + \dot{\mathbf R}_s\!\cdot\nabla\langle\chi\rangle\Bigr),$$4, $$\partial_t\,\delta f_s + \dot{\mathbf R}_s\!\cdot\nabla \delta f_s + \dot v_{\parallel s}\,\partial_{v_\parallel}\delta f_s = C_{\rm Sugama}[\delta f_s] + \frac{Z_s e}{T_s}\,F_{0s}\Bigl(\partial_t\langle\chi\rangle + \dot{\mathbf R}_s\!\cdot\nabla\langle\chi\rangle\Bigr),$$5, $$\partial_t\,\delta f_s + \dot{\mathbf R}_s\!\cdot\nabla \delta f_s + \dot v_{\parallel s}\,\partial_{v_\parallel}\delta f_s = C_{\rm Sugama}[\delta f_s] + \frac{Z_s e}{T_s}\,F_{0s}\Bigl(\partial_t\langle\chi\rangle + \dot{\mathbf R}_s\!\cdot\nabla\langle\chi\rangle\Bigr),$$6, $$\partial_t\,\delta f_s + \dot{\mathbf R}_s\!\cdot\nabla \delta f_s + \dot v_{\parallel s}\,\partial_{v_\parallel}\delta f_s = C_{\rm Sugama}[\delta f_s] + \frac{Z_s e}{T_s}\,F_{0s}\Bigl(\partial_t\langle\chi\rangle + \dot{\mathbf R}_s\!\cdot\nabla\langle\chi\rangle\Bigr),$$7, $$\partial_t\,\delta f_s + \dot{\mathbf R}_s\!\cdot\nabla \delta f_s + \dot v_{\parallel s}\,\partial_{v_\parallel}\delta f_s = C_{\rm Sugama}[\delta f_s] + \frac{Z_s e}{T_s}\,F_{0s}\Bigl(\partial_t\langle\chi\rangle + \dot{\mathbf R}_s\!\cdot\nabla\langle\chi\rangle\Bigr),$$8, $$\partial_t\,\delta f_s + \dot{\mathbf R}_s\!\cdot\nabla \delta f_s + \dot v_{\parallel s}\,\partial_{v_\parallel}\delta f_s = C_{\rm Sugama}[\delta f_s] + \frac{Z_s e}{T_s}\,F_{0s}\Bigl(\partial_t\langle\chi\rangle + \dot{\mathbf R}_s\!\cdot\nabla\langle\chi\rangle\Bigr),$$9, and $$\chi=\delta\phi-v_\parallel\delta A_\parallel+(\mu/B)\delta B_\parallel$$0, with shifted-periodic shear boundary conditions and $$\chi=\delta\phi-v_\parallel\delta A_\parallel+(\mu/B)\delta B_\parallel$$1 (Rodriguez-Fernandez et al., 2024). In ITER baseline scenario studies, typical local domains were $$\chi=\delta\phi-v_\parallel\delta A_\parallel+(\mu/B)\delta B_\parallel$$2 with $$\chi=\delta\phi-v_\parallel\delta A_\parallel+(\mu/B)\delta B_\parallel$$3, 24–32 points in the global field-aligned $$\chi=\delta\phi-v_\parallel\delta A_\parallel+(\mu/B)\delta B_\parallel$$4 grid, and 8–12 energy nodes with 24 pitch angles (Howard et al., 2024). Pedestal linear calculations in NSTX employed $$\chi=\delta\phi-v_\parallel\delta A_\parallel+(\mu/B)\delta B_\parallel$$5, eight radial modes, $$\chi=\delta\phi-v_\parallel\delta A_\parallel+(\mu/B)\delta B_\parallel$$6, and $$\chi=\delta\phi-v_\parallel\delta A_\parallel+(\mu/B)\delta B_\parallel$$7 at fixed high-fidelity resolution (Yang et al., 16 Sep 2025).

3. GPU acceleration, communication bottlenecks, and portability

CGYRO’s recent computational history is closely tied to GPU systems. OpenACC directives such as $\delta f$05 were inserted around nested loops in the stream, shear, field, and collision modules, and the nonlinear kernel offloads large numbers of 2D complex-real FFTs to cuFFT on NVIDIA systems. When ported to AMD MI250X GPUs on Frontier, source-level compilation was almost unchanged after switching from nvfortran to cce, with a handful of optional directive clauses made explicit and cuFFT calls wrapped so that AMD targets dispatch to hipFFT built on rocFFT (Sfiligoi et al., 2023).

Source portability, however, did not imply performance portability. A straight compile on AMD exhibited overall runtimes 20–40% slower than on the A100 despite higher MI250X peak throughput. Profiling showed that the collision kernel was already faster on AMD, while the field and stream kernels ran significantly below expected performance and hipFFT lagged behind cuFFT. The reported remedies were explicit gang and vector_length prescriptions, loop reordering and temporary elimination, and FFT padding chosen to factor into small primes $$\chi=\delta\phi-v_\parallel\delta A_\parallel+(\mu/B)\delta B_\parallel$$8. These changes accelerated the field kernel by nearly $$\chi=\delta\phi-v_\parallel\delta A_\parallel+(\mu/B)\delta B_\parallel$$9 on MI250X, reduced some loop runtimes by 20–30%, and cut MI250X FFT time by up to 50%. On the sh03b benchmark, AMD runtime went from roughly 200 s to 120 s on 24 GPUs, while the A100 improved from 125 s to about 115 s; for the larger em04b multiscale case, AMD beat NVIDIA by 8–12% in total wall-clock time (Sfiligoi et al., 2023).

Communication-heavy behavior is a recurrent theme. In the nl03 benchmark, a single Google Cloud instance with 16 A100 GPUs connected by NVLINK/SWITCH achieved 520 ms per step, of which 300 ms was compute and 170 ms MPI communication, with 6.8 GB exchanged per GPU per step and no off-node network. The same study reported that this single-node run outperformed 8 Perlmutter Phase1 nodes, 16 Summit nodes, and 256 Cori nodes for the average-scale case, while larger problems still required multi-node execution because the nl03 case requires $F_{0s}$0 GB total GPU memory and the collision matrix must reside entirely on GPU (Belli et al., 2022).

Cloud and ensemble studies reinforce the same point.

Study	Configuration	Reported result
Single-node vs multi-node nl03 (Belli et al., 2022)	16× A100 on one NVLINK/SWITCH-connected node	520 ms/step; acceptable speed on a single node; outperformed 8 Perlmutter Phase1 nodes, 16 Summit nodes, and 256 Cori nodes
NVIDIA-to-AMD port (Sfiligoi et al., 2023)	sh03b, 24 GPUs	AMD runtime reduced from roughly 200 s to 120 s; optimized Frontier code became ~5–10% faster overall than on Perlmutter
Azure A100 benchmark (Sfiligoi et al., 2021)	16 ND A100 v4 nodes, tuned OpenMPI 4.1.1	67 s total time and 60 s communication time for the n103 benchmark; communication remained ~80–90% of runtime
XGYRO ensemble execution (Sfiligoi et al., 29 Jul 2025)	8 variations of nl03c on 32 Frontier nodes	1.5× end-to-end speedup by sharing the collisional constant tensor structure

The XGYRO work makes explicit why ensemble optimization matters. The collisional constant tensor $F_{0s}$1 dominates memory consumption, and in typical parameter sweeps it is bit-wise identical across runs. By distributing one shared copy across an ensemble rather than one copy per simulation, XGYRO reduces the memory footprint from $F_{0s}$2 to $F_{0s}$3, saving $F_{0s}$4, and shrinks the dominant latency-bound AllReduce cost by operating each simulation on smaller communicators (Sfiligoi et al., 29 Jul 2025). This suggests that, for CGYRO, workload orchestration across ensembles can be as important as per-run kernel tuning.

4. Flux matching, profile prediction, and burning-plasma studies

A major contemporary use of CGYRO is flux-matched profile prediction. In PORTALS, the steady-state problem is posed through

$F_{0s}$5

with $F_{0s}$6. CGYRO is treated as a black box: it takes normalized logarithmic gradients $F_{0s}$7 and derived parameters $F_{0s}$8, and returns gyro-Bohm-normalized particle, heat, momentum, and exchange fluxes. Gaussian-process surrogates are then refined iteratively to reduce the number of full nonlinear evaluations (Rodriguez-Fernandez et al., 2023).

The performance gains reported for this workflow are substantial. In a DIII-D ITER Similar Shape case with 5 channels and $F_{0s}$9, PORTALS converged in 26 profile-evaluations, corresponding to 130 CGYRO simulations, whereas a TGYRO-Newton style loop using CGYRO was estimated to require $\dot{\mathbf R}_s$0–$\dot{\mathbf R}_s$1 simulations. In a SPARC PRD L-mode parameter scan with 12 scenarios, 6 radii, and 3 channels, the first scenario required 12 evaluations and the remaining 11 reused the surrogates so that the whole 12-case sweep required 71 evaluations or 426 simulations, versus $\dot{\mathbf R}_s$2 simulations naively (Rodriguez-Fernandez et al., 2023).

In projected SPARC first-campaign plasmas, the PORTALS-CGYRO framework found stiff core transport conditions, particularly when fusion gain approaches unity, and a nearly linear $\dot{\mathbf R}_s$3 trend with $\dot{\mathbf R}_s$4 attained at $\dot{\mathbf R}_s$5–$\dot{\mathbf R}_s$6 kPa. The same study reported density peaking $\dot{\mathbf R}_s$7–$\dot{\mathbf R}_s$8 for $\dot{\mathbf R}_s$9–0.2, impurity optimization effects that raised $\dot v_{\parallel s}$0 from approximately 0.29 to approximately 0.40 as $\dot v_{\parallel s}$1 increased from 1.5 to 2.5, and validation of TGLF SAT3 with global metrics reproduced to roughly 10–20% depending on quantity and radius (Rodriguez-Fernandez et al., 2024).

In ITER baseline scenario conditions, nonlinear gyrokinetic profile prediction using CGYRO and surrogate modeling produced $\dot v_{\parallel s}$2 keV, $\dot v_{\parallel s}$3 keV, $\dot v_{\parallel s}$4, $\dot v_{\parallel s}$5 MW, $\dot v_{\parallel s}$6, $\dot v_{\parallel s}$7 s, and $\dot v_{\parallel s}$8. The study identified ITG-dominated core turbulence with $\dot v_{\parallel s}$9 at all radii and concluded that pure D and 50/50 D-T produced statistically identical fluxes and profiles, implying no core isotope effect around the ITER baseline conditions that were examined (Howard et al., 2024).

ARC-class simulations extend the same framework to fast alpha particles. At $\mathbf E\times\mathbf B$0, CGYRO predicted a low-$\mathbf E\times\mathbf B$1 fast-ion-destabilized electromagnetic mode, a reduction of the saturated DT ion heat flux from $\mathbf E\times\mathbf B$2 to $\mathbf E\times\mathbf B$3, and a nonlinear upshift of the ITG zero-flux gradient by approximately 25% relative to a thermalized-alpha control. The reported mechanism combines a slow-growing alpha-driven mode, enhanced zonal-flow energy $\mathbf E\times\mathbf B$4, and suppression of ITG eddies, with the effect confined to the radial region containing significant fast-particle density (Hall et al., 14 Jun 2026).

5. Experimental interpretation, pedestal physics, and verification

CGYRO is routinely used to identify experimentally observed modes and to connect fluctuation measurements to transport mechanisms. In nitrogen-seeded EAST H-mode with a full metal wall, linear CGYRO simulations at $\mathbf E\times\mathbf B$5 identified the coherent pedestal-foot mode as a dissipative trapped-electron mode. The dominant simulated mode peaked at $\mathbf E\times\mathbf B$6, corresponding to $\mathbf E\times\mathbf B$7, with $\mathbf E\times\mathbf B$8 and predominantly electrostatic character, and it matched the measured $\mathbf E\times\mathbf B$9 and $h_s$0–$h_s$1 kHz frequency range within the reported uncertainties (Hu et al., 29 Apr 2026).

Pedestal applications also include large linear databases. In NSTX, CGYRO served as the ground-truth solver for 864 local equilibria in an ELMy-H pedestal database, with 706 KBM-dominant, 89 MTM-dominant, and 69 TEM-dominant cases identified from growth rates, real frequencies, and field parity. The same study used fixed high-fidelity CGYRO resolution while Bayesian optimization tuned the reduced GFS model, thereby making CGYRO the reference against which reduced pedestal eigenmode solvers were evaluated (Yang et al., 16 Sep 2025). In a separate integrated pedestal workflow, linear CGYRO and QLGYRO were combined with ELITE and GATO to map KBM transport and ideal-MHD limits; for SPARC-like scans, the quasilinear electron heat flux at the top of the KBM-unstable band increased from approximately $h_s$2 at $h_s$3 T to $h_s$4 at $h_s$5 T, suggesting that access to KBM second stability may become more difficult in high-field devices (McClenaghan et al., 24 Jun 2026).

Diagnostic design is another prominent use. In MAST-U, local electron-scale CGYRO simulations at $h_s$6 supplied the fluctuation spectra for a synthetic high-$h_s$7 scattering diagnostic. The simulations showed ETG growth peaking around $h_s$8 and nonlinear electron heat fluxes $h_s$9–0.5, while the synthetic diagnostic predicted measurable $\phi$0 in the core and $\phi$1 near the pedestal, with localization lengths from $\phi$2 m to $\phi$3 m depending on radius and wavenumber (Speirs et al., 2024).

Cross-code verification has exposed both agreement and unresolved discrepancies. In KSTAR L-mode, CGYRO and GENE agreed closely on linear frequencies, growth rates, fractional changes in energy flux due to fast ions, movement of the dominant flux hot spot toward $\phi$4, and the appearance of an even-parity X-shaped phase-angle structure in $\phi$5 space. At the same time, the absolute level of thermal transport in the no-fast-ion, no-flow case differed by about a factor of three, and the effect of rotation on thermal transport was qualitatively different between the two codes in some scenarios (Kim et al., 2024). In DIII-D impurity-transport studies, nonlinear CGYRO correctly predicted that aluminum diffusion increases dramatically below the linear ITG-TEM threshold, but underpredicted the magnitude of the experimental variation: the measured midradius diffusion changed from $\phi$6 to $\phi$7, whereas CGYRO predicted approximately 4 to 10 $\phi$8 across the heating scan (Odstrcil et al., 2020). These results caution against treating CGYRO outputs as uniformly accurate across all observables and regimes.

6. Reduced models, machine learning, and the CGYRO data ecosystem

CGYRO increasingly functions as a generator of high-fidelity reference data for reduced models. In pedestal studies, QLGYRO wraps linear CGYRO outputs with the TGLF SAT1 saturation rule to produce quasilinear heat-flux maps over pedestal height-width scans (McClenaghan et al., 24 Jun 2026). In NSTX pedestal validation, GFS used CGYRO eigenvalues as the benchmark and, after Bayesian optimization of its own resolution parameters, achieved 15.91% RMS deviation in growth rate and 21.11% in real frequency relative to CGYRO, compared with 26.00% and 32.60% for TGLF in the same validation framework (Yang et al., 16 Sep 2025).

Machine-learned saturation rules are a direct extension of this program. SAT3-NN was trained on a high-resolution nonlinear CGYRO database originally generated for SAT3 development. Relative to SAT3, the neural model reduced peak-location RMS percentage error from 16.98% to 13.06% and peak-amplitude RMS percentage error from 24.26% to 8.54% for the 1D saturated potential spectrum, while also yielding smaller average and RMS errors in ion- and electron-energy flux predictions and preserving the anti-gyroBohm scaling seen in TEM-dominated cases (Sar et al., 1 Apr 2026).

A more ambitious surrogate is GKFieldFlow, which uses evolving CGYRO electrostatic potential fields in $\phi$9 to predict both future fields and transport channels. The reported architecture combines a multi-resolution 3D U-Net encoder-decoder with a dilated temporal convolutional network, trained on a DIII-D L-mode ion-scale electrostatic CGYRO run with $\delta f$00 resolution and $\delta f$01 saturated-state snapshots. Reported errors were typically $\delta f$02–3% for flux RMSE, $\delta f$03 for field amplitude NRMSE, and bounded rollout errors out to $\delta f$04, with spectral content and phase coherence preserved during autoregressive inference (Ashourvan, 6 Jan 2026).

The accumulation of such studies suggests that CGYRO now occupies two linked roles. It remains a production nonlinear and linear gyrokinetic solver for first-principles transport and stability, but it also serves as a reference standard for quasilinear closures, pedestal eigenmode models, neural saturation rules, and spatio-temporal surrogates. A plausible implication is that future CGYRO-centered workflows will increasingly combine a smaller number of expensive, high-fidelity calculations with a larger layer of surrogate-assisted orchestration, while retaining CGYRO as the calibration and validation anchor.