cuGMEC: A High-Performance Code for Gyrokinetic-MHD Hybrid Simulation on GPUs with CUDA C++
cuGMEC is the fully reconstructed GPU version of the GMEC (Gyrokinetic-MHD Energetic-particle Code) family:
-
P. Y. Jiang, Z. Y. Liu, S. Y. Liu, J. Bao, and G. Y. Fu
Development of a gyrokinetic-MHD energetic particle simulation code. I. MHD version
Physics of Plasmas 31, 073904 (2024)
DOI: 10.1063/5.0203252 -
Z. Y. Liu, P. Y. Jiang, S. Y. Liu, L. L. Zhang, and G. Y. Fu
Development of a gyrokinetic-MHD energetic particle simulation code. II. Linear simulations of Alfvén eigenmodes driven by energetic particles
Physics of Plasmas 31, 073905 (2024)
DOI: 10.1063/5.0206762 -
S. Y. Liu, P. Y. Jiang, and G. Y. Fu
cuGMEC: A High-Performance Code for Gyrokinetic-MHD Hybrid Simulation on GPUs with CUDA C++
Computer Physics Communications 327, 110249 (2026)
DOI: 10.1016/j.cpc.2026.110249
cuGMEC solves a nonlinear gyrokinetic-MHD hybrid model with fluid electrons and gyrokinetic ions, including both thermal ions and energetic particles.
The code can be used to study energetic-particle-driven Alfvén eigenmodes, such as TAE, RSAE, and BAE. It can also be applied to drift-wave and electromagnetic microinstabilities, such as ITG and KBM.
For details, please refer to our CPC paper.
The MHD component includes:
- Gyrokinetic vorticity equation
- Gyrokinetic Poisson equation
- Parallel Ampère’s law
- Parallel Ohm’s law
- Electron continuity equation
- Electron isothermal condition
Thermal ions and energetic particles are modeled gyrokinetically and advanced using the delta-f particle-in-cell method. PIC-computed pressures enter the gyrokinetic vorticity equation through the pressure-curvature term.
cuGMEC uses:
- Shifted metric coordinates
- Delta-f method for particle simulation
- Five-point central finite-difference scheme for spatial discretization
- Fourth-order Runge-Kutta method for time integration
For details, please refer to our CPC paper.
The performance reported in our CPC paper no longer represents cuGMEC's best performance, as both the MHD and PIC components have since been optimized. The updated single-GPU speed tests below are based on the ITPA TAE benchmark case, using a 256x64x16 grid and dt = 0.025 Alfvén time (major radius divided by Alfvén velocity), and we believe the results are quite fast among codes of the same type. The MHD component uses double precision; the PIC component uses 400 keV fast ions with a maximum velocity of about 1.2 times the Alfvén velocity and is tested in both double and float precision. Note that particle deposition is performed in double precision in both double- and float-precision PIC runs. GYRO denotes the number of gyro-average points, and P/G denotes the average number of particles per grid. For convenience, all GPUs are tested using the same gridDim and blockDim, so the timings may differ from the absolute optimum.
Click the triangle below for results.
NVIDIA GeForce RTX 4090 D
The MHD per-step time is 14.6ms. The PIC per-step times are shown in the table.
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
NVIDIA A800-SXM4-80GB
The MHD per-step time is 13.8ms. The PIC per-step times are shown in the table.
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
NVIDIA RTX PRO 6000 Blackwell Server Edition
The MHD per-step time is 9.76ms. The PIC per-step times are shown in the table.
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
NVIDIA B200
The MHD per-step time is 9.89ms. The PIC per-step times are shown in the table.
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
The total runtime of a simulation can be estimated as (MHD per-step time + PIC per-step time) × the total number of time steps. For other numerical parameters or multi-GPU runs, the runtime can be estimated by simple scaling.
For example, an ITER steady-state full-torus case up to toroidal mode number 36 uses a 384x32x576 grid (about 7.1 million points), P/G = 240 (about 5.1 billion particles across three ion species), GYRO = 4, dt = 0.025 Alfvén time, 20,000 steps, double-precision MHD, and float-precision PIC. Using only linear scaling from the tables above gives 5.31 hours on 8 NVIDIA A800-SXM4-80GB GPUs, compared with a measured 4.97 hours, a 6% difference; on 8 NVIDIA B200 GPUs, it gives 2.32 hours, compared with a measured 2.23 hours, a 4% difference. This agreement indicates that the tables above capture performance from the simplest benchmarks to more demanding cases.
A typical cuGMEC workflow is:
- Compute the tokamak equilibrium with scripts/equilibrium/compute2D_0170.ipynb, and export it with scripts/equilibrium/output2D_0170.ipynb.
- Generate the cuGMEC input files with scripts/preprocess/generateInput2D.m. If phase-space diagnostics are needed, also run scripts/preprocess/generatePhaseSpaceMapping2D.m.
- Configure src/cuGMEC_param.h according to docs/en-US/parameters.md.
- Configure the repository-root Makefile for the target environment according to docs/en-US/environment.md, and then compile cuGMEC.
- Run the simulation locally or submit it to a computing cluster.
- Visualize the simulation output with scripts/postprocess/visualizeMHD.m and scripts/postprocess/visualizePIC.m.
Part of our effort is going into extending cuGMEC to support non-axisymmetric devices, such as stellarators.
Other future directions remain open.
cuGMEC is released under the GNU General Public License v3.0.
Using cuGMEC generally requires nontrivial preprocessing and postprocessing, including equilibrium preparation, metric conversion, input-file generation, and output analysis. If you are interested in using cuGMEC, it is strongly recommended to contact the developer first.
For questions, suggestions, or collaboration, please contact: