Development and testing of an unstructured mesh method for whole plasma gyrokinetic simulations in realistic tokamak geometry

TL;DR

This paper introduces a novel unstructured mesh finite element method for gyrokinetic simulations in realistic tokamak geometries, demonstrating significant speed improvements and successful applications to ITG mode modeling with experimental equilibria.

Contribution

The paper develops the TRIMEG code utilizing an unstructured mesh scheme and particle-triangle mapping, enabling efficient and accurate gyrokinetic simulations in complex tokamak geometries.

Findings

Speed-up factor of ~30 over brute force schemes.

Reasonable agreement with established codes like ORB5.

Successful simulations with realistic experimental equilibria.

Abstract

In this work, we have formulated and implemented a mixed unstructured mesh-based finite element (FE)-Fourier decomposition scheme for gyrokinetic simulations in realistic tokamak geometry. An efficient particle positioning (particle-triangle mapping) scheme for the charge deposition and field scattering using an intermediate grid as the search index for triangles has been implemented and a significant speed-up by a factor of $\sim30$ is observed as compared with the brute force scheme for a medium-size simulation. The TRIMEG (TRIangular MEsh based Gyrokinetic) code has been developed. As an application, the ion temperature gradient (ITG) mode is simulated using the simplified gyrokinetic Vlasov-Poisson model. Our simulation and that using the ORB5 code for the DIII-D Cyclone case show reasonable agreement. As an additional application, ITG simulations using an ASDEX Upgrade equilibrium have been performed with density and temperature gradient profiles similar to the Cyclone case. Capabilities of the TRIMEG code for simulations with realistic experimental equilibria in the plasma core and in the whole plasma volume with open field lines are demonstrated.

Tables

Table 1. Table 1: Parametes for ITG with adiabatic electrons (same as those in Ref. 30 ), where R 0 = 1.67 ​ m subscript 𝑅 0 1.67 𝑚 R_{0}=1.67m , ρ s = c s / Ω i subscript 𝜌 𝑠 subscript 𝑐 𝑠 subscript Ω 𝑖 \rho_{s}=c_{s}/\Omega_{i} , c s = T e / m i subscript 𝑐 𝑠 subscript 𝑇 𝑒 subscript 𝑚 𝑖 c_{s}=\sqrt{T_{e}/m_{i}} and the value of ρ s subscript 𝜌 𝑠 \rho_{s} is calculated using T e ​ ( r / a = 0.5 ) subscript 𝑇 𝑒 𝑟 𝑎 0.5 T_{e}(r/a=0.5) and B ​ ( r = 0 ) 𝐵 𝑟 0 B(r=0) .

$r_c/a$	$a/R_0$	$a/{\rho }_s$	$T_e/T_i$	$R_0/L_{ref}$	${\kappa }_{Ti}$	${\kappa }_n$	$W_n=W_{Ti}$
0.5	0.36	180	1	1	6.69	2.23	0.3

Table 2. Table 2: Parameters for the three ITG simulations using the AUG equilibrium. OFL refers to the Open Field Line region. For all cases, κ n = 2.23 subscript 𝜅 𝑛 2.23 \kappa_{n}=2.23 , κ T = 6.96 subscript 𝜅 𝑇 6.96 \kappa_{T}=6.96 .

Case	OFL	${\rho }_{p,c}$	$W_{p,c}$
A	w/o	0.5	0.3
B	w/	0.5	0.3
C	w/	1.0	0.1

Equations

$$\mathbf{B}=\nabla \psi \times \nabla \phi +F\nabla \phi ,$$

$$B_R=-\frac{1}{R}\frac{\partial \psi }{\partial Z},B_Z=+\frac{1}{R}\frac{\partial \psi }{\partial R},B_{\phi }=\frac{F}{R}.$$

$${\partial }_{\parallel }\delta \varphi =\left(b_R{\partial }_R+b_Z{\partial }_Z+\frac{b_{\varphi }}{R}{\partial }_{\varphi }\right)\delta \varphi ,$$

$${\frac{dR(\chi ,\xi ,l)}{d\varphi }|}_{\chi ,\xi }$$

$${\frac{dZ(\chi ,\xi ,l)}{d\varphi }|}_{\chi ,\xi }$$

$${\frac{dl}{d\varphi }|}_{\chi ,\xi }$$

$${\partial }_{\parallel }\delta \varphi ={\frac{\partial }{\partial l}|}_{\chi ,\xi }\delta \varphi .$$

$$\dot{\mathbf{R}}$$

$${\dot{\rho }}_{\parallel }$$

$$\dot{\mathbf{x}}$$

$${\dot{v}}_{\parallel }$$

$$\dot{\mathbf{R}}$$

$${\dot{v}}_{\parallel }$$

$$v_{d0}^R$$

$$v_{d0}^Z$$

$$v_{d0}^{\phi }$$

$${\mathbf{v}}_{d0}\approx \frac{v_{\perp }^2+2v_{\parallel }^2}{2{\omega }_{c}B}\mathbf{b}\times \nabla B,$$

$$v_{d0}^R\approx 0,v_{d0}^Z\approx -b_{\varphi }\frac{v_{\parallel }^2+v_{\perp }^2/2}{{\omega }_{c}B}{\partial }_R{B}_{\varphi }.$$

$$\delta v_E^R$$

$$\delta v_E^Z$$

$$\delta v_E^{\phi }$$

$$\delta {\mathbf{v}}_E=\frac{\mathbf{b}\times ⟨\nabla \delta \phi ⟩}{B}=\frac{1}{B}⟨[\hat{R}{\partial }_Z\delta \phi -\hat{Z}{\partial }_R\delta \phi ]⟩,$$

$${\dot{v}}_{\parallel 0}=-\frac{\mu }{B_{\parallel }^{\ast }}(B_R^{\ast }{\partial }_{R}B+B_Z^{\ast }{\partial }_{Z}B).$$

$${\dot{v}}_{\parallel 0}\approx -\mu (b_R{\partial }_{R}B+b_Z{\partial }_{Z}B).$$

$$\frac{d\delta f}{dt}=\tau (\mathbf{E}),$$

$$\tau (\mathbf{E})=-f_0\delta \mathbf{R}\cdot \nabla \ln f_0+f_0\frac{e}{m}⟨\mathbf{E}⟩\cdot {\frac{d\mathbf{R}}{dt}|}_0.$$

$$\delta f(\mathbf{R},v_{\parallel },\mu )$$

$$\delta f({\mathbf{R}}_p,v_{\parallel ,p},{\mu }_p){\Omega }_p=\frac{N_{ph}}{N}{w}_p(t),$$

$$\frac{d{w}_p}{dt}=\frac{N}{N_{ph}}{\Omega }_p\tau (\mathbf{E}).$$

$$f_{\text{mark}}=\frac{N}{N_{ph}}f,$$

$${\Omega }_p\equiv \frac{B_{\parallel }^{\ast }d\mathbf{R}d\mu d\alpha }{dN}=\frac{1}{f_{\text{mark}}},$$

$$\frac{d{w}_p}{dt}=-\delta \mathbf{R}\cdot \nabla \ln f_0+\frac{e}{m}⟨\mathbf{E}⟩\cdot {\frac{d\mathbf{R}}{dt}|}_0.$$

$$\frac{\delta n}{⟨n⟩}=\frac{V_{tot}}{N\Delta V}\sum _{p\in \Delta V}{w}_p,$$

$$\frac{\delta n_n(R,Z)}{⟨n⟩}=\frac{V_{tot}}{N\Delta V}\sum _{(R_p,Z_p)\in \Delta S}{w}_p{e}^{-in{\phi }_p},$$

$$-{\nabla }_{\perp }\frac{n_0}{{\omega }_{c}B}\cdot {\nabla }_{\perp }\delta \varphi =\delta n_i-\delta n_e.$$

Full text

I Introduction

Gyrokinetic simulations play an important role in predicting the transport level due to neoclassical physics or turbulence lee1983gyrokinetic ; lin1995gyrokinetic ; lin1998turbulent ; wang2015distinct ; ku2009full ; jenko2000electron ; dorland2000electron . One of the leading methods is the particle-in-cell (PIC) method. Numerous PIC codes, such as GTC lin1995gyrokinetic , GEM parker2000electromagnetic , ORB5 jolliet2007global , have been developed for simulations in the core of the tokamak plasmas. Edge plasma simulations have attracted significant attention in recent years due to their connection to, e.g., the high confinement regime of tokamak plasmas; the prediction of the divertor heat-flux width of ITER chang2017gyrokinetic ; edge localized mode (ELM) control huysmans2007mhd . In order to simulate the edge physics, besides comprehensive physics models qin2006general ; scott2006edge , numerical schemes such as finite element methods for unstructured meshes in XGC ku2009full ; ferraro2010ideal , GTS wang2015distinct and GTC/GTC-X nishimura2006finite ; bao2019global and multiple patches of structured meshes in JOREK huysmans2007mhd have been developed in order to treat the open field line (OFL) region. While whole plasma simulations for neoclassical transport, ELMs and micro-turbulence have been reported and various numerical schemes have been developed for treating the OFL geometry ku2009full ; chang2017gyrokinetic ; huysmans2007mhd ; nishimura2006finite ; de2019kinetic ; vugt2019kinetic , there is still space to understand the features of different schemes, such as the particle-in-Fourier method (cf. [​​19; 20] and references therein), and thus to optimize the efficiency and the fidelity of the whole volume simulation.

In this work, we developed the mixed unstructured mesh based finite element–Fourier decomposition scheme, i.e., the mixed particle-in-cell-particle-in-Fourier (PIC-PIF) scheme, for gyrokinetic simulations in general tokamak geometry, including the OFL region. In addition, an efficient particle positioning scheme for the charge deposition and the field gathering using an intermediate grid as the search index for triangles has been implemented. This work is organized as follows. In Section II, the physics model and numerical schemes are given. In Section III, we perform convergence/scaling studies and simulations of ion temperature gradient (ITG) modes using (I) the concentric circular magnetic geometry and the DIII-D cyclone parameters and (II) the AUG realistic magnetic equilibrium and analytical density and temperature profiles. In IV, conclusions and an outlook are given.

II Physics model and numerical methods

In the following, we will define the normalization in Section II.1. In Sections II.2–II.4, we will describe the three basic classes in the code, namely, the equilibrium, particle and field classes. We then describe the numerical methods in Section II.5. While the equations of motion for the guiding center, the weight equation, and the gyrokinetic Poisson equation (in the long wavelength limit) in this work are the same as or are the simplified version of other codes such as ORB5 jolliet2009gyrokinetic ; lanti2019orb5arxiv , the mixed unstructured mesh-based finite element-Fourier decomposition scheme makes our work different from ORB5. In ORB5, the OFL region is not included in either the particle pusher or the field solver. The mixed approach in this work thus also serves as a potential candidate for the extension of the present codes such as ORB5 towards whole plasma volume simulations. The details related to the finite element and unstructured meshes are described in another work for circular tokamak geometry lu2019mixed and will be omitted in this work.

II.1 Normalization

Normalization units are defined and physics quantities are normalized to the normalization units. The length unit is $R_{N}=1m$. The velocity unit is $v_{N}=v_{th,hy}$, where $v_{th,hy}=\sqrt{2T_{hy}/m_{hy}}$, $T_{hy}$ is the reference temperature, $m_{hy}$ is the the mass of hydrogen, the subscripts ‘$N$’ and ‘$hy$’ indicate “normalization” and “hydrogen” respectively. For each particle species s, $v_{th,s}$ is used as the velocity unit for the particle initialization according to a Maxwellian distribution function while unit conversion to the normalization unit $R_{N}$ is performed in the equations of motion for the guiding center and the field equation. The time unit is $R_{N}/v_{th,hy}$. The magnetic field unit is $B_{N}=1T$.

II.2 The coordinates and the equilibrium

In the right-handed coordinates $(R,\varphi,Z)$ and $(\psi,\varphi,\theta)$, where $\psi$ is the poloidal flux function, using the EFIT convention, the magnetic field is represented as

[TABLE]

where $F$ is the poloidal current function. In the $(\psi,\varphi,\theta)$ coordinates, the safety factor is defined as $q={{\bf B}\cdot\nabla\varphi}/({{\bf B}\cdot\nabla\theta})={JF}/{R^{2}}$, where $J=\{\nabla\psi\times\nabla\varphi\cdot\nabla\theta\}^{-1}$. The equilibrium variables are constructed using B-splines in the $(R,Z)$ plane of the $(R,\varphi,Z)$ coordinates. Equilibrium variables such as $B$, $B_{R}$, $B_{Z}$ and their derivatives in $R$ and $Z$ directions can be obtained using the B-spline subroutines.

The hybrid coordinates $(R,\varphi,Z)$ and $(\psi,\varphi,\theta)$ are used. On the one hand, the user specified computation grids are aligned along the magnetic flux surface using the $\psi$ coordinate, or, if the OFL region is included, along the plasma boundary; the refinement grids are generated when the Delaunay refinement algorithm is called for the generation of the unstructured meshes shewchuk2002delaunay . The refinement grids using the Delaunay algorithm are not necessarily along the magnetic flux surface but it is a widely used technique for the improvement of the mesh quality. Two cases of the grids and the unstructured meshes for the simulations in Section III are shown in Fig. 1. More details about the finite element for the unstructured meshes are in another work lu2019mixed . On the other hand, in order to treat the open field line geometry, the equations for the field and particles are solved in $(R,\varphi,Z)$ coordinates. In $(R,\varphi,Z)$ coordinates, the equilibrium magnetic field is expressed as

[TABLE]

TRIMEG relies on the B-Spline subroutines for the interpolation of $B_{R}$ and $B_{Z}$ accurately and keeps the error in $\nabla\cdot{\bf B}$ at a low level. For simulations in this manuscript, we use $N_{R}=N_{Z}=129$ as the grid numbers in R and Z directions for the equilibrium construction. In the relevant simulation domain, the error in $\nabla\cdot{\bf B}/B$ is well below $10^{-4}$ with the upper limit of $10^{-3}$. The parallel derivative is

[TABLE]

where $\delta\phi$ is the perturbed electrostatic scalar potential, $b_{R}=B_{R}/B$, $b_{Z}=B_{Z}/B$ and $b_{\phi}=B_{\phi}/B$.

The field-aligned coordinates are constructed as another option for the calculation of the parallel derivatives in addition to Eq. 3. Along the magnetic field line $\bf B$, the auxiliary Clebsch coordinates $(\chi,\xi,l)$ are determined by

[TABLE]

where $(\chi,\xi)$ labels a magnetic field line and can be taken as $(R,Z)$ at $\varphi=0$ for convenience and $l$ is the coordinate along $\bf B$. The parallel derivative using these Clebsch coordinates is

[TABLE]

For parallel derivatives of high $n$ (toroidal mode number) field-aligned modes, Eq. 7 can provide high accuracy compared with Eq. 3.

Since $(R,\varphi,Z)$ or $(\chi,\xi,l)$ coordinates are used for solving all equations in the code except when the grids are initialized according to $\psi$, and the safety factor defined in $(\psi,\varphi,\theta)$ does not appear explicitly but $B_{R}$, $B_{Z}$ and $B_{\varphi}$ are used directly in our equations, the singularity of the $q$ at the X point does not appear in any equation solved in TRIMEG. The auxiliary coordinates $(\chi,\xi,l)$ maintain the high accuracy of the parallel derivative calculation for high $n$ field-aligned modes and are similar to the flux-coordinate independent approach hariri2013flux ; stegmeir2016field . Discussions related to the flux coordinates and singularity of the safety factor at the X point can be found in [​​27] and references therein.

II.3 Particles

II.3.1 Equations of motion

In this work, we focus on electrostatic simulations. In order to describe the particle guiding center motion, we follow the canonical Hamiltonian equations white1990canonical . Using $({\bf{R}},\rho_{\parallel},\mu)$ as the guiding center coordinates, the equations of motion are as follows,

[TABLE]

where ${\bf B}_{\parallel}^{*}={\bf B}+\rho_{\parallel}\nabla\times{\bf B}$, $\rho_{\parallel}=mv_{\parallel}/(ZeB)$, $Z$ is the charge number, $H=Z^{2}e^{2}\rho_{\parallel}^{2}B^{2}/(2m)+m\mu B+Ze\delta\phi$ is the Hamiltonian, and $\mu=v_{\perp}^{2}/(2B)$. The equations of motion above are equivalent to those adopted in XGC ku2009full ,

[TABLE]

where $D\equiv 1+v_{\parallel}{\bf b}\cdot\nabla\times{\bf b}/B$ is related to the higher order corrections.

The variables ${\bf\dot{R}},\dot{v}_{\parallel}$ can be written as

[TABLE]

where the subscripts [math] and $E$ indicate the motion of the guiding center in equilibrium and that due to the electrostatic field. In $(R,\varphi,Z)$ coordinates, the contravariant components of the velocity $v^{\alpha}\equiv{\bf v}\cdot\nabla\alpha$ are calculated for different terms.

Magnetic drift ${\bf v}_{d0}$

[TABLE]

The dominant terms of the equations of motion are obtained by omitting the terms of the order of $\rho*=\rho_{th}/a$ and of $B_{\theta}/B$, where $\rho_{th}=mv_{th}/(ZeB)$, $v_{th}=\sqrt{2T/m}$. In $(R,\varphi,Z)$ coordinates, the dominant terms in the equations of motion are

[TABLE]

which, noticing that ${\bf b}\approx b_{\phi}\hat{\phi}$, where $\hat{\phi}$ is the unit vector in the toroidal direction, and $B\approx B_{\phi}$, can be further reduced to,

[TABLE] 2. 2.

$E\times B$ drift ${\bf v}_{E}$. Generally, $E\times B$ drift ${\bf v}_{E}$ is contributed from the equilibrium scalar potential $\delta\phi_{eq}$ and the perturbed scalar potential $\delta\phi$, i.e., ${\bf v}_{E}={\bf v}_{E,eq}+\delta{\bf v}_{E}$. In this work, we only consider the latter one, i.e.,

[TABLE]

where $\langle\ldots\rangle$ indicates a gyro average. In TRIMEG, the four point gyro average scheme is adopted. For the calculation of $\langle A(R,Z)\rangle$, where $A$ indicates $\delta\phi$ or its derivative in $R$ or $Z$ direction, $A(R\pm\rho_{\perp},Z)$ and $A(R,Z\pm\rho_{\perp})$ are calculated and the average value is obtained. The dominant term is

[TABLE]

where $\hat{R}$ and $\hat{Z}$ are the unit vectors in the $R$ and $Z$ directions. 3. 3.

Parallel acceleration $\dot{v}_{\parallel 0}$

[TABLE]

The dominant term is

[TABLE]

The simplification of other terms such as $v_{\parallel}$ and $\dot{v}_{\parallel E}$ is trivial and is omitted.

II.3.2 Weight equation

The gyrokinetic Vlasov equation for the evolution of the perturbed distribution function $\delta f({\bf R},v_{\parallel},\mu)$ is

[TABLE]

where

[TABLE]

The weight of the markers is defined to represent the perturbed distribution function,

[TABLE]

where $N_{ph}$ and $N$ are the total numbers of the physical particles and numerical markers respectively and the subscript $p$ is the marker index. Defining $\Omega_{p}$ as the phase space volume occupied by the marker $p$, Eq. 25 yields

[TABLE]

and Eq. 24 gives

[TABLE]

The above definition is the same as that in ORB5 jolliet2007global ; lanti2019orb5arxiv . However, in this work, for the sake of simplicity, we load the markers with the distribution function the same as that of the physical particles, i.e.,

[TABLE]

then

[TABLE]

where $\alpha$ is the gyro angle. Equations 27, 28 and 29 yield

[TABLE]

The perturbed density in a small volume $\Delta V$ is calculated from the marker weight in $\Delta V$

[TABLE]

where $\langle n\rangle$ is the volume averaged density. For unstructured meshes, the volume $\Delta V$ is centered around a vertex and $\Delta V=\Delta V_{\text{vert}}$ is calculated using $\Delta V_{\text{tria}}$, where $\Delta V_{\text{tria}}$ is a triangular prism which extends along the $\varphi$ direction. Using the particle-in-Fourier method in the toroidal direction, we have for each toroidal mode number $n$,

[TABLE]

where $\Delta S$ is the projection of $\Delta V$ in the $(R,Z)$ plane.

II.4 Field equation

The gyrokinetic Poisson equation with the long wavelength approximation is adopted in this work, i.e.,

[TABLE]

Generally, the electron response is dominated by the adiabatic response. Thus, the electron response can be decomposed into the adiabatic and non adiabatic (NA) parts, i.e.,

[TABLE]

where $\delta\tilde{\phi}$ is the non-zonal component, i.e., $\delta\tilde{\phi}=\delta{\phi}-\delta{\phi}_{0,0}$, $\delta\phi_{0,0}$ is the poloidal harmonic with $n=0,m=0$, where $m$ is the poloidal mode number. Notice that the Fourier decomposition is used in the $\varphi$ direction but the finite element method is used in the $(R,Z)$ plane.

For $n\neq 0$, with the subscript $n$ omitted,

[TABLE]

For $n=0$,

[TABLE]

In this work, we focus on the $n\neq 0$ modes while the studies involving $n=0$ components such as the geodesic acoustic mode will be reported in another separate work lu2019mixed .

II.5 Numerical methods

II.5.1 General description

This gyrokinetic Poisson-Vlasov system is implemented in Fortran. The field equation is solved using the finite element method for unstructured meshes. The sparse matrix corresponding to the gyrokinetic Poisson equation is solved using PETSc (Portable, Extensible Toolkit for Scientific Computation) balay2019petsc . The Runge-Kutta fourth order integrator is implemented for particles and coupled to the field solver. The Runge–Kutta fourth-order method is given by the following steps,

[TABLE]

where $d_{t}{\bf X}=d{\bf X}/dt$, ${\bf X}=({\bf R},v_{\parallel})$, $\Delta t$ is the time interval and $\delta\phi^{(i)}$ is solved from Eq. 33 with $\delta n(t)$ obtained using ${\bf X}_{t}^{(i)}$ for $i=1,2,3$.

II.5.2 Particle positioning (deposition/gathering) scheme

When calculating the toroidal component of the charge density perturbation $\delta n_{n}(R,Z)$ in Eq. 32 using marker weights $w_{p}$ in the so-called “charge deposition” stage, or when interpolating the field value at the particle position using the grid field value during the “field gathering” stage, the marker-triangle mapping, i.e., the particle positioning, needs to be treated. For the brute force particle position scheme, each triangle is checked for each marker whether the triangle contains the marker, which leads to a $NN_{t}$ scale computational cost, where $N$ and $N_{t}$ are the marker and triangle numbers respectively. In this work, rectangular grids (“boxes”) are constructed in $(R,Z)$ space and the box-triangle index $\{\text{Box}_{i}\mapsto\text{Triangle}_{j}\}$ is built when there is overlap between a box $i$ and a triangle $j$. The mapping $\{\text{Box}_{i}\mapsto\text{Triangle}_{j}\}$, $j=1,\ldots,N_{\text{i,tri}}$ is stored in the dynamically growing arrays for each box $i$. For a given marker $p$, the box which contains Marker $p$ is first found, i.e., the $\{p\mapsto\text{Box}_{i}\}$ mapping is identified. Then using the box-triangle mapping $\{\text{Box}_{i}\mapsto\text{Triangle}_{j}\}$, the corresponding triangle is identified. The computational cost is $\alpha N$ for $N$ markers, where $\alpha$ is a constant number.

The intermediate boxes are generated in the simulation domain with given $N_{x}$ and $N_{y}$, where $N_{x}$ and $N_{y}$ are the rectangular grid numbers in $R$ and $Z$ directions respectively. One limit is $N_{x}\ll N_{r}$, where $N_{r}$ is the radial grid number of the unstructured meshes. For $N_{x}=2$ (the box number is one), the positioning scheme is identical to the brute force scheme. The other limit is $N_{x}\gg N_{r}$, for which the box size is much smaller than the triangle size. A typical case between these two limit cases is that with $N_{x}\approx N_{r}$. These three cases are shown in Fig. 2. The computational speed-up versus the box size or $N_{x}$ will be studied in Sec. III.

III Numerical results

III.1 Parameters and simplifications for the simulation

In this section, two experimental cases are discussed. For numerical studies and benchmarks in Sections III.2 and III.3, the DIII-D Cyclone case is adopted, for which the parameters are the same as those in the benchmark work gorler2016intercode . The geometry with concentric circular magnetic flux surfaces is assumed. The nominal safety factor profile isgorler2016intercode

[TABLE]

where $\bar{r}\equiv r/a$. In Sections III.2 and III.3, an ad hoc equilibrium model is adopted. By assuming the following form of the safety factor profile,

[TABLE]

where $\bar{q}=q/\sqrt{1-r^{2}/R_{0}^{2}}$, the poloidal flux function can be obtained analytically bottino2007nonlinear ; lu2019mixed . The values of $q$ and the magnetic shear $\hat{s}$ are matched to Eq. 37 at $r_{c}=0.5a$, i.e., $q_{2}=q_{c}\hat{s}_{c}/(2r_{c}^{2})$, $q_{0}=(2/\hat{s}_{c}-1)\bar{r}_{c}^{2}q_{2}$, where $q_{c}$ and $\hat{s}_{c}$ are calculated at $r=r_{c}$ using Eq. 37. The temperature and density profiles indicated by $A(r)$ and the corresponding normalized logarithmic gradients indicated by $L_{ref}/L_{A}$ are given by

[TABLE]

where the subscript ‘c’ denotes the center of the gradient and the values of $r_{c}$, $W_{A}$ etc are in Table 1.

For the studies using the realistic geometry in Section III.4, the ASDEX Upgrade (AUG) case with shot number 34924 at 3.600s is chosen. This is a typical discharge for the study of energetic particle and turbulence physics lauber2018strongly . In the simulation, we use the experimental equilibrium but use the analytical density and temperature profiles in Eq. 40, with the radial coordinate replaced with $\rho_{p}=\sqrt{(\psi-\psi_{0})/(\psi_{b}-\psi_{0})}$, where $\psi_{0}$ and $\psi_{b}$ are the poloidal magnetic flux function at the magnetic axis and at the last closed surface respectively. The purpose of this study is to test the capability of treating the realistic geometry with minimum technical complexity. The fully self consistent treatment of the density/temperature profile and the equilibrium will be addressed in another work.

Since our purpose is to study the mixed PIC-PIF scheme and the particle search scheme in this work and address the basic ITG mode problem in the whole plasma geometry with minimum complexity, we have made the following simplifications.

Only the dominant terms in the equations of motion, Eqs. 17, 21 and 23, are solved. 2. 2.

The equilibrium variation of $n$, $B$ and $T$ in the gyrokinetic Poisson equation, Eq. 33, is ignored. 3. 3.

The ITG instability drive in Eq. 40 for the weight equation is kept but the equilibrium variation in $n$, $T$ and $B$ is omitted. 4. 4.

A single toroidal harmonic is simulated without the nonlinear terms, even though the dominant nonlinear term $\delta v_{E}\cdot\nabla\delta f$ for the ITG saturation is implemented in TRIMEG. 5. 5.

The Dirichlet boundary condition is adopted for the gyrokinetic Poisson equation with $\delta\phi=0$ at the boundary. The “absorbing boundary condition” for markers are adopted, i.e., the markers hitting the boundary are removed from the system. 6. 6.

Adiabatic electron approximation is adopted, i.e., $\delta n_{e}^{NA}=0$ in Eq. 34. 7. 7.

As the initial condition, markers with Maxwellian distribution are loaded in the simulation domain. Markers hitting the wall are removed (absorbing boundary condition). Since in this work, we only performed linear simulations, the marker distribution does not change after all absorbed markers are removed.

This simplified model can be replaced with a more comprehensive one by either future development of the TRIMEG code, or by implementing the finite element solver for unstructured meshes and the PIC-PIF scheme in other codes such as ORB5 and GTC.

III.2 Convergence and scaling studies

The effects of the rectangular box size on the computational cost in the particle positioning scheme are studied using a medium size case whose radial grid number is $N_{r}=90$ and the total marker number $N$ is 25.6 million. The brute force scheme ($N_{x}=2$) serves as the baseline and the speed-up for other values of $N_{x}=N_{y}=4,8,\ldots 4096,8192$ is shown in Fig. 3. There is an optimal value of the box size with respect to the triangle size in the range of $1<N_{x}/(2N_{r})<10$. The speed-up for $N_{x}=256,512,1024$ are $35.5,35.7,36.3$; larger than those for other values of $N_{x}$. For $N_{x}/(2N_{r})\ll 1$, each box contains a large number of triangles, as shown in Fig. 2 (left) and identifying the particle-triangle mapping consumes a significant amount of computing time. For $N_{x}=2$, the particle positioning in the charge deposition and field gathering can cost $>95\%$ of the total computing time. As $N_{x}/(2N_{r})$ increases and becomes larger than 1, the particle positioning consumption is reduced and the cost of the charge deposition and the field gathering is comparable to the particle pusher. For $N_{x}/(2N_{r})\gg 1$, the memory cost for storing the box-triangle mapping increases but without significant CPU cost, and thus only slows down the simulation slightly. Even as $N_{x}$ changes from 1024 to 8192, the speed-up decreases from 36.3 to 32.7, by only around $10\%$.

The convergence of the simulation results in terms of growth rate with respect to the radial grid number, the marker number per triangle and the time step size is studied. The convergence with respect to $N_{r}$, marker number $N$ and time step size $\Delta t$ is shown in Fig. 4. For this $n=20$ mode, as shown in the left frame, from $N_{r}=64$, the simulation starts to converge. The corresponding poloidal grid number per wave length $\approx 7$. Note that the poloidal grid spacing is set to be close to the radial grid spacing. For the linear studies in this work, the mode structure is elongated along the radial direction and the minimum value of $N_{r}$ is determined by the grid number per wave length in poloidal direction, i.e., $N_{\theta}/m\approx 2\pi N_{r}(r/a)/(nq)\gg 1$. In the middle frame, the results start to converge when the marker number per triangle $N/N_{t}>4$. In the right frame, the simulation starts to converge for $dt\leq 0.5$ and becomes numerically unstable for $dt>1$.

The parallel performance is tested for evaluating the scaling properties. The speed-up for the simulation with $N=25.6$ million markers with different numbers of cores is analyzed and shown in Fig. 5. Its comparison with the ideal scaling shows the good strong scaling for small to moderate core numbers (core number $\leq 1280$). For even larger core numbers, the consumption of the field solver parallel communication can increase since the field solver is distributed over all cores. As a result, the deviation of the speed-up curve away from the ideal scaling becomes significant as the core number is larger than 1280.

III.3 ITG simulation using Cyclone parameters

Following the convergence and scaling studies in Section III.2, ITG mode simulations using the Cyclone case parameters are performed. The growth rate and the frequency are shown in Fig. 6 and are compared with the ORB5 results noticing $v_{ti}/R_{N}=\sqrt{2T_{i}/T_{e}}(R_{0}/R_{N})(c_{s}/R_{0})$, where $R_{0}$ is the major radius. The agreement between the TRIMEG results and the ORB5 results is reasonable, noting the simplifications in TRIMEG as discussed in Sec. III.1. Nevertheless, considering the spatial scale separation between the equilibrium profile variation spatial scale $L_{E}$, the mode structure radial envelope width $L_{A}$ and the single poloidal harmonic width $1/(ndq/dr)$, i.e., $L_{E}\gg L_{A}\gg 1/(ndq/dr)$ lu2012theoretical , the simulation from TRIMEG already captures the leading order solution. More comprehensive physics models will be implemented in the future.

III.4 ITG simulation using AUG equilibrium in the core plasma and in the whole plasma geometry

In this section, we perform simulations of ITG modes using the AUG equilibrium described in Section III.1. The main purpose is to demonstrate the capability of treating the realistic magnetic equilibrium from an experiment using TRIMEG. Three cases are defined in Table 2. The simulations without and with the open field line region are performed and compared as shown in Fig. 7 and Fig. 8. For Cases (A) and (B), the open field line plays a weak role on the core ITG mode due to the narrow envelope of the radial mode structure. As a result, the 2D mode structures and the growth rate are almost identical between these two cases. For Case C, the local variables at $\psi=\psi_{b}$ such as $q$, $\hat{s}$, and the minor radius that determine the mode growth rate are different from (A) and (B). As a result, the growth rate curve is shifted due to the change in the finite Larmor radius / finite orbit width effects ($k_{\theta}\rho_{th}\approx nq\rho_{th}/r$) and other effects. While the convergence of the simulation for the whole plasma volume is achieved, a benchmark with other codes with the treatment of the whole plasma geometry will be studied in the future. Missing physics in the TRIMEG code includes, but is not limited to, the fully nonlinear collision operator hager2016fully ; a more realistic boundary condition such as a sheath boundary condition cohen2004sheath ; boesl2019gyrokinetic ; the radial electric field consistent with neoclassical physics; and zonal flow physics. In addition, more comprehensive gyrokinetic/gyrofluid models for the edge also need to be considered scott2006edge ; qin2006general .

IV Conclusion

In this work, the TRIMEG code has been developed based on the mixed unstructured mesh based FEM-Fourier decomposition scheme and the intermediate grid for the particle position. The parallel scalability of charge deposition and field gathering has been achieved and strong scaling up to moderate core numbers has been demonstrated. The benchmark with ORB5 using the DIII-D Cyclone test case shows reasonable agreement in terms of growth rate and frequency. The capability of treating the whole plasma volume is demonstrated by using an AUG magnetic equilibrium with an X-point and analytical density and temperature profiles. Futher development of the TRIMEG code for specific physics studies such as the mode structure symmetry breaking lu2017symmetry ; lu2018mode ; lu2018kinetic ; lu2019theoretical and for more comprehensive simulations including the wave-particle and wave-wave nonlinearities and multiple species will be explored in the future.

Acknowledgments

Simulations were performed on MPCDF computing systems. Suggestions from ORB5, EUTERPE, HMGC, GTC and GTS groups, discussion with B.D. Scott, D. Coster, A. Bierwage and support by National Natural Science Foundation of China under Grant No. 11605186 are appreciated by ZL. ZL thanks X. Wang for naming the code TRIMEG. Support by the ENR projects “NAT” and “MET” is acknowledged. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.