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 Reference Scott2001), 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 Reference McMillan2017).

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. Reference Hatzky, Kleiber, Könies, Mishchenko, Borchardt, Bottino and Sonnendrücker2019; Kleiber et al. Reference Kleiber2024; Lanti et al. Reference Lanti2020; Mishchenko et al. Reference Mishchenko, Borchardt, Hatzky, Kleiber, Könies, Nührenberg, Xanthopoulos, Roberg-Clark and Plunk2023, Reference Mishchenko2019), 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. Reference Jost, Tran, Cooper, Villard and Appert2001). 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}

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. Reference Lu, Lauber, Hayward-Schneider, Bottino and Hoelzl2019 Reference Lu, Lauber, Hayward-Schneider, Bottino and Hoelzla , 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.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 Reference Rosenbluth and Hinton1998). As drift kinetic electrons are included in our simulation model, the high-frequency oscillation (‘ $\omega _H$ mode’) exists (Lee Reference Lee1987) 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.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. Reference Hatzky, Kleiber, Könies, Mishchenko, Borchardt, Bottino and Sonnendrücker2019). 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. Reference Chen, Qiu and Zonca2024 Reference Chen, Qiu and Zoncaa , Reference Chen, Chen, Zonca and Qiub). Compared with the previous simulation for the studies of the zonal flow generation due to single toroidal harmonics (Wang et al. Reference Wang, Dai and Wang2022), 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. Reference Wang, Hahm, Ethier, Zakharov and Diamond2011).

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. Reference Mishchenko, Borchardt, Hatzky, Kleiber, Könies, Nührenberg, Xanthopoulos, Roberg-Clark and Plunk2023). 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. Reference Bottino, Vernay, Scott, Brunner, Hatzky, Jolliet, McMillan, Tran and Villard2011) rely on the complete form of the gyrocentres’ equations of motion and the noise control schemes such as the coarse-graining algorithm (Chen & Parker Reference Chen and Parker2007 Reference Chen and Parkera) and the weight smoothing operator (Sonnendrücker et al. Reference Sonnendrücker, Wacher, Hatzky and Kleiber2015).