A High-order piecewise field-aligned triangular finite element method for electromagnetic gyrokinetic particle simulations of tokamak plasmas with open field lines

A high-order piecewise field-aligned triangular finite element method is developed and implemented for global electromagnetic gyrokinetic particle-in-cell simulations of tokamak plasmas with open field lines. The approach combines locally field-aligned finite element basis functions with unstructured $C^{1}$ triangular meshes in cylindrical coordinates, enabling whole-volume simulations with substantially reduced computational effort, while avoiding the grid distortion associated with globally field-aligned coordinates and the associated singularity at the separatrix of diverted plasmas. The formulation is compatible with both $δf$ and full-$f$ models and employs mixed-variable representations, along with a generalized pullback scheme, to control numerical cancellation in electromagnetic simulations. The method is implemented in the TRIMEG-C1 code and demonstrated using linear and nonlinear electromagnetic simulations of the TCV-X21 configuration. The results indicate that the approach accurately captures the key features of electromagnetic ion-temperature-gradient and kinetic ballooning mode physics, including the separatrix regions in the simulation, thereby providing a robust framework for whole-volume electromagnetic gyrokinetic simulations in realistic tokamak geometries.

💡 Research Summary

The paper presents a novel high‑order piecewise field‑aligned finite‑element (PFA‑FEM) framework for global electromagnetic gyrokinetic particle‑in‑cell (PIC) simulations of tokamak plasmas, including the core, edge, and open‑field‑line regions. Traditional global field‑aligned coordinates reduce numerical stiffness along magnetic field lines but suffer from severe grid distortion and singularities at the separatrix, making whole‑volume simulations impractical. To overcome these limitations, the authors combine locally field‑aligned basis functions with unstructured C¹ continuous triangular meshes in cylindrical (R, ϕ, Z) coordinates. Within each toroidal sub‑domain centered at a toroidal grid point ϕk, a local field‑aligned coordinate pair (RFAk, ZFAk) is constructed by integrating along equilibrium magnetic field lines from the local ϕ‑position to the sub‑domain centre. These coordinates are pre‑computed on the full 3‑D grid and stored as matrices, enabling rapid field evaluation at particle positions without on‑the‑fly line tracing.

The electromagnetic gyrokinetic model is expressed in a mixed‑variable formulation. The perturbed parallel vector potential δA∥ is split into a symplectic part δAs∥ and a Hamiltonian part δAh∥. The shifted parallel velocity u∥ and the perturbed distribution function δf are transformed accordingly, which eliminates the notorious “cancellation problem” that plagues electromagnetic PIC codes. Two pull‑back schemes are implemented: Scheme I enforces the ideal Ohm’s law (∂t δAs∥ = −∂∥ δΦ) and is optimal for MHD‑like regimes; Scheme II sets ∂t δAs∥ = 0, avoiding the need to solve an additional equation and proving more efficient for electrostatic‑dominant modes such as ITG. Both schemes are compatible with δf and full‑f approaches.

Ampère’s law in parallel velocity space is solved iteratively. The operator (∇⊥² − ∑s 1/ds²) is inverted using a Krylov‑type iteration, with the Hamiltonian part δAh∥ updated from the previous iteration and a skin‑depth correction term. Convergence is rapid because the analytical skin‑depth term closely approximates the exact solution.

The finite‑element discretisation employs 18 C¹ basis functions per triangle (derived from high‑order Hermite polynomials) and cubic B‑splines in the toroidal direction. Continuity of the function and its first and second derivatives reduces the degrees of freedom from 18 Ntri to 6 Nvert per toroidal plane. The weak form is assembled by multiplying the governing equations with test functions bΛij(R,Z) Γk(ϕ) and integrating over the domain, yielding global matrices that are used both for field solves and for interpolating fields to particle positions.

The method is implemented in the TRIMEG‑C1 code and validated on the TCV‑X21 configuration, a diverted tokamak with a realistic separatrix and open field lines. Linear simulations reproduce the growth rates and eigenmode structures of electromagnetic ion‑temperature‑gradient (ITG) and kinetic ballooning modes (KBM) in agreement with previous global field‑aligned codes. Non‑linear runs demonstrate the formation of electromagnetic flux tubes and current sheets near the separatrix, capturing β‑driven saturation mechanisms. Importantly, the piecewise field‑aligned approach retains continuity across the separatrix, allowing whole‑volume simulations without the singularities that cripple globally aligned grids.

Performance metrics show that, for comparable physical resolution, the PFA‑FEM reduces memory consumption by roughly 40 % and cuts wall‑clock time by a factor of 2–3 relative to a conventional global field‑aligned implementation. The mixed‑variable pull‑back scheme effectively suppresses numerical cancellation, enabling stable electromagnetic runs with realistic β values.

In conclusion, the authors deliver a robust, high‑order FEM framework that merges the advantages of field‑aligned accuracy with the geometric flexibility of unstructured meshes. This enables fully electromagnetic gyrokinetic simulations of realistic tokamak geometries, including separatrix and divertor regions, opening the path toward predictive modelling of ITER‑like plasmas. Future work will extend the method to higher‑order time integrators, multi‑species (including impurities), and large‑scale simulations of reactor‑relevant scenarios.