Gyrokinetic Toroidal Code (GTC)

• GTC is a gyrokinetic particle-in-cell simulation tool designed for first-principles modeling of turbulence and wave-particle interactions in tokamaks and stellarators.
• It employs global electromagnetic gyrokinetic formulations with advanced numerical schemes, MPI/OpenMP parallelization, and flexible coordinate systems to accurately capture core and edge plasma dynamics.
• The code is rigorously validated against analytical theories and experimental benchmarks, and continues evolving with improvements in GPU acceleration and surrogate modeling for fusion research.

The Gyrokinetic Toroidal Code (GTC) is a highly specialized particle-in-cell (PIC) simulation tool designed for first-principles modeling of turbulence, wave–particle interactions, and stability phenomena in toroidal magnetic confinement devices, including tokamaks and stellarators. GTC implements global electromagnetic gyrokinetic formulations for multiple species and equilibrium geometries, combining advanced numerical algorithms, parallelization strategies, and novel physics models to address both core and edge plasma physics in fusion devices (Bao et al., 2017, Wang et al., 2015, Wang et al., 2017, Li et al., 17 May 2024, De et al., 2017, Wei et al., 2021).

1. Mathematical Formulation and Core Physics Models

GTC solves the five-dimensional (5D) gyrokinetic Vlasov–Maxwell system for ions—reduced from the six-dimensional full-f distribution—using the gyrocenter phase-space coordinates $(\mathbf{R}, v_\parallel, \mu)$ under strong magnetic field ordering $(\omega/\Omega_i \ll 1,\ k_\perp\rho_i \ll 1)$. Electrons are treated via a drift-kinetic approximation $$(k_\perp\rho_e \ll 1,\ \omega/k_\parallel v_{te} \ll 1)$$, split into adiabatic and non-adiabatic components. The evolution equation for species $a$ is: $$\frac{df_a}{dt} = \frac{\partial f_a}{\partial t} + \dot{\mathbf{R}}\cdot\nabla f_a + \dot{v}_\parallel\frac{\partial f_a}{\partial v_\parallel} = 0$$ with

$$\dot{\mathbf{R}} = \frac{1}{B_\parallel^*} \left[ v_\parallel \mathbf{B}^* + \frac{c}{q_a} \mathbf{b} \times ( \mu \nabla B_0 + q_a \nabla \phi ) \right ], \qquad \dot{v}_\parallel = -\frac{1}{m_a B_\parallel^*} \mathbf{B}^* \cdot \nabla \phi - \mu \frac{\mathbf{b}\cdot \nabla B_0}{m_a}$$

where $\phi$ is the electrostatic potential, $A_\parallel$ the parallel component of the vector potential, and $$\mathbf{B}^* = \mathbf{B}_0 + (c/q_a) \nabla \times A_\parallel \mathbf{b}$$ (Bao et al., 2017, Wang et al., 2015).

Gyrokinetic Poisson and parallel Ampère field equations are solved self-consistently: $$-\nabla_\perp \cdot \left( n_{0i} \frac{q_i^2}{T_i} \rho_i^2 \nabla_\perp \phi \right ) = \sum_a q_a \langle \delta f_a \rangle$$

$$-\nabla_\perp^2 A_\parallel = \frac{4\pi}{c} \sum_a q_a \int d^3v\, v_\parallel\, \delta f_a$$

Electron kinetics, including non-adiabatic moments $(\delta p_e^\parallel, \delta p_e^\perp)$, enter via conservative split-weight schemes that avoid skin-depth grid constraints (Bao et al., 2017, Wei et al., 2021).

2. Electromagnetic Formulation and Parallelization

GTC incorporates both electrostatic and fully electromagnetic models, permitting rigorous simulation of kinetic Alfvén waves, lower hybrid (LH) waves, tearing modes, ballooning modes, and infernal modes. The electromagnetic system handles compressional magnetic perturbations $(\delta B_\parallel)$, equilibrium current $(J_{||0})$, and advanced split-weight representations to address cancellation phenomena in high-$\beta$ regimes (Wei et al., 2021, Bao et al., 2016, Wang et al., 2015).

Spatial discretization uses staggered cylindrical or Boozer coordinates, with magnetic surfaces defined either analytically, from EFIT, or via VMEC in 3D non-axisymmetric geometries. Field quantities are allocated on poloidal–toroidal grids and split among MPI domains. GTC is highly scalable: particle–field communication uses a combination of MPI for inter-node exchange and OpenMP for intra-node threading, with demonstrated efficiency to $10^5$ cores; GPU acceleration is under development (Wang et al., 2015, Wei et al., 2021, Bao et al., 2017).

3. Advanced Numerical Schemes and PIC Implementation

GTC employs a PIC approach where the $\delta f$ perturbation is tracked as a weight $w_a = \delta f_a / f_{0a}$ on each marker. A combination of leapfrog, Runge–Kutta, and predictor–corrector schemes is used for time integration, ensuring second-order accuracy. Charge deposition and field updates utilize volume-weighted scatter/gather operations, often with "importance sampling" for efficient particle loading in regions of high field amplitude (Wang et al., 2015, Bao et al., 2016, Kuley et al., 2017).

Field solvers use finite-difference or Fourier methods for Laplacians and polarization operators. In the presence of electron skin-depth constraints, the conservative split-weight electron scheme completely eliminates the need for grid size restrictions $(\Delta x < d_e)$, without sacrificing physical fidelity (Bao et al., 2017).

4. Equilibrium and Magnetic Geometry: Axisymmetric and 3D Stellarators

Equilibria can be axisymmetric (tokamaks) or fully three-dimensional (stellarators), with data imported from EFIT, IPREQ, or VMEC. In Boozer coordinates $(\psi, \theta, \zeta)$, magnetic field and geometrical tensors are expressed as Fourier series over $\zeta$, capturing field-periodicity and discrete harmonics as required for stellarator symmetry. Metric coefficients, Jacobians, and curvature tensors are constructed via spline interpolation, with poloidal and toroidal mode resolution controlled by FFT and mode filtering (Wang et al., 2017, Wei et al., 2021).

GTC’s non-axisymmetric extension enables direct global simulation of ITG and microturbulence modes in devices such as W7-X and LHD, with benchmarking of toroidal and poloidal spectral resolution, mode–number convergence, and spatial mode structure as a function of device symmetry (Wang et al., 2017).

5. Physical Validation and Benchmark Cases

GTC has undergone extensive verification against analytical theory and experiment. Examples include:

• Kinetic Alfvén wave (KAW): frequencies and damping rates matched theory within $<1\%$ across $\Delta x/d_e$ ratios from $0.1$ to $2.0$, and for $\beta_e$ up to 50\% (Bao et al., 2017).
• Collisionless tearing mode: reproduction of Drake–Lee and Liu–Chen growth rates in cylindrical geometry, confirming both fluid and kinetic-$e$ regimes (Bao et al., 2017).
• RF wave processes (LH and ion Bernstein): recovery of analytic dispersion curves for ion plasma oscillation, ion Bernstein mode harmonics, and lower hybrid wave frequencies within $1-4\%$ (Kuley et al., 2017, Bao et al., 2016).
• ITG turbulence: full-torus simulations in W7-X and LHD showing mode symmetry, growth rates, and amplitude modulations in alignment with theory (Wang et al., 2017).
• Internal kink and kinetic-MHD modes: benchmarking against M3D-C1, NOVA-K, XTOR-K, yielding growth rates and frequencies within $10\%$ (Wei et al., 2021).
• Edge orbit loss: guiding center and fully kinetic pusher validation against analytic and experimental ion loss fractions at the DIII-D separatrix (De et al., 2017).
• Zonal flows: correct damping and residual level accordance with Rosenbluth–Hinton theory (De et al., 2017).
• Kinetic infernal modes: identification of mode localization, growth-rate scaling, and transition with magnetic shear modifications (Li et al., 17 May 2024).

6. Performance, Scalability, and Surrogate Modeling

GTC has demonstrated robust scaling on diverse supercomputing architectures—CPU, GPU, Intel Xeon Phi—leveraging 2D domain decomposition, private grid buffers, NUMA-aware particle sorting, and vectorized field/particle operations. On platforms such as Mira, Edison, and Piz Daint, full ITER-scale ITG runs with $>10^9$ particles and $30,000$ steps have maintained efficiency within $10-20\%$ of ideal scaling (Wang et al., 2015). GPU implementations achieve up to $4.7\times$ speedup for key kernels.

A database encompassing $>5000$ linear MHD equilibria has enabled the training of deep learning surrogate models for rapid prediction of growth rates and frequencies, with $<8\%$ mean relative error and $<1$ ms inference time—contrasting with the multi-hour runtime of full PIC simulations (Wei et al., 2021).

7. Future Extensions and Applications

Planned developments for GTC include:

• Incorporation of fully electromagnetic compressional perturbations (perpendicular Ampère’s law) (Bao et al., 2017).
• Addition of Coulomb collisions for comprehensive neoclassical and transport studies.
• GPU-based solver kernels for further performance gains (Bao et al., 2017, Wang et al., 2015).
• Monte-Carlo neutral–plasma coupling for SOL-edge simulations (De et al., 2017).
• Nonlinear validation against experimental turbulence and MHD-mode characteristics (Li et al., 17 May 2024).

Anticipated applications span micro-tearing and neoclassical tearing modes, radio-frequency nonlinear interaction studies, global turbulence and transport prediction, and real-time surrogate modeling for experimental validation and discharge planning.

---

GTC’s rigorously validated global gyrokinetic PIC framework, advanced parallel/computational techniques, and comprehensive physics models position it as a reference tool for predictive simulation of multiscale plasma phenomena in fusion-relevant toroidal geometries (Bao et al., 2017, Wang et al., 2015, Wei et al., 2021, Wang et al., 2017).