Piecewise field-aligned finite element method for multi-mode nonlinear particle simulations in tokamak plasmas

Abstract

This paper presents a novel approach for simulating plasma instabilities in tokamak plasmas using the piecewise field-aligned finite element method in combination with the particle-in-cell method. Our method traditionally aligns the computational grid, but defines the basis functions in piecewise field-aligned coordinates to avoid grid deformation while naturally representing the field-aligned mode structures. This scheme is formulated and implemented numerically. It also applies to the unstructured triangular meshes in principle. We have conducted linear benchmark tests, which agree well with previous results and traditional schemes. Furthermore, multiple-$n$ simulations are also carried out as a proof of principle, demonstrating the efficiency of this scheme in nonlinear turbulence simulations within the framework of the finite element method.

1. Introduction

A variety of waves and instabilities in tokamak plasmas are characterised by scale separation between the parallel and perpendicular directions to the magnetic field $\textbf { B}$ , namely, $k_{||}\ll k_\perp$ , where $k_{||}$ and $k_\perp$ are the wave vectors in the parallel and perpendicular directions, respectively. To simulate these waves and instabilities, especially in the high- $n$ limit, where $n$ is the toroidal mode number, theoretical and numerical schemes have been developed including the ballooning representation (Connor et al. 1979; Cheng et al. 1985), the mode structure decomposition (MSD) approach (Zonca & Chen 2014; Lu et al. 2012, 2013, 2017 ) and the flux-coordinate independent (FCI) scheme (Hariri & Ottaviani 2013; Stegmeir et al. 2018; Michels et al. 2021). As a specific numerical treatment, the field-aligned coordinates and flux tube method have been developed for gyrofluid turbulence simulations (Beer et al. 1995) and have been shown to be an efficient scheme in gyrokinetic simulations (Lin et al. 1998; Wang et al. 2011; Chen & Parker 2007 b ). The metric procedure for flux tube treatments of toroidal geometry has been developed previously to avoid grid deformation (Scott 2001).

While the field-aligned coordinates used in most codes have been based on the finite difference scheme, their application in the framework of the finite element method (FEM) is still not available in particle simulations of magnetically confined plasmas even though many hybrid or gyrokinetic codes are based on FEM (Kraus et al. 2017; Kleiber et al. 2024; Lanti et al. 2020; Hatzky et al. 2019; Huijsmans et al. 2023; Hoelzl & Huijsmans 2023; Holderied et al. 2021). Previous studies have introduced the field-aligned discontinuous Galerkin method for anisotropic wave equations (Dingfelder & Hindenlang 2020). A partially mesh-free approach has also been proposed for representing anisotropic spatial variations along field lines (McMillan 2017). FEM has advantages in terms of conservation properties and high accuracy in particle simulations. Gyrokinetic particle codes have been developed for the studies of turbulence, Alfvén waves and energetic particle physics in both tokamaks (Lanti et al. 2020; Huijsmans et al. 2023; Lu et al. 2019 a ; Hoelzl & Huijsmans 2023) and stellarators (Jost et al. 2001; Kleiber et al. 2024). While various schemes have been developed to efficiently represent the field-aligned mode structures such as the Fourier filter in ORB5 (Lanti et al. 2020), the phase factor in EUTERPE and the high-order differential operator-based spatial filter in JOREK (Huijsmans et al. 2023), a field-aligned FEM approach has not yet been developed.

In this work, we focus on the formulation and implementation of the field-aligned finite element method. The scheme we propose is characterised by two key features:

1. (i) the computational grids are aligned in a traditional pattern without any shift;

2. (ii) the finite element basis functions are defined on piecewise field-aligned coordinates, with each basis function being continuous along the magnetic field line. When the cubic spline is used, the $C^2$ continuity is maintained.

The rest of the paper is organised as follows. In § 2, the models and equations are detailed. The numerical schemes are described in § 3. In § 4, the simulation results are listed. In § 5, we give the conclusions and outlook.

2. Models and equations

2.1. Formulation in 3-D straight field line coordinates

2.1.1. Coordinates and magnetic field

In the curvilinear coordinates $(X,Y,Z)$ , considering the magnetic field pointing in the $(Y,Z)$ direction, the field is expressed as $\textbf { B}=\nabla \psi _Y\times \nabla Z-\nabla \psi _Z\times \nabla Y$ , where $\psi _Y$ and $\psi _Z$ are flux functions and are independent of $Y$ and $Z$ , namely, $\psi _Y=\psi _Y(X)$ , $\psi _Z=\psi _Z(X)$ . The twisting of the magnetic field along $Y$ and $Z$ is described by

(2.1) \begin{eqnarray} \bar {q}\equiv \frac {B^Z}{B^Y}=\frac {\partial \psi _Z}{\partial \psi _Y}, \end{eqnarray}

where $\bar {q}=\bar {q}(X)$ , $B^Y\equiv \textbf { B}\cdot \nabla Y$ , $B^Z\equiv \textbf { B}\cdot \nabla Z$ . Defining a Clebsch coordinate

(2.2) \begin{eqnarray} A\equiv Y-\frac {Z}{\bar {q}}, \end{eqnarray}

the magnetic field can be rewritten as

(2.3) \begin{eqnarray} \textbf { B}=-\bar {q} (\nabla \psi _Y\times \nabla A) =-\bar {q}(\partial _X\psi _Y)\nabla X\times \nabla A . \end{eqnarray}

Two Clebsch coordinates systems can be chosen as follows:

1. (i) $(X,A\equiv Y-Z/\bar {q},Z)$ , which is well defined if $\bar {q}\ne 0$ ;

2. (ii) $(X,K\equiv \bar {q}Y-Z,Z)$ , which is well defined if $\bar {q}\ne \infty$ .

Other choices are possible by choosing the third coordinate as $Y$ , giving the $(X,Y,A)$ or $(X,Y,K)$ coordinates. In the following discussion, we adopt the first choice $(X,A,Z)$ since we consider the studies of the tokamak plasmas for which $\bar {q}\ne 0$ but $\bar {q}=\infty$ at the so-called ‘X’ point.

Figure 1. Grids for simulations and the local supports where the basis functions are non-zero. The linear basis functions are adopted. Grey lines represent the $Y$ and $Z$ grids, while the dashed parallelograms, which match the colour of basis function $N_{j,k}$ , represent the local supports. The overlaps between different basis functions are shown. The basis function $N_{j,k}$ overlaps with itself and two other functions $N_{j^{\prime},k^{\prime}}$ when $k^{\prime}=k$ , as shown by the fact that $N_{0,0}$ overlaps with $N_{1,0}$ and $N_{-1,0}$ . For $k^{\prime}=k-1$ , $N_{j,k}$ overlaps with the other four basis functions (the overlap between $N_{1,1}$ with $N_{j^{\prime}=-1,0,1,2,k^{\prime}=0}$ ).

2.1.2. Piecewise field-aligned finite elements and partition of unity

The simulation grids are written as $(X_i,Y_j,Z_k)$ , where $i,j,k$ are the grid indices along $X,Y,Z$ , respectively. We show the construction of the piecewise field-aligned linear basis functions in figure 1, where the grey lines indicate the $Y,Z$ grids. If cubic B-spline basis functions are adopted, each basis function covers four intervals in each direction. The red arrow indicates the magnetic field. Instead of constructing a global Clebsch coordinate $A$ , the piecewise Clebsch coordinate $A_k$ depends on the index of the subdomain in the $Z$ direction,

(2.4) \begin{eqnarray} A_k=Y-\frac {Z-Z_k}{\bar {q}}. \end{eqnarray}

This definition is similar to the shifted metric procedure for the finite difference scheme (Scott 2001), but our method is adapted for the finite element method. The piecewise field-aligned basis functions are defined in the coordinates $X,A_k,Z$ . The function $\delta \Phi$ in configuration space is represented as a linear combination of the basis functions,

(2.5) \begin{eqnarray} \delta \Phi =\sum _{i,j,k}\delta \Phi _{i,j,k} N_i(X)N_j(A_k)N_k(Z), \end{eqnarray}

where $N$ is the basis function. In figure 1, we show the linear basis functions in $(Y,Z)$ space for a given $X$ , where $N_{j,k}\equiv N_j(A_k)N_k(Z)$ . In the $Y$ direction, each basis function $N_{j,k}$ overlaps with itself and two other functions (in total three) $N_{j^{\prime},k^{\prime}}$ when $k^{\prime}=k$ , as shown by the overlap of $N_{0,0}$ with $N_{j^{\prime}=1,-1;\;k^{\prime}=0}$ . However, when $k^{\prime}\neq k$ , each basis function $N_{j,k}$ overlaps with four other basis functions $N_{j^{\prime},k^{\prime}=k-1}$ , as shown by the overlap between $N_{1,1}$ and $N_{j^{\prime}=2,1,0,-1;\;k^{\prime}=0}$ . Additionally, $N_{j,k}$ overlaps with four other basis functions when $k^{\prime}=k+1$ ; for instance, $N_{1,1}$ overlaps with $N_{j^{\prime}=2,1,0,-1;\;k^{\prime}=2}$ . In figure 1, we illustrate the basic principles using a uniform magnetic field and linear basis functions. This simplified set-up is intended to demonstrate the foundational concepts. A more general case is considered in our code implementation, allowing for varying equilibrium magnetic fields, with cubic B-splines adopted. Nevertheless, the essential components are already demonstrated in figures 1 and 2 (to be discussed in the next section), and the more general implementation can be readily derived from this set-up.

Figure 2. Tokamak geometry with the piecewise field-aligned basis functions. The directions of the magnetic field and the alignment of the basis functions are different from (left) the inner surface to (right) the outer surface. Along the toroidal direction, the red lines follow the magnetic field. The centre of local support is at the node of the traditional grid. Only the right/left half of the local support of the cubic basis function (the two middle sections in the poloidal direction) on the outer/inner surface is shown.

For each direction, the partition of unity is satisfied,

(2.6) \begin{eqnarray} \sum _i N_i(\alpha =\alpha _p)=1,\quad \alpha \in [X,A,Z], \end{eqnarray}

where $\alpha _p$ indicates a given position of a marker or a sampling point. In the three-dimensional case,

(2.7) \begin{eqnarray} &&\sum _i\sum _j\sum _k N_i(X_p)N_j(A_{p,k})N_k(Z_p) \nonumber \\ &=&\sum _j\sum _k N_j(A_{p,k})N_k(Z_p)\sum _iN_i(X_p) \nonumber \\ &=&\sum _j\sum _k N_j(A_{p,k})N_k(Z_p) \nonumber \\ &=&\sum _kN_k(Z_p)\sum _j N_j(A_{p,k}=Y_p-(Z_p-Z_k)/\bar {q}(X_p)) \nonumber \\ &=&\sum _kN_k(Z_p)=1, \end{eqnarray}

where $A_{p,k}$ denotes the coordinate $A$ at the particle location and its dependence on $k$ is specified in the subscript. We have demonstrated the partition of unity for this scheme, which has also been verified numerically in our implementation. An equivalent and intuitive proof of the partition of unity can also be found in a mesh-free scheme (McMillan 2017).

2.2. Formulation in TRIMEG-GKX for tokamak plasmas

We have modified the previous version of the TRIMEG-GKX code (Lu et al. 2021, 2023) to deal with experimental equilibria starting from the EQDSK file and to incorporate the piecewise field-aligned FEM. In a tokamak geometry with nested magnetic surfaces, the magnetic field is specified by

(2.8) \begin{eqnarray} \textbf { B}=\nabla \psi \times \nabla \phi +F\nabla \phi, \end{eqnarray}

where $\psi$ is the poloidal magnetic flux function and $F$ is the poloidal current function. We formulate this method using tokamak coordinates $(r,\phi, \theta )$ since TRIMEG-GKX deals with the core plasma. This method can be readily extended to the whole device simulation using the unstructured triangular meshes in $(R,\phi, Z)$ coordinates, which enables the studies of the edge and scrape-off layer physics. A possible application is the combination of this piecewise field-aligned FEM and the unstructured triangular mesh in TRIMEG-C1 (Lu et al. 2019 a , 2024), which is beyond the scope of this work and will be reported in the future. In tokamak coordinates, the radial-like coordinate is defined as $r=\sqrt {(\psi -\psi _{{edge}})/(\psi _{{axis}}-\psi _{{edge}})}$ and the poloidal-like coordinate as $\theta =\arctan ((Z-Z_0)/(R-R_0))$ for $\theta \in [{-}\pi /2,\pi /2)$ . The Jacobian of $r,\phi, \theta$ is defined as $J=(\nabla r\times \nabla \phi \cdot \nabla \theta )^{-1}$ . In addition, we implemented an analytical ad hoc equilibrium with concentric circular magnetic flux surfaces using a radial-like coordinate $\rho =\sqrt {(R-R_{{axis}})^2+(Z-Z_{{axis}})^2}$ , which not only recovers the ad hoc equilibrium in ORB5 (Lanti et al. 2020) but also applies to a reversed shear safety factor profile (Meng et al. 2022). Another choice for a radial-like coordinate with $r=(\psi -\psi _{{edge}})/(\psi _{{axis}}-\psi _{{edge}})$ provides an improved resolution near the plasma edge to capture the edge instabilities and will be considered in the future. In $r,\theta, \phi$ directions, we use uniform grid spacing.

The safety factor is

(2.9) \begin{eqnarray} q(r,\theta )\equiv \frac {\textbf { B}\cdot \nabla \phi }{\textbf { B}\cdot \nabla \theta } = \frac {JF}{R^2\partial _r\psi }. \end{eqnarray}

The piecewise field-aligned coordinate $\eta _{k}$ in the subdomain centred at the $\phi _k$ grid is obtained as follows:

(2.10) \begin{eqnarray} \eta _{k}(r,\theta, \phi )= \theta -\int _{\phi _k}^\phi \textrm {d}\phi ^{\prime}\frac {1}{q(r,\theta ^{\prime}(r,\phi ^{\prime}),\phi ^{\prime})}, \end{eqnarray}

where the integral is along the magnetic field line and thus $q=q(r,\theta ^{\prime},\phi ^{\prime})$ , $\theta ^{\prime}$ is determined by following the magnetic field while varying $\phi ^{\prime}$ , namely, ${\rm d}\theta ^{\prime}/{\rm d}\phi ^{\prime}=1/q(r,\theta ^{\prime},\phi ^{\prime})$ , $\phi _k$ and $\phi$ denote the integral’s starting and end points, respectively. Equation (2.10) is general and has been implemented in this work. Specifically, in the straight field line coordinates $r,\bar \theta, \bar \phi$ , we readily get

(2.11) \begin{eqnarray} \eta _{k}(r,\bar \theta, \bar \phi )=\bar \theta -\frac {\bar \phi -\bar \phi _k}{\bar q}, \end{eqnarray}

which can be constructed in other codes that use straight field line coordinates ( $q$ is constant in a magnetic flux surface) such as EUTERPE (Kleiber et al. 2024).

The piecewise field-aligned FEM can be constructed for three-dimensional geometry with a similar treatment as that shown in figure 1, but with different directions of the magnetic field and a consequent different alignment of the basis functions at various radial locations, as indicated in figure 2.

In tokamak plasmas, the pros and cons of the two choices of the field-aligned coordinates are shown in Table 1. Using the $r,\eta, \phi$ coordinates, it is possible to treat the ‘X’ point at the separatrix of the tokamak plasma where $q=\infty$ , which is one motivation of the TRIMEG code (Lu et al. 2019 a ). By choosing $(r,\eta =\bar \theta -(\phi -\bar {\phi }_k)/\bar {q},\phi )$ instead of $(r,\bar \theta, \eta = \phi -\phi _k-\bar q \bar \theta )$ , $\eta$ is more orthogonal to the other coordinate $\phi$ than to $\bar \theta$ since the magnetic field is mainly along the toroidal direction. It should be noted that the equilibrium variables along $\phi$ with fixed $r,\eta =\bar \theta -(\phi -\phi _k)/\bar {q}$ are not constant anymore. The symmetry along the coordinate $\eta =\phi -\phi _k-\bar q \bar \theta$ remains if choosing $(r,\bar \theta, \eta =\phi -\phi _k-\bar q \bar \theta )$ . Some detailed discussions can be found in the previous work and the references therein (Lu et al. 2012; Scott 2001; Zonca & Chen 2014).

Table 1. Different choices of the Clebsch coordinates. For the sake of simplicity, we choose $\phi$ as the toroidal angle but adapt $\theta$ to $\bar \theta$ so that $(r,\bar \theta, \phi )$ is a straight field line coordinate system.

2.3. Discretisation of the distribution function and marker initialisation

Following the formulation in the previous $\delta f$ work (Hatzky et al. 2019; Lanti et al. 2020; Kleiber et al. 2024; Mishchenko et al. 2019), $N_g$ markers are used with a given distribution,

(2.12) \begin{align} g(z,t)\approx \sum _{p=1}^{N_g} \frac {\delta [z_p-z_p(t)]}{J_z}, \end{align}

where $z$ is the phase space coordinate, $N_g$ is the marker number, $\delta$ is the Dirac delta function, $J_z$ is the corresponding Jacobian, $z=({\boldsymbol{R}},v_\|,\mu \equiv v_\perp ^2/(2B))$ and $\boldsymbol{R}$ is the real space coordinate. For the full $f$ model, the total distribution of particles is represented by the markers,

(2.13) \begin{align} f(z,t)=C_{\textrm{g2f}} P_{\textrm{tot}}(z,t)g(z,t) \approx C_{\textrm{g2f}} \sum _{p=1}^{N_g} p_{p,\textrm{tot}}(t) \frac {\delta [z_p-z_p(t)]}{J_z}, \end{align}

where, for the constant $C_{\textrm{g2f}}\equiv N_f/N_g$ , $N_{f/g}$ is the number of particles/markers, and $g$ and $f$ indicate the markers and physical particles, respectively. For each marker,

(2.14) \begin{align} p_{p,\textrm{tot}}(t)=\frac {1}{C_{\textrm{g2f}}}\frac {f(z_p,t)}{g(z_p,t)}=\text {const}, \end{align}

for collisionless plasmas since

(2.15) \begin{align} \frac {\textrm {d}g(z,t)}{\textrm {d}t}=0, \quad \frac {\textrm {d}f(z,t)}{\textrm {d}t}=0 . \end{align}

The expression of $P_{\textrm{tot}}(z,t)$ (and consequently, $p_{p,\textrm{tot}}$ ) can be readily obtained as

(2.16) \begin{align} P_{\textrm{tot}}(z,t) = \frac {1}{C_{\textrm{g2f}}}\frac { f(z,t)}{g(z,t)}=\frac {n_f}{\langle n_f\rangle _V}\frac {\langle n_g \rangle _V}{n_g} \frac {f_v}{g_v}, \end{align}

where $n_f$ is the density profile and $f_v$ is the distribution in velocity space, namely, the particle distribution function $f=n_f({\boldsymbol{R}})f_v(v_\parallel, \mu )$ , and $\langle \ldots \rangle _V$ indicates the volume average. There are different choices of the marker distribution functions as discussed previously (Hatzky et al. 2019; Lanti et al. 2020). In our previous work (Lu et al. 2023), the markers are randomly distributed in the toroidal direction and in the $(R,Z)$ plane, but the distribution in velocity space is identical to that of the physical particles, which leads to

(2.17) \begin{align} p_{p,\textrm{tot}}(z,t) =\frac {\phi _{\textrm {wid}}SR}{V_{\textrm{tot}}}, \end{align}

where $\phi _{ \textrm{wid}}$ is the width of the simulation domain in the toroidal direction, $S$ is the area of the poloidal cross-section, $V_{\textrm{tot}}$ is the total volume. Equation (2.17) is reduced to $p_{p,\textrm{tot}}(z,t)={n_f}{R}/({\langle n_f\rangle _V}{R_0})$ for the tokamak equilibrium with concentric circular flux surfaces. In this work, as the first scheme, markers are loaded uniformly in $r^2$ , $\theta$ and $\phi$ directions, which yields

(2.18) \begin{align} P_{\textrm{tot}}(z,t) = \frac {n_f}{\langle n_f\rangle _V V}Jr_w\theta _w\phi _w\frac {r_{ {mid}}}{r}, \end{align}

where $r_{ {mid}}=(r_{ {min}}+r_{ {max}})/2$ is the middle location of the simulation domain, $r_w$ , $\theta _w$ and $\phi _w$ are the widths in the radial, poloidal and toroidal directions. As the second scheme, markers are loaded uniformly in $r$ , $\theta$ and $\phi$ directions, which yields

(2.19) \begin{align} P_{\textrm{tot}}(z,t) = \frac {n_f}{\langle n_f\rangle _V V}Jr_w\theta _w\phi _w . \end{align}

In TRIMEG-GKX code, the scheme in (2.18) is used more often since the poloidal Fourier filter is not used and thus we use the uniform loading in space to prevent a significant drop in the number of markers per cell near the plasma edge.

2.4. Electrostatic gyrokinetic quasi-neutrality system

In this work, we adopt the electrostatic gyrokinetic model to minimise the technical complexity related to the electromagnetic model. The time step size is set small enough so that the so-called omega-H mode (Lee 1987) does not make the simulation crash.

2.4.1. Gyrocentre’s equations of motion

The gyrocentre’s equations of motion are decomposed into the equilibrium part, corresponding to that in the equilibrium magnetic field, and the perturbed part due to the perturbed field

(2.20) \begin{align} \dot{\boldsymbol{R}} & = \dot{\boldsymbol{R}}_0 + \dot{\boldsymbol{R}}_1, \\[-10pt] \nonumber \end{align}

(2.21) \begin{align} {\dot{v}}_{\|} & = {\dot{v}}_{ \|,0}+ {\dot {v}} _{ \|,1} . \\[12pt] \nonumber \end{align}

The gyrocentre’s equations of motion are consistent with previous work (Hatzky et al. 2019; Kleiber et al. 2024; Lanti et al. 2020; Mishchenko et al. 2023, 2019), but are reduced to the electrostatic limit in this work,

(2.22) \begin{align} \dot {\boldsymbol{R}}_0 &= u_\| {\boldsymbol{b}}^* + \frac {m\mu }{qB^*_\|} {\boldsymbol{b}}\times \nabla B, \\[-10pt] \nonumber \end{align}

(2.23) \begin{align} \dot v_{\|,0} &= -\mu {\boldsymbol{b}}^*\cdot \nabla B, \\[-10pt] \nonumber \end{align}

(2.24) \begin{align} \dot {\boldsymbol{R}}_1 &= \frac {{\boldsymbol{b}}}{B^*_\|}\times \nabla \langle \delta \Phi \rangle, \\[-10pt] \nonumber \end{align}

(2.25) \begin{align} \dot v_{\|,1} &= -\frac {q_s}{m_s} {\boldsymbol{b}}^*\cdot \nabla \langle \delta \Phi \rangle, \\[12pt] \nonumber \end{align}

where ${\boldsymbol{b}}^*={\boldsymbol{b}}+(m_s/q_s)v_\|\nabla \times {\boldsymbol{b}}/B_\|^*$ , ${\boldsymbol{b}}={\boldsymbol{B}}/B$ , $B_\|^*=B+(m_s/q_s)v_\|{\boldsymbol{b}}\cdot (\nabla \times {\boldsymbol{b}})$ .

The gyrocentres’ equations of motion are expressed in $(r,\phi, \theta )$ coordinates by keeping the dominant terms similar to the early treatment in the EUTERPE code (Jost et al. 2001). The normalised equations are written where the velocity is normalised to $v_N$ , $v_N=\sqrt {2T_N/m_N}$ , where $T_N$ is the temperature unit. The length is normalised to $R_N=1$ m and the time is normalised to $t_N=R_N/v_N$ . In addition, we define the reference Larmor radius as $\rho _N=m_Nv_N/(eB_{{ref}})$ since it appears with specific physics meaning related to the magnetic drift velocity and the finite Larmor radius effect, where $m_N$ is the proton mass and $B_{{ref}}$ is the reference magnetic field. The parallel velocity and the magnetic drift velocity in equilibrium are given by

(2.26) \begin{align} \dot {r}_0&= C_d\frac {F}{J}\partial _\theta B, \\[-10pt] \nonumber \end{align}

(2.27) \begin{align} \dot \theta _0&= \frac {B^\theta }{B}v_\| -C_d\frac {F}{J}\partial _r B, \\[-10pt] \nonumber \end{align}

(2.28) \begin{align} \dot \phi _0&= \frac {B^\phi }{B}v_\| +C_d\partial _r\psi (g^{rr}g^{\phi \phi }\partial _rB+g^{r\theta }g^{\phi \phi }\partial _\theta B), \\[-10pt] \nonumber \end{align}

(2.29) \begin{align} C_d &= \frac {\bar {m}_s}{\bar {e}_s}\rho _N\frac {B_{\textrm {ref}}}{B^3}(v_\|^2+\mu B), \\[6pt] \nonumber \end{align}

where $B^\alpha \equiv {\boldsymbol{B}}\cdot \nabla \alpha$ is the contra-variant component of the equilibrium magnetic field and $g^{\alpha \beta }\equiv \nabla \alpha \cdot \nabla \beta$ is the metric tensor. The equation due to the mirror force is given by

(2.30) \begin{eqnarray} \dot {v}_{\|,0}=-\frac {\mu }{JB}\partial _r\psi \partial _\theta B. \end{eqnarray}

Regarding the perturbed equations of motion, the equilibrium variables are calculated in $(r,\phi, \theta )$ coordinates, but $\partial _r\langle \delta \Phi \rangle$ , $\partial _\theta \langle \delta \Phi \rangle$ and $\partial _\phi \langle \delta \Phi \rangle$ are calculated in $r,\phi, \eta$ coordinates. Especially, ${\boldsymbol{b}}\cdot \nabla \langle \delta \Phi \rangle$ is directly calculated in $(r,\phi, \eta )$ using ${\boldsymbol{b}}\cdot \nabla ={\boldsymbol{b}}\cdot \nabla \phi \partial /\partial \phi |_{r,\eta }$ .

The $E\times B$ velocity is given by

(2.31) \begin{align} \dot {r}_1 &= -\frac {\partial _r\psi }{B^2}\frac {g^{rr}}{R^2}\partial _\phi \langle \delta \Phi \rangle +\frac {F}{JB^2}\partial _\theta \langle \delta \Phi \rangle, \end{align}

(2.32) \begin{align} \dot \theta _1 &= -\frac {\partial _r\psi }{B^2} \frac {g^{r\theta }}{R^2} \partial _\phi \langle \delta \Phi \rangle -\frac {F}{JB^2}\partial _r\langle \delta \Phi \rangle, \end{align}

(2.33) \begin{align} \dot \phi _1&= \frac {\partial _r\psi }{B^2R^2} \left (g^{rr}\partial _r\langle \delta \Phi \rangle +g^{r\theta }\partial _\theta \langle \delta \Phi \rangle \right ), \\[12pt] \nonumber \end{align}

where $\partial _r\langle \delta \Phi \rangle$ , $\partial _\theta \langle \delta \Phi \rangle$ and $\partial _\phi \langle \delta \Phi \rangle$ are in $(r,\phi, \theta )$ coordinates and need to be calculated from the Clebsch coordinates $r,\phi, \eta$ . Using the chain rule, we readily have

(2.34) \begin{align} \partial _r|_{\phi, \theta } &= \partial _r|_{\phi, \eta } +(\partial _r\eta )\partial _\eta, \\[-10pt] \nonumber \end{align}

(2.35) \begin{align} \partial _\theta |_{r,\phi } &= (\partial _\theta \eta )\partial _\eta, \\[-10pt] \nonumber \end{align}

(2.36) \begin{align} \partial _\phi |_{r,\theta } &= \partial _\phi |_{r,\eta } +(\partial _\phi \eta )\partial _\eta . \\[12pt] \nonumber \end{align}

The gyro-average $\langle \delta \Phi \rangle$ is calculated in $(r,\theta, \phi )$ coordinates. The equation due to the parallel perturbed field is

(2.37) \begin{eqnarray} \dot v_{\|,1}=-\frac {e_s}{m_s}\frac {B^\phi }{B}\partial _\phi |_{r,\eta }\langle \delta \Phi \rangle . \end{eqnarray}

2.5. Field equation

Since we focus on the electrostatic drift waves, the quasi-neutrality equation is adopted. In addition, the polarisation density is expressed in the long wavelength limit. Then, the quasi-neutrality equation is given by

(2.38) \begin{align} \nabla \cdot \left ( G_{\textrm {P}}\nabla _\perp \delta \Phi \right ) =-\sum _s e_s\delta n_s, \quad G_P=\sum _s G_s,\quad G_s=\frac {|e_s|n_s}{\omega _{c,s} B}, \end{align}

where $\omega _{c,s}=|e_s| B/m_s$ is the cyclotron frequency.

The rigorous expression in $(r,\theta, \phi )$ or $(r,\eta, \phi )$ is given by

(2.39) \begin{eqnarray} \frac {1}{J}\partial _\alpha (JG_s g^{\alpha \beta } \partial _\beta \delta \Phi ) -\frac {1}{J}\partial _\alpha (JG_s\nabla \alpha \cdot {\boldsymbol{b}}{\boldsymbol{b}}\cdot \nabla \beta \partial _\beta \delta \Phi ) =-\sum _s e_s\delta n_s, \end{eqnarray}

where $\alpha, \beta \in (r,\theta, \phi )$ for the previous scheme we have implemented (Lu et al. 2023) or $\alpha, \beta \in (r,\eta, \phi )$ for the piecewise field-aligned FEM and the Einstein summation rule is adopted. Since the poloidal cross-section is a good approximation of the perpendicular surface, the simplified expression is given by

(2.40) \begin{eqnarray} \frac {1}{J}\partial _\gamma (JG_s g^{\gamma \iota } \partial _\iota \delta \Phi ) =-\sum _s e_s\delta n_s, \end{eqnarray}

where $\gamma, \iota \in (r,\theta )$ or $\gamma, \iota \in (r,\eta )$ . In this work, both expressions are implemented numerically, but no significant differences are found in the growth rate for the base case with $n=5$ in § 4. Thus, the following studies are all based on the simplified expression of the field equation.

When the piecewise field-aligned FEM is used, the metric tensor $g^{\alpha \beta }$ is calculated using the chain rule as follows (note $\eta =\eta _k$ ):

(2.41) \begin{align} g^{r\eta }&= (\partial _r\eta ) g^{rr} + (\partial _\theta \eta ) g^{r\theta } + (\partial _\phi \eta ) g^{\phi \theta }, \\[-12pt] \nonumber \end{align}

(2.42) \begin{align} g^{\eta \eta }&= (\partial _r\eta )^2 g^{rr} + (\partial _r\eta ) (\partial _\theta \eta ) g^{r\theta } + (\partial _\theta \eta ) (\partial _r\eta ) g^{\theta r} +(\partial _\theta \eta )^2 g^{\theta \theta }, \\[-12pt] \nonumber \end{align}

(2.43) \begin{align} g^{\phi \eta }&= (\partial _r\eta ) g^{r\phi } + (\partial _\theta \eta ) g^{\phi \theta } + (\partial _\phi \eta ) g^{\phi \phi } . \\[12pt] \nonumber \end{align}

The metric tensors $g^{rr}$ , $g^{r\phi }$ and $g^{\phi r}$ in $(r,\eta, \phi )$ are the same as those in $(r,\theta, \phi )$ .

3. Numerical implementation

3.1. General description

The TRIMEG-GKX code is based on structured meshes for studying the core plasmas in tokamaks (Lu et al. 2021, 2023). It is written in Fortran. Object Oriented Programming (OOP) concepts are considered with a similar structure to the TRIMEG-C0/C1 code based on the unstructured meshes (Lu et al. 2019 a , 2024). The gyrokinetic field-particle system is decomposed into different classes, namely, equilibrium, particle, field, solver and B-spline classes. The application of the gyrokinetic field-particle classes is constructed by other basic classes. The kernel of the Fortran code is approximately14 000 lines. The PETSc library is adopted to solve linear field equations using the KSP solver. The shared memory in the MPI3 standard is used to store the three-dimensional (3-D) field with affordable memory consumption. The equilibrium variables are represented using the B-splines (Williams 2008). The FEM is implemented using cubic splines with the details provided in our previous work (Lu et al. 2023). The wedge number of the torus, $ L_{\phi, {mult}}$ , is introduced to facilitate simulations by restricting the simulated region in the toroidal direction to $ 1/{L_{\phi, {mult}}}$ of the entire torus. For simulations involving a single toroidal harmonic, $ L_{\phi, {mult}}$ can be set equal to the toroidal mode number $ n$ . In contrast, for nonlinear simulations with multiple $ n$ values, $ L_{\phi, {mult}}$ is typically set to 1 to include toroidal harmonics $ n = 0, 1, 2, 3, \dots$ . However, other positive integer values of $ L_{\phi, {mult}}$ can also be used to simulate harmonics $ n = 0, L_{\phi, {mult}}, 2L_{\phi, {mult}}, 3L_{\phi, {mult}}, \dots$ .

3.2. The 3-D field-aligned FEM solver and the mixed 2-D1F solver

For the comparison with the 3-D field-aligned FEM solver, a mixed 2-D1F solver is also developed in this work, following the previous work (Lu et al. 2019 a ), but using structured meshes and cubic splines. For the 3-D solver, the size of the grids is $(N_r,N_\theta, N_\phi )$ and $(N_{r,\textrm{FEM}},N_{\theta, \textrm{FEM}},N_{\phi, \textrm{FEM}})$ basis functions are adopted to represent functions in the simulation domain, where $N_{r,\textrm{FEM}}=N_r+\Delta N$ , $N_{\theta, \textrm{FEM}}=N_\theta$ , $N_{\phi, \textrm{FEM}}=N_\phi$ , which are consistent with the boundary conditions, where $\Delta N=2$ since cubic splines are adopted. We apply the periodic boundary conditions in the $(\theta, \phi )$ directions and implement the Dirichlet and Neumann boundary conditions in the $r$ direction. In poloidal and toroidal directions, the cubic finite element basis functions $N(x)$ are as follows:

(3.1) \begin{align} N_{\textrm{cubic}}^{ \textrm{inner}}(x) = \begin{cases} 4/3+2x+x^2+x^3/6 & \text {if $x\in [{-}2,-1),$}\\ 2/3-x^2-x^3/2 & \text {if $x\in [{-}1,0), $}\\ 2/3-x^2+x^3/2 & \text {if $x\in [0,1), $} \\ 4/3-2x+x^2-x^3/6 & \text {if $x\in [1,2) $}. \end{cases} \end{align}

Along $\theta$ and $\phi$ , the $i$ th basis function is $N_i=N_{\textrm {cubic}}(x+1-i)$ . In the radial direction, $N_i$ is the same as those in poloidal/toroidal directions as $i\geqslant 4$ or $i\leqslant N_{r,{FEM}}-3$ . The first basis function is

(3.2) \begin{align} N_{\textrm{cubic}}^{ \textrm{outer},1}(x) = \begin{cases} 0 & \text {if $x\in [{-}2,1)$,}\\ -x^3+6x^2-12x+8 \;\; & \text {if $x\in [1,2) $}. \end{cases} \end{align}

The second basis function is

(3.3) \begin{align} N_{\textrm{cubic}}^{\textrm{outer},2}(x) = \begin{cases} 0 & \text {if $x\in [{-}2,0)$},\\ 7x^3/6-3x^2+2x & \text {if $x\in [0,1) $}, \\ 4/3-2x+x^2-x^3/6 & \text {if $x\in [1,2) $} . \end{cases} \end{align}

The third basis function is

(3.4) \begin{align} N_{\textrm{cubic}}^{\textrm{outer},3}(x) = \begin{cases} 0 & \text {if $x\in [{-}2,-1)$,}\\ -x^3/3-x^2+2/3 & \text {if $x\in [{-}1,0) $,}\\ x^3/2-x^2+2/3 & \text {if $x\in [0,1) $,} \\ -x^3/6+x^2-2x+4/3 & \text {if $x\in [1,2) $} . \end{cases} \end{align}

The last three basis functions are a symmetric mapping of the first three basis functions to the middle point of the simulation domain. All radial basis functions are constructed according to $N_i=N_{ \textrm{cubic}}(x+1-i)$ , where $i\in [1,N_{r,{FEM}}]$ .

The matrix form of the quasi-neutrality equations is

(3.5) \begin{eqnarray} \bar {\bar {M}}_{\mathrm {P},L,ii^{\prime},jj^{\prime},kk^{\prime}} \cdot \delta \Phi _{i^{\prime}j^{\prime}k^{\prime}} &&= C_{\mathrm P} \delta N^{i,j,k}, \end{eqnarray}

where $\delta \Phi$ is normalised to $m_Nv_N^2/e$ and $C_P=1/\rho _N^2$ .

For the 3-D field-aligned FEM solver,

(3.6) \begin{eqnarray} \bar {\bar {M}}_{\mathrm {P},L,ii^{\prime},jj^{\prime},kk^{\prime}} &&=\sum \bar {m}_s \int \textrm {d} r\, \textrm {d}\theta \, \textrm {d}\phi \, J \tilde {N}_{ijk} \nabla \cdot \left [ n_{0s}\frac {B^2_{{ref}}}{B^2} \nabla _\perp \tilde {N}_{i^{\prime}j^{\prime}k^{\prime}}\right ], \nonumber \\ \delta N^{i,j,k} &&= \sum _s C_{ {p2g},s}\sum _{p=1}^{N_g} w_p \tilde {N}_{ijk}(r_p,\eta _{p,k},\phi _p), \end{eqnarray}

where $\tilde {N}_{ijk}(r,\eta, \phi )=N_i(r)N_j(\eta _k)N_k(\phi )$ , $C_{\textrm {p2g},s}=-\bar {q}_s \langle n\rangle _V V_{\textrm{tot}}/N_g$ , $ \langle \ldots \rangle _V$ indicates the volume average and $V_{\textrm{tot}}$ is the total volume.

For the 2-D1F solver, we have

(3.7) \begin{eqnarray} \bar {\bar {M}}_{\mathrm {P},L,ii^{\prime},jj^{\prime},kk^{\prime}} &&=\sum \bar {m}_s\int \textrm {d} r\, \textrm {d}\theta \, \textrm {d}\phi \, J N_i N_j \textrm {e}^{{\text i}k\phi } \nabla \cdot \left [ n_{0s}\frac {B^2_{\textrm {ref}}}{B^2} \nabla _\perp (N_{i^{\prime}}N_{j^{\prime}}\textrm {e}^{-{\text i}k\phi })\delta _{k,-k^{\prime}}\right ] \;\;, \nonumber \\ \delta N^{i,j,k} &&=\sum _s C_{ {p2g},s}\sum _{p=1}^{N_g} w_p N_i(r_p)N_j(\theta _p) \textrm {e}^{-{\text i}k\phi _p}, \end{eqnarray}

where $N_i=N_i(r)$ , $N_j=N_j(\theta )$ , $\delta _{i,j}=1$ if $i=j$ and $\delta _{i,j}=0$ if $i\ne j$ . The grid variables $\delta \Phi _{ijk}$ and $\delta {N}^{ijk}$ and the particles interact in the ‘scatter’ and ‘gather’ processes. The scatter operation assigns charge density and current density back to the grid. In the gather operation, field values are interpolated at the position of the particles. In the scattering process using the 3-D field-aligned FEM, the values of the basis function $\tilde {N}_{ijk}(r,\eta, \phi )$ is calculated at the particle location, but in $(r,\eta, \phi )$ coordinates, as shown in (3.6). Thus, the grid value $\delta {N}^{ijk}$ is obtained by taking into account all particles where $\tilde {N}_{ijk}(r,\eta, \phi )$ does not vanish. In the gathering process using the 3-D field-aligned FEM, as the specific application of (2.5), the field at the particle location is

(3.8) \begin{eqnarray} \delta \Phi (r_p,\theta _p,\phi _p)=\sum _{i,j,k}\delta \Phi _{i,j,k} N_i(r_p)N_j(\eta _{k,p})N_k(\phi _p). \end{eqnarray}

The B-spline coefficients of $\eta (r,\theta, \phi )$ are stored in a three-dimensional matrix in the TRIMEG code as it also applies to the open field line region when the $R,\phi, Z$ coordinates are used and the three-dimensional interpolation is evoked in both the ‘gather’ and ‘scatter’ processes. As a future optimisation for core plasma simulations using the numerical equilibrium, $\eta$ can be reduced to $\eta =\bar \theta -\phi /\bar {q}$ in the straight field line coordinate $(r,\phi, \bar \theta )$ where $\bar {q}\equiv \textbf { B}\cdot \nabla \phi /\textbf { B}\cdot \nabla \bar \theta$ is independent of $\bar \theta$ , which is more efficient as only a one-dimensional interpolation of $\bar {q}(r_p)$ is needed.

3.3. Matrix construction using the piecewise field-aligned FEM

There are different ways to calculate the matrix and stiffness matrices using the piecewise field-aligned FEM. Note that the partition of unity is always satisfied. The main difference is the numerical efficiency and accuracy in the matrix generation.

3.3.1. The rigorous numerical integral

The general form of the element of the mass/stiffness matrix is

(3.9) \begin{eqnarray} M_{ii^{\prime},jj^{\prime},kk^{\prime}}&=&\int \textrm {d}\!\mathop x^{\rightarrow}\, \partial _{d_1}N_i(X)\partial _{d^{\prime}_1}N_{i^{\prime}}(X) \nonumber \\ &\times & \partial _{d_2}N_j(A_k)\partial _{d^{\prime}_2}N_{j^{\prime}}(A_{k^{\prime}}) \partial _{d_3}N_k(Z)\partial _{d^{\prime}_3}N_{k^{\prime}}(Z), \end{eqnarray}

where $i,j,k$ indicate the row indices, $d_1$ , $d_2$ and $d_3$ indicate the differential orders in the three directions, and $^{\prime}$ indicates the variables of the matrix column. We use the cubic B-spline basis functions. Each basis function $N_i$ is defined in four sections of the meshes. In $X$ and $Z$ directions, the overlaps of two basis functions can be one, two, three or four sections. However, in the $Y$ direction, the overlap can be a fractional length of the sections, as shown in figure 1. As a result, the integral region needs to be identified first, as shown in figure 3, where $\partial _{d_2}N_j(A_k)\partial _{d^{\prime}_2}N_{j^{\prime}}(A_{k^{\prime}})$ needs to be integrated between the vertical red dashed line and the vertical blue dashed line. The Gauss–Legendre quadrature points are generated in this overlapping region and the integral is numerically calculated from the summation. This scheme is challenging to extend to unstructured triangular meshes since it is more complicated to identify the overlapping regions.

3.3.2. The mixed particle-wise–basis-wise construction

In addition to the rigorous numerical integral in § 3.3.1, the mixed particle-wise–basis-wise construction is also possible with a minor loss of accuracy. As shown in figure 3, the Gauss–Legendre points are generated along the basis function of the row in the matrix, and the corresponding values of the basis function of the column in the matrix are calculated at the same $Y$ locations. Then the integral of the term $\partial _{d_2}N_j(A_k)\partial _{d^{\prime}_2}N_{j^{\prime}}(A_{k^{\prime}})$ is calculated using the values at these points, indicated by the red circles in figure 3. This scheme can be less accurate than the rigorous integral, but is easier to extend to unstructured triangular meshes in TRIMEG-C0/C1 (Lu et al. 2019 a , 2024), as will be developed in the future.

Figure 3. The numerical integral using the Gauss–Legendre quadrature when $k^{\prime}\ne k$ , where $k$ and $k^{\prime}$ are the indices in the toroidal direction.

3.3.3. Particle-wise construction using the Monte Carlo integration and importance sampling

The mass and stiffness matrices are calculated using the Monte Carlo integration, instead of the Gaussian quadrature, to avoid the complexities of identifying the non-conforming intersection surfaces of different volumes. The integral can be generally calculated using the Monte Carlo integration in connection with the importance sampling,

(3.10) \begin{eqnarray} I\equiv \int \textrm{d}\!\mathop x^{\rightarrow} \, K(\!\mathop x^{\rightarrow})\approx \sum _i\frac {K(\!\mathop x^{\rightarrow})}{g(\!\mathop x^{\rightarrow})}, \end{eqnarray}

where $K$ is the integrand and $g(\!\mathop x^{\rightarrow})$ is the distribution of the sampling points. For uniformly distributed sampling points,

(3.11) \begin{eqnarray} I\approx \frac {V}{N}\sum _iK(\!\mathop x^{\rightarrow}). \end{eqnarray}

In general geometry, uniform loading can be adopted along three coordinates $(X,Y,Z)$ and thus,

(3.12) \begin{eqnarray} g(\!\mathop x^{\rightarrow})=\frac {N}{JX_wY_wZ_w}, \end{eqnarray}

where $J$ is the Jacobian, and $X_w$ , $Y_w$ and $Z_w$ are the width along $X$ , $Y$ and $Z$ , respectively. This scheme can be less accurate than the rigorous integral and the mixed scheme, but is easier to extend to two-dimensional unstructured triangular meshes or three-dimensional unstructured meshes.

The properties of the three methods of matrix construction are listed in Table 2. While the rigorous scheme gives high accuracy, it cannot be directly applied to the unstructured meshes adopted in TRIMEG-C1 since it is complicated to identify the overlapping regions of different triangles. The mixed and the Monte Carlo schemes apply to the unstructured meshes, but are not as accurate as the rigorous scheme.

Table 2. Properties of different schemes for matrix construction.

4. Simulation results

In this work, the GA-STD case (core) parameters are adopted as reported in the previous benchmark work (Görler et al. 2016). Concentric circular magnetic flux surfaces are adopted. The equilibrium density and temperature profiles, denoted as $H(\tilde r)$ , and the normalised logarithmic gradients $R_0/L_{\mathrm {H}}$ are given by

(4.1) \begin{eqnarray} \frac {H(\tilde r)}{H(\tilde r_0)} =\exp \left [{-}\kappa _{\mathrm {H}}w_{\mathrm {H}}\frac {a}{R_0}\tanh \left (\frac {\tilde r-\tilde r_0}{w_{\mathrm {H}}a}\right )\right ], \end{eqnarray}

where $\tilde r=\sqrt {(R-R_{{axis}})^2+(Z-Z_{{axis}})^2}$ , $L_{\mathrm {H}}=-\left [\textrm {d}\ln H(\tilde r)/\textrm {d}\tilde r\right ]^{-1}$ is the characteristic length of profile $H(\tilde r)$ and $\tilde r_0=a/2$ . The ion-to-electron mass ratio is $m_{\textrm {i}}/m_{\textrm {e}}=1836$ and the deuterium ( $m_{\textrm {i}}/m_{\textrm {p}}=2$ ) is the only ion species where $m_{\textrm {p}}$ is the mass of a proton. The on-axis magnetic field $B_0=2\textrm {T}$ , $\rho ^*=\rho _{\textrm {i}}/a=1/180$ , $\rho _{\textrm {i}}=\sqrt {2T_{\textrm {i}}m_{\textrm {i}}}/(eB_0)$ , aspect ratio $\epsilon =a/R_0=0.36$ , $T_{\textrm {e}}/T_{\textrm {i}}=1$ , characteristic length of temperature and density profiles $R_0/L_{T_{\textrm {i}}}=-(\textrm {d}\,\textrm {ln }T_{\textrm {i}}/\textrm {d}\tilde r)^{-1}=6.96$ , $R_0/L_{T_{\textrm {e}}}=-(\textrm {d}\ln T_{\textrm {e}}/\textrm {d}\tilde r)^{-1}=6.96$ , $R_0/L_{n}=-(\textrm {d}\ln n/\textrm {d}\tilde r)^{-1}=2.23$ , and collision frequency $\nu _{{coll}}=0$ .

Figure 4. (left) Two-dimensional mode structure and (right) the radial structure of the ITG mode for $n=25$ .

4.1. Single toroidal harmonic simulations for benchmark

The 2-D1F field solver is used in the benchmark with the GENE results in the previous work (Görler et al. 2016). The simulation for the $n=25$ mode is carried out using $8\times 10^6$ electrons and $10^6$ ions. For electrons, the gyro-average is switched off while the 4-point gyro-average is adopted for ions. The mass ratio is $m_{\textrm {i}}/m_{\textrm {e}}=1836$ . The reference Larmor radius is $\rho _N=0.0033422$ m. The time step size is $\Delta t=10^{-4} R_N/v_N$ for $n=25$ . The 2-D mode structure is shown in figure 4.

The growth rate of the ITG mode for different values of the toroidal mode number $n$ is studied. The growth rate is measured using the total field energy. The total field energy is defined as

(4.2) \begin{eqnarray} E_\Phi =-C_{\textrm {P}}\int \textrm {d}V \frac {1}{G_{\textrm {P}}}\delta \bar \Phi \delta \bar N, \;\; \end{eqnarray}

where $G_{\textrm {P}}$ is defined in (2.38) and $E_{ {\Phi }}$ is an approximate value of the field energy in $\delta \Phi$ in the limit $|k_\|/k_\perp |\ll 1$ and $|\nabla _\perp \ln G_{\textrm {P}}|/|\nabla \ln \delta \bar \Phi |\ll 1$ . A reasonably good agreement is observed between the results from the TRIEMG-GKX code and the GENE code, as shown in figure 5. The time step size is small enough to stabilise the omega-H mode and is at least smaller than $10^{-3}$ of the ITG period, sufficient to resolve the ITG mode accurately. The simulations with $4\times 10^6$ and $8\times 10^6$ electron markers give similar growth rates as an indicator of the convergence concerning the marker number. Note that the quasi-neutrality equation in the long wavelength is adopted in this work, while in GENE, there is no truncation in the quasi-neutrality equation, which can lead to the discrepancy for high- $n$ modes as shown in figure 5. The discrepancy at $n\approx 32$ (namely, $k_y\rho _s\approx 0.5$ using the notation in the previous work by Görler et al. 2016) is expected.

4.2. Comparison of the field-aligned FEM solver and the 2-D1F solver

4.2.1. Simulation of the ITG mode

The 3-D piecewise field-aligned FEM solver is verified by comparing it with the 2-D1F solver using the economical CBC parameters. Since the grids are still arranged in a traditional pattern without any shift, the toroidal Fourier filter can be readily applied in the 3-D solver. The parameters are the same as the nominal ones except the mass ratio $m_{\textrm {i}}/m_{\textrm {e}}=100$ and the reference Larmor radius $\rho _N=0.01$ m. From our simulations using various values of $\rho _{\textrm {N}}=0.02, 0.01, 0.005$ , the computational cost to avoid a numerical crash depends on several physics parameters as follows:

(4.3) \begin{eqnarray} C_{ {comp}}\propto C_{\Delta t}C_{ {marker}} C_{m_{\textrm {e}}} \approx \rho _{\textrm {N}}^{-1} \rho _{\textrm {N}}^{-2} m_{\textrm {e}}^{-1/2}, \end{eqnarray}

where $C_{\Delta t}$ is due to the reduction of the time step size $\Delta t$ as $\rho _{\textrm {N}}$ decreases, $C_{ {marker}}$ is due to the increment of the marker number as $\rho _{\textrm {N}}$ decreases, and $C_{m_{\textrm {e}}}$ is due to the reduction of the time step size $\Delta t$ as $m_{\textrm {e}}$ decreases. As a result, larger $\rho _N$ is adopted to make the simulation less costly. For the 3-D field-aligned FEM solver, the toroidal Fourier filter is applied to simulate a single toroidal harmonic. The comparison is shown in figure 6. A good agreement is observed.

Figure 5. Growth rate of ITG mode for different values of the toroidal mode number $n$ .

Figure 6. Growth rate and frequency of ITG mode for different values of the toroidal mode number $n$ , using the 2-D1F scheme and the 3-D field-aligned FEM.

4.2.2. Simulations of the zonal flow residual

For collisionless zonal flow dynamics, the Rosenbluth–Hinton residual flow test is widely used to evaluate the validity and accuracy of gyrokinetic simulations. The flat $q$ profile, and the uniform density and temperature profiles are adopted. The dimensionless parameters are $\rho _N=0.01$ (i.e. $a/\rho ^*\approx 60$ ) and the mass ratio $m_{\textrm {i}}/m_{\textrm {e}}=100$ . The on-axis magnetic field and the major and minor radii are the same as the GA-STD case in the previous sections. The radial, poloidal and toroidal grid numbers are $(32,64,8)$ for the 3-D field-aligned solver, with $n=0$ toroidal Fourier filter. For the case with the traditional 2-D1F solver, the radial and the poloidal grid numbers are $(32,64)$ . The ion and electron marker numbers are $4\times 10^6\ \text{and}\ 10^6$ , respectively. The time step size is $\Delta t=0.01 R_N/v_N$ . As an initial condition, the perturbed density is introduced by setting the marker weight according to $\delta f/f_0=A_w \exp [{-}(r-r_{w,\textrm {c}})^2/\Delta _{w}^2]$ , where $A_w=10^{-8}$ , $r_{w,\textrm {c}}=0.5$ , $\Delta _{w}/a=0.1$ and $r=\rho /a$ . The whole torus is simulated with the wedge number $ L_{\phi, {mult}} =1$ .

The time evolution of the scalar potential simulated using the 3-D field-aligned FEM solver for $q=2.1$ is shown in the left frame of figure 7. A probe at $r_{\textrm {prob}}=0.5$ is used to measure the time evolution of the scalar potential $\delta \Phi _{ {prob}}$ . The scalar potential oscillates and reaches the residual level given by $R_{\textrm {ZF}}=1/(1+q^2\sqrt {r/R_0})$ (Rosenbluth & Hinton 1998). As drift kinetic electrons are included in our simulation model, the high-frequency oscillation (‘ $\omega _H$ mode’) exists (Lee 1987) in the presence of the non-zonal component with finite $k_\|$ , which is observed in the time evolution in the left frame of the figure (the small magnitude high-frequency oscillations). The dashed red line denotes the theoretical value of the Rosenbluth–Hinton (R-H) zonal flow residual. The zonal flow residual from the simulation agrees with the theoretical value at the end of the simulation. Using the 2-D1f and the new 3-D solver, the zonal flow residual is simulated with different safety factor values, as shown in the right frame of figure 7. The zonal flow residual is calculated by averaging $\delta \Phi _{{prob}}$ in $\sim 30 R_N/v_N$ at the end of each simulation. The right frame shows a good agreement between the simulation results and the theoretical values. However, the field-aligned FEM solver does not give any computational benefit in simulating the zonal flow residual as the memory consumption and the computational cost are both higher than with the 2-D1F solver. The field-aligned FEM solver is expected to be more useful in multi- $n$ simulations.

Figure 7. (Left) Time evolution of the scalar potential at the probe position simulated using the 3-D field-aligned FEM solver. (Right) Zonal flow residual for different safety factor values $q$ calculated using the 2-D1F solver and the 3-D field-aligned FEM solver compared with the theoretical result (red dashed line).

4.3. Simulations of the nonlinear evolution

The nonlinear simulation is performed using the 2-D1F and the 3-D solvers with a diagnostic that traces the energy flux. To make the simulation affordable and the $\rho _N$ value relatively small, we choose $\rho _N=0.0066844$ , namely $a/\rho _N=90$ . The mass ratio is $m_{\textrm {i}}/m_{\textrm {e}}=100$ . For the case of a 2-D1F solver, we simulate the most unstable toroidal harmonic, using $10^6$ ion markers and $2.5\times 10^5$ electron markers, with the time step size $\Delta t=0.002 R_N/v_N$ . For the case of a 3-D solver, we use $8\times 10^6$ ion markers and $10^6$ electron markers, with $\Delta t=0.005 R_N/v_N$ . The time evolution of the total field energy is shown in the left frame of figure 8. For this case, the $n=10$ toriodal harmonic is the most unstable mode with a growth rate $\gamma _{n=10}=0.248 v_N/R_N$ . The time evolution of the field energy is very similar during the exponentially growing stage at $t\lt 40 R_N/v_N$ between the multi- $n$ simulation with the 3-D field-aligned FEM solver and the single- $n$ simulation with the traditional 2-D1F solver. The nonlinear saturation level in the multi- $n$ simulation is higher than that in the single $n$ simulation as other $n$ components also contribute to the field energy. The energy flux is calculated according to $Q=\int \textrm {d}v^3 (mv^2/2)\delta f \dot {r}_1/|\nabla r|$ , where $\dot {r}_1$ is the perturbed radial velocity in (2.31), consistent with our previous studies (Hatzky et al. 2002; Lu et al. 2019 b ). The time evolution of the energy flux is shown in the right frame of figure. 8. The radial averaged value of the energy flux is calculated to improve the statistic reliability in $0.60895\lt \rho /a\lt 0.75105$ with $\rho =\sqrt {(R-R_{{axis}})^2+(Z-Z_{{axis}})^2}$ , where the magnitude of the energy flux is significant. The energy flux increases in the exponentially growing stage $t\lt 50 R_N/v_N$ and decreases in the nonlinear stage after $t\approx 100 R_N/v_N$ , as there is no particle or energy source in the simulation. Longer time scale nonlinear simulations with smaller $\rho _N$ and particle and heat sources merit more effort in the future.

Figure 8. Time evolution of (left) the total field energy and (right) the energy flux.

4.3.1. Multi-toroidal harmonic simulations with nominal $\rho _N/a$

The multi- $n$ nonlinear simulations are carried out without any Fourier filter. Buffer regions are applied in the inner and outer boundaries to minimise the noise near the simulation boundary. The nominal parameters are adopted except $m_{\textrm {i}}/m_{\textrm {e}}=100$ . The grid numbers along the radial, poloidal and toroidal (parallel) directions are $N_r=96, N_\theta =192, N_\phi =8$ , respectively. The time step size is $\Delta t=0.0025 R_N/v_N$ . Here, $32\times 10^6$ electrons and $4\times 10^6$ ions are used. The simulation is run on 16 nodes (AMD EPYC Genoa 9554) of the Viper supercomputer at MPCDF, with 128 CPU cores on each node, with a processor base frequency of $3.1$ GHz and a max turbo frequency of $3.75$ GHz. It takes $\sim 1.136$ hours to simulate one normalised time unit $t_{\textrm {N}}=R_{\textrm {N}}/v_{\textrm {N}}$ .

The statistic error is indicated by the unbiased estimator of the variance (Hatzky et al. 2019). Generally, the statistical error for $N_g$ markers is given as

(4.4) \begin{eqnarray} \varepsilon _E=\frac {\sigma }{\sqrt {N_g}}, \end{eqnarray}

where $\sigma$ is the standard deviation. The standard deviation of the total particle number is given by

(4.5) \begin{eqnarray} \sigma _n = \sqrt {\sum _s\frac {1}{N_{ {mk},s}-1} \left [ \sum _{p=1}^{N_g} w_{s,p}^2 -\frac {1}{N_{ {mk},s}} \left (\sum _{p=1}^{N_g} w_{s,p}\right )^2 +\varepsilon _{ {Var},n} \right ] }, \end{eqnarray}

where $\varepsilon _{ {Var},n}$ is the statistical error of the variance of the total perturbed density, which is a minor correction to the first term when the marker number is sufficiently large. The standard deviation of the total current is given by

(4.6) \begin{eqnarray} \sigma _j = \sqrt {\sum _s\frac {q_s^2}{N_{ {mk},s}-1} \left [ \sum _{p=1}^{N_g} (v_{\|,p}w_{s,p})^2 -\frac {1}{N_{ {mk},s}} \left (\sum _{p=1}^{N_g} v_{\|,p}w_{s,p}\right )^2+\varepsilon _{{Var},j} \right ] } \;\;, \end{eqnarray}

where $\varepsilon _{ {Var},j}$ is the statistical error of the variance of the total perturbed current. The errors of the total perturbed density and current are shown in figure 9. The standard deviation is reasonably low during the nonlinear stage, suggesting the quality of the simulation. The two cases with different numbers of markers give a similar evolution of $\sigma _n$ and $\sigma _j$ , which indicates that the noise level is decreasing according to (4.4). In addition, the error increases slowly, but does not change dramatically in the nonlinear stage, indicating that it does not significantly worsen.

Figure 9. Standard deviation of the total perturbed density $\sigma _{\delta n}$ (red line) and current $\sigma _{\delta j}$ (black dashed line) ( $32\times 10^6$ electrons and $4\times 10^6$ ions). A case with half the marker numbers ( $16\times 10^6$ electrons and $2\times 10^6$ ions) is also shown, where $\sigma _{\delta n}$ is given by the blue dashed line and $\sigma _{\delta j}$ is given by the dashed magenta line.

Figure 10. Time evolution of field energy for multiple $n$ simulation using 3-D field-aligned FEM for nominal parameters ( $\rho _N=0.0033422$ m) except $m_i/m_{\textrm {e}}=100$ . The growth rates $\gamma$ of the zonal component ( $n=0$ ) and the turbulent part ( $n\ne 0$ ) of the total field energy are calculated using the linear fit of the time evolution of the logarithmic total field energy during the stable linear stage indicated by the two dashed red lines.

Figure 11. The 2-D mode structure (left) in the linear stage at $t=30R_N/v_N$ and (right) in the nonlinear saturated stage at $t=50R_N/v_N$ .

Figure 12. (left) Mode structure $\delta \Phi (r_c,\theta, \phi )$ in the magnetic flux surface at $r=r_c$ and (right) the Fourier spectrum $\delta \Phi _{m,n}(r_c)$ . The linear and nonlinear structures are shown in the upper frame and lower frame, respectively. Note that the maximum toroidal mode number $n$ in the simulation is limited by the poloidal grid number with the given profile of the safety factor $q$ . However, for the Fourier decomposition of the mode structure, the field value $\delta \Phi$ is calculated in a denser grid in $(r,\theta, \phi )$ coordinates with 512 and 1024 grid points along $\phi$ and $\theta$ , respectively. Thus, the Fourier component $\delta \Phi _{m,n}$ is less reliable for $n\gt 35$ , but physically, its magnitude should be small due to the finite Larmor radius effect.

Figure 13. (left) Radial mode structure and (right) spectrum of the toroidal harmonics in the linear stage ( $t=30R_N/v_N$ ) and nonlinear stage ( $t=50R_N/v_N$ ).

The time evolution of the total field energy is shown in figure 10. The zonal part ( $n=0$ ) and the turbulent part ( $n\ne 0$ ) are separated, and the corresponding total field energy is calculated. The growth rate of the zonal part is close to twice that of the turbulent part, supporting that the zonal component is due to the beat-driven excitation (Chen et al. 2024 a , b). Compared with the previous simulation for the studies of the zonal flow generation due to single toroidal harmonics (Wang et al. 2022), we have included all the toroidal harmonics $0\leqslant n\leqslant 35$ and thus $\gamma _{ {ZF}}/\gamma _{ {turb}}$ can deviate slightly from $2$ , which can also be due to the error in the fitting procedure related to the finite duration of the exponentially growing linear stage in our case. Our purpose is mainly to demonstrate the capability of the piecewise field-aligned FEM for multi- $n$ simulations and more dedicated studies on the zonal flow generation will be our future work.

The turbulent part of the 2-D mode structure is also visualised. The most unstable mode appears in the linear stage and becomes dominant, as shown in the left frame of figure 11. Note that the exponentially growing stage is not very long since we start the simulation with the initial perturbed level of $\delta f/f_0\sim 10^{-3}$ for nonlinear cases and the 2-D mode structure is not a pure state of a single toroidal harmonic, but a state with several unstable toroidal harmonics. As a result, the 2-D mode structure is not as pure as the single $n$ linear result in figure 4. In the nonlinear stage, the most unstable mode reaches its saturation level, and the turbulence spreads in the radial direction, as shown in the right frame where the zonal component is extracted to illustrate the features of the turbulent part.

The mode structure in the magnetic flux surface at $r=r_{\mathrm {c}}=0.5a$ is visualised in the left frames of figure 12. The mode structures are aligned along the magnetic field lines. Note that although the parallel grid number is small ( $N_\phi =8$ ), the construction of the field in the toroidal direction is possible since the parallel mode structure is smooth. The Fourier components are calculated and the logarithmic amplitude is shown in the right frames. In the linear stage at $t=30R_N/v_N$ , the peak of the spectrum is at $n\approx 25$ , consistent with the linear results in figure 5. In the nonlinear stage at $t=50R_N/v_N$ , other toroidal harmonics also grow. Specifically, the low- $n$ harmonics have a significant amplitude, consistent with previous particle-in-cell gyrokinetic turbulence simulations (Wang et al. 2011).

The features of the nonlinear evolution are summarised in figure 13. The radial profile of the turbulence intensity is calculated by extracting the zonal component ( $n=m=0$ ) and integrating it along the toroidal and the poloidal directions, as shown in the left frame. In the linear stage ( $t=30R_N/v_N$ ), the mode is localised near $r=0.5a$ , as indicated by the red dashed line. In the nonlinear stage ( $t=50R_N/v_N$ ), the nonlinear spreading occurs and the radial structure of $\langle \delta \Phi ^2\rangle _\theta$ is broader than that of the linear structure, consistent with that in figure 11 and the previous simulations (Mishchenko et al. 2023). The spectrum of the toroidal harmonics $\delta \Phi _n\equiv \sum _m|\delta \Phi _{mn}|^2$ is calculated by summing up the poloidal harmonics (over $m$ ) for each $n$ . Consistent with the spectrum pattern in the right frame of figure 12, in the linear stage ( $t=30R_N/v_N$ ), the spectrum of the toroidal harmonics peaks at $n\approx 25$ , as shown in the right frame of figure 13. In the nonlinear stage ( $t=50R_N/v_N$ ), other toroidal harmonics, especially the low- $n$ harmonics, also grow to a certain magnitude as indicated by the red line. Note that the simulation time in our nonlinear simulations is limited to $50R_N/v_N$ merely to demonstrate the capability of the piecewise field-aligned FEM by simulating the early nonlinear evolution. Longer simulations of the turbulent time scale (Bottino et al. 2011) rely on the complete form of the gyrocentres’ equations of motion and the noise control schemes such as the coarse-graining algorithm (Chen & Parker 2007 a) and the weight smoothing operator (Sonnendrücker et al. 2015).

5. Conclusions

We have formulated and implemented a piecewise field-aligned finite element method in tokamak geometry with magnetic shear. On one hand, the computational grids are aligned in the traditional pattern and the centres of the basis functions are located at the grid vertices. On the other hand, the finite element basis functions are defined on piecewise field-aligned coordinates. The linear and nonlinear simulations demonstrate the features of this scheme. Our discussions and results are in the framework of the finite element method (FEM), consistent with the previous theoretical and numerical work (Connor, Hastie & Taylor 1979; Cheng, Chen & Chance 1985; Beer, Cowley & Hammett 1995; Scott 2001; Lu, Zonca & Cardinali 2012; Zonca & Chen 2014), but extending it to FEM. The good properties of this scheme are as follows.

1. (i) It is formulated in the framework of the finite element method and conservation properties are inherited since the partition of unity is satisfied.

2. (ii) Strong grid deformation is avoided due to the piecewise treatment in defining the finite element basis functions.

3. (iii) High efficiency is expected in multi- $n$ nonlinear simulations due to the reduced grid number in the parallel direction.

4. (iv) It applies to particle (Lagrangian) and Eulerian schemes using the finite element method.

5. (v) The separation of the parallel direction from the two perpendicular directions enables the optimised treatment in the parallel direction for solving the ideal Ohm’s law (Hatzky et al. 2019).

6. (vi) The combination with a toroidal Fourier filter, partial torus simulations and the application with open field lines is possible.

The piecewise field-aligned FEM has been implemented in TRIMEG-GKX code for the electrostatic particle simulations with drift kinetic electrons and gyrokinetic ions. The Cyclone base case benchmark shows reasonable agreement between the TRIMEG-GKX results and the GENE results regarding the simplification in TRIMEG-GKX. The single toroidal harmonic simulations show agreement between the traditional 2-D1F solver (FEM in the poloidal plane and the particle-in-Fourier in the toroidal direction) and the 3-D field-aligned FEM solver. The multi- $n$ simulations demonstrate the capability of the field-aligned FEM in nonlinear turbulence studies. The nonlinear evolution of the ITG turbulence is simulated. The radial intensity profile and the spectrum of the Fourier modes demonstrate the nonlinear spreading of the ITG turbulence in real space and spectrum space.

As further applications, this scheme can be applied to the gyrokinetic simulations using unstructured meshes (Lu et al. 2019 a , 2024) for electromagnetic particle simulations (Hatzky et al. 2019; Lanti et al. 2020; Lu et al. 2023; Mishchenko et al. 2023), for the whole device modelling. It can also be useful for particle simulations in stellarators (Kleiber et al. 2024), especially for multi- $n$ nonlinear simulations.