TRIMEG Gyrokinetic Code for Tokamak Plasmas

• The TRIMEG code is a gyrokinetic PIC simulation tool that uses unstructured triangular meshes to resolve complex tokamak geometries and capture plasma instabilities.
• It employs high-order finite-element discretizations and a mixed-variable scheme to accurately simulate electromagnetic and electrostatic phenomena.
• Advanced parallelization and GPU acceleration in TRIMEG yield significant speedups in particle push and grid interpolation operations, enhancing simulation efficiency.

The TRIangular MEsh based Gyrokinetic (TRIMEG) code is a suite of gyrokinetic particle-in-cell (PIC) simulation tools developed for investigating plasma microturbulence, transport, and instability phenomena in realistic tokamak geometries using high-order finite-element discretizations on unstructured triangular meshes. TRIMEG and its modern extensions address the numerical challenges of simulating electromagnetic and electrostatic instabilities, including those driven by energetic particles, across both core and edge plasma regions while accommodating complex magnetic topologies such as separatrix and X-points. These codes incorporate advanced parallelization, GPU offloading, and rigorous field-particle coupling algorithms, and have been benchmarked against international standards and experimental equilibria (Rekhviashvili et al., 2023, Lu et al., 2019, Lu et al., 30 Apr 2025, Lu et al., 2024, Lu et al., 2022, Daneri, 17 Jan 2026).

1. Triangular Mesh-Based Discretization and Code Architecture

TRIMEG employs a fully unstructured triangular mesh in the poloidal plane (R, Z) generated via Delaunay triangulation to resolve arbitrary plasma boundaries, X-points, and divertor leg geometries. The toroidal direction (φ) is discretized uniformly, supporting periodic boundary conditions. Field variables—such as poloidal flux ψ(R, Z), safety factor q(ψ), and magnetic field B(R, Z)—are interpolated to the mesh nodes, typically using 2D B-splines or high-order Hermite-type C¹ basis functions for enhanced smoothness and accuracy (Lu et al., 2024). Each mesh triangle is equipped with a local coordinate mapping and geometric data, and mesh generation supports both axisymmetric and up–down asymmetric domains, including EFIT (EQDSK) equilibria.

The code architecture follows a modular, object-oriented paradigm (e.g., TRIMEG-GKX), with separate core classes handling mesh geometry, field interpolation, marker management, and parallel data decomposition (Lu et al., 30 Apr 2025). Field equations are assembled into sparse matrices, solved using PETSc iterative solvers in distributed-memory parallel environments.

2. Gyrokinetic Model Formulations and Field Equations

TRIMEG solves the (drift-)gyrokinetic Vlasov–Maxwell system in both δf and full-f PIC schemes (Lu et al., 2022). In the δf approach, the distribution function is split as f(z, t) = f₀(z) + δf(z, t), and the marker weights w_p evolve according to

$$\frac{dw_p}{dt} = -p_p(t)\frac{d}{dt}\ln f_0(z_p(t)),$$

where p_p relates the background Maxwellian f₀ to the marker sampling density (Rekhviashvili et al., 2023, Lu et al., 2022).

Guiding-center equations of motion:

• The gyrocenter evolves under contributions from equilibrium flows (parallel streaming, drifts) and perturbed electromagnetic fields (δΦ, δA_∥), with explicit gyro-averaging.
• For electromagnetic simulations, a mixed-variable/pullback scheme splits the parallel vector potential as

$\delta A_{\|} = \delta A_{\|}^s + \delta A_{\|}^h,$

with δA_∥^s (symplectic) evolving via

$$\partial_t \delta A_{\|}^s + \partial_{\|}\delta\Phi = 0,$$

maintaining the ideal Ohm's law constraint, and δA_∥^h (Hamiltonian) iteratively solving for kinetic corrections in Ampère’s law.

The field equations include quasi-neutrality in the long-wavelength limit,

$$- \nabla \cdot \left( \sum_s \frac{q_s n_{0s} B}{\omega_{cs}} \nabla_{\perp} \delta\Phi \right) = \sum_s q_s \delta n_{s,v},$$

and Ampère's law for the parallel vector potential,

$$- \nabla^2_{\perp} \delta A_{\|} = \mu_0 \delta j_{\|,v},$$

with corrections for finite Larmor radius (FLR) and electron skin-depth, especially vital in electromagnetic regimes (Lu et al., 2024, Lu et al., 30 Apr 2025).

3. Finite-Element Schemes and Filtering Strategies

Early TRIMEG implementations use globally C¹-continuous Hermite triangular finite-element bases, with 18 DOFs per triangle and cubic-spline expansions in φ. This design ensures 5th-order convergence in field variables, facilitating accurate and low-noise electromagnetic simulations—essential for resolving multiscale interactions and minimizing unphysical field-particle cancellation (Lu et al., 2024).

The advanced TRIMEG-GKX extension incorporates a "Piecewise Field-Aligned Finite Element Method" (PFAFEM), where high-order B-spline basis functions are aligned with local magnetic field lines. Field-aligned coordinates (local Clebsch-type) are constructed per toroidal sector, allowing representational efficiency for modes with strong anisotropy along the field (Lu et al., 30 Apr 2025).

A key feature is the avoidance of Fourier/poloidal mode filtering in the poloidal direction: TRIMEG-C1 and TRIMEG-GKX operate in real space, supporting arbitrary geometry and open field lines. This filter-free approach eliminates spectral leakage and boundary artifacts, but results in larger sparse systems and stricter timestep constraints due to the need to resolve the highest k_∥ on the unstructured mesh (Lu et al., 2024).

4. Particle-in-Cell Algorithms and Parallelization Methods

The PIC algorithm in TRIMEG maintains marker arrays with five-dimensional phase space (R, φ, Z, v_∥, μ). Weight evolution, scatter/gather operations (particle-to-grid and grid-to-particle), and gyrocenter trajectory integration are implemented using high-order numerical methods (typically 4th-order Runge–Kutta) (Lu et al., 2019, Lu et al., 2024).

For fast spatial lookup, especially important for charge deposition, TRIMEG employs optimized particle-triangle mapping via bounding-box or box–triangle intermediate grids, reducing point-in-triangle search cost from O(N_p·N_tri) to O(N_p), with empirical speed-ups up to 30-36× over brute-force (Lu et al., 2019).

Parallelization combines hybrid MPI+OpenMP strategies. The “domain cloning” model is used in TRIMEG-GKX, where the grid is replicated across MPI ranks, and marker decomposition controls data distribution. Markers are advanced independently with minimal communication, and only a global reduction is needed for field solves (Lu et al., 30 Apr 2025). Large-scale parallel runs, tested up to >1,000 cores, demonstrate near-ideal scaling in marker push and field operations (Lu et al., 2019).

5. GPU Acceleration and Code Portability

Recent work benchmarks the acceleration of TRIMEG by offloading computationally intensive PIC kernels—specifically particle push and grid-to-particle interpolation—to GPUs using OpenMP target directives. The approach delivers up to ~30× kernel speedup on AMD MI300A (with ~20× at high concurrency), and up to ~400× on NVIDIA A100/H100 relative to CPU-only baselines for the particle-push kernel. End-to-end speedup in multi-node hybrid MPI–OpenMP configurations reaches 2–5× for full simulations (Daneri, 17 Jan 2026).

Implementation challenges (e.g., handling polymorphic arguments and device-resident data structures) are addressed via inlining, memory flattening, conditional compilation, and hybrid memory management. Limitations remain in offloading the particle-to-grid (scatter) operation efficiently, due to irregular memory writes, motivating further work on GPU-centric atomics and memory reductions. Portability across AMD (amdflang + USM) and NVIDIA (nvfortran) platforms is demonstrated, with compiler maturity being the principal restricting factor.

6. Physics Validation and Benchmark Results

TRIMEG and its extensions have been benchmarked extensively:

• Electrostatic ITG (Ion Temperature Gradient) mode: Core and edge simulations for DIII-D Cyclone and ASDEX Upgrade equilibria demonstrate growth rate and mode structure agreement within 10% of ORB5, with experimental boundary data handled seamlessly (Lu et al., 2019).
• Electromagnetic TAE (Toroidicity-induced Alfvén Eigenmode): Linear and nonlinear simulations in the challenging small electron skin-depth regime (e.g., d_e ~ 1.18 × 10⁻³ m) confirm growth rates and mode structure in agreement with ORB5 and GYGLES, validating the efficacy of the mixed-variable/pullback scheme (Lu et al., 2022, Lu et al., 2024, Lu et al., 30 Apr 2025).
• Multi-species, multi-mode, and realistic geometry: TRIMEG-GKX achieves linear scaling with number of species and demonstrates high-fidelity simulations of ITG, KBM (kinetic ballooning mode), TCV ITG, RSAE (reverse shear Alfvén eigenmode) in JET and ASDEX-Upgrade geometries, matching established codes (e.g., GENE, ORB5, EUTERPE) (Lu et al., 30 Apr 2025).
• Neoclassical transport and bootstrap current: The neoclassical module accurately recovers analytic transport scalings for banana, plateau, and Pfirsch-Schlüter regimes across a range of aspect ratios, with flux-surface averages computed directly on the unstructured mesh, independent of poloidal angle (Rekhviashvili et al., 2023).
• Nonlinear and mixed-scheme capability: Mixed full-f/δf schemes allow reduced noise and computational cost for energetic particle studies, maintaining correct global profile evolution and resolving transient phenomena (Lu et al., 2022).

7. Strengths, Limitations, and Future Prospects

Strengths:

• Fully unstructured poloidal meshing supports separatrix, X-point, and scrape-off layer (SOL) physics.
• High-order finite elements (C¹ and field-aligned splines) confer spectral-like accuracy crucial for electromagnetic and turbulence studies.
• Consistent implementation of the pullback/mixed-variable scheme mitigates numerical cancellation and enables robust electromagnetic operations.
• Flexible hybrid CPU/GPU parallelization leverages exascale platforms.

Limitations:

• Present neoclassical and collision physics are limited (e.g., basic electron-ion Lorentz operator, lack of full Fokker-Planck terms).
• Real-time performance bottleneck remains in the particle-to-grid deposition not yet offloaded to GPU accelerators.
• The treatment of the self-consistent radial electric field (E_nc), ion kinetic dynamics, and nonlinear multi-mode turbulence remains an avenue for ongoing development (Rekhviashvili et al., 2023, Daneri, 17 Jan 2026).

Future Directions:

• Implementation of advanced collision operators (Fokker-Planck), sheath-resolving boundary conditions, and multi-species coupling.
• Adaptive mesh refinement and noise reduction strategies targeted at steep-gradient and edge/SOL regions.
• Extension of field-aligned finite element frameworks (e.g., PFAFEM, global C¹ triangular elements) for comprehensive tokamak edge modelling.
• Complete offload of all core PIC operations to GPU and support for multi-GPU, load-balanced architectures (Daneri, 17 Jan 2026).
• Direct validation with experimental turbulence and transport data in next-generation devices using true axis-to-edge geometries.

---

References:

• (Rekhviashvili et al., 2023)
• (Lu et al., 2024)
• (Lu et al., 30 Apr 2025)
• (Lu et al., 2022)
• (Daneri, 17 Jan 2026)
• (Lu et al., 2019)