Verification of the global gyrokinetic stellarator code XGC-S for linear ion temperature gradient driven modes

TL;DR

This paper benchmarks the extended XGC-S gyrokinetic code for stellarator configurations, demonstrating its accuracy in modeling ion temperature gradient-driven microinstabilities and showing good agreement with established codes in tokamak and stellarator geometries.

Contribution

It introduces and validates XGC-S for stellarator modeling, extending the original XGC code to 3D configurations and benchmarking its performance against other gyrokinetic codes.

Findings

XGC-S shows good agreement with XGC1, ORB5, and GENE in tokamak geometry.

XGC-S accurately reproduces mode spectra and growth rates in Wendelstein 7-X stellarator.

Benchmark results confirm XGC-S's reliability for stellarator turbulence studies.

Abstract

XGC (X-point Gyrokinetic Code) is a whole-volume, total-f gyrokinetic particle-in-cell code developed for modelling tokamaks. In recent work, XGC has been extended to model more general 3D toroidal magnetic configurations, such as stellarators. These improvements have resulted in the XGC-S version. In this paper, XGC-S is benchmarked for linear electrostatic ion temperature gradient-driven microinstabilities, which can underlie turbulent transport in stellarators. An initial benchmark of XGC-S in tokamak geometry shows good agreement with the XGC1, ORB5, and global GENE codes. A benchmark between XGC-S and the EUTERPE global gyrokinetic code for stellarators has also been performed, this time in geometry of the optimised stellarator Wendelstein 7-X. Good agreement has been found for the mode number spectrum, mode structure, and growth rate.

Equations

$$\frac{\partial \delta f}{\partial t}+\dot{\vec{X}}\cdot \frac{\partial \delta f}{\partial \vec{X}}+{\dot{v}}_{\Vert }\frac{\partial \delta f}{\partial v_{\Vert }}=-{\dot{\vec{X}}}_1\cdot \frac{\partial f_0}{\partial \vec{X}}-{\dot{v}}_{\Vert 1}\frac{\partial f_0}{\partial v_{\Vert }}$$

$$\dot{\vec{X}}=\frac{1}{G}\left[v_{\Vert }\vec{b}+\frac{m{v}_{\Vert }^2}{qB}\nabla \times \vec{b}+\frac{1}{q{B}^2}\vec{B}\times \left(\mu \nabla B+q\nabla ⟨\varphi ⟩\right)\right],$$

$${\dot{v}}_{\Vert }=-\frac{1}{mG}\left(\vec{b}+\frac{m{v}_{\Vert }}{qB}\nabla \times \vec{b}\right)\cdot \left(\mu \nabla B+q\nabla ⟨\varphi ⟩\right)$$

$$G=1+\frac{m{v}_{\Vert }}{qB}\vec{b}\cdot \left(\nabla \times \vec{b}\right),$$

$$-{\nabla }_{\perp }\cdot \frac{{\rho }_i^2}{{\lambda }_{Di}^2}{\nabla }_{\perp }\varphi =e\left(1-{\nabla }_{\perp }\cdot {\rho }_i^2{\nabla }_{\perp }\right)\left(⟨{\tilde{n}}_i⟩-{\tilde{n}}_e\right),$$

$$\frac{\mathrm{d}\ln T}{\mathrm{d}s}=-\sqrt{2}\left(\frac{1}{2}-\mid s-\frac{1}{2}\mid \right)\frac{a}{L_T},$$

Full text

I Introduction

Anomalous transport, believed to be principally turbulent, is known to be the main cause of heat and particle loss in high performance tokamaks. Anomalous transport has also been observed in stellarator experiments, including Wendelstein 7-AS Stroth (1998), Large Helical Device (LHD) Motojima et al. (2003) and Helically Symmetric eXperiment Deng et al. (2015). In the second operations phase of the Wendelstein 7-X (W7-X) stellarator Nührenberg and Zille (1986); Wolf et al. (2017) (OP 1.2), turbulence has been observed to be the dominant driver of transport Klinger et al. (2019). Theoretical work has shown that turbulent transport in stellarators can be qualitatively different to that in tokamaks Proll et al. (2012); Xanthopoulos et al. (2014); Plunk et al. (2019). Understanding turbulent transport in stellarators is therefore increasingly relevant.

Numerical tools are useful for investigating turbulent transport. One of the most effective physical models for numerical modelling of plasma turbulence has been 5D gyrokinetics, which combines high physical completeness with reduced computational resource requirements compared to a 6D kinetic model. Gyrokinetics is an approximation to the Vlasov-Poisson (or Vlasov-Maxwell) system that is valid for phenomena of low frequency relative to the ion gyrofrequency Brizard and Hahm (2007). The fast gyromotion is decoupled from the slower guiding centre motion. The 6D Vlasov equation is thus reduced to a 5D gyrokinetic equation.

In stellarator geometry, gyrokinetic simulation of ion-scale microinstabilities has been carried out with local models Proll, Xanthopoulos, and Helander (2013), flux-surface global models Xanthopoulos et al. (2014, 2016), and models fully global to the last closed flux surface Helander et al. (2015); Riemann, Kleiber, and Borchardt (2016). In the former two categories, nonlinear simulations have been performed to calculate saturated turbulent states and fluxes. Microinstability calculations global to the last closed flux surface have only been demonstrated with the EUTERPE code, linearly with a delta-$f$ model Kornilov et al. (2004); Helander et al. (2015); Riemann, Kleiber, and Borchardt (2016). Neoclassical transport has been modelled with a global full-$f$ gyrokinetic model to the last closed flux surface with the GT5D code Matsuoka, Idomura, and Satake (2018). The GTC Spong et al. (2017) and MEGA Todo et al. (2017) codes have also been used to perform global simulations to the last closed flux surface in LHD geometry. In these codes, an energetic particle species has been treated with a gyrokinetic or drift kinetic model, interacting with a background plasma red treated as an MHD fluid in MEGA and with a gyrokinetic ion/fluid electron hybrid model in GTC.

The gyrokinetic total-$f$ code XGC Ku, Chang, and Diamond (2009); Ku et al. (2015, 2018) is global to the first wall, and has so far been applied to tokamaks including the pedestal and scape-off layer in the simulation volume Churchill et al. (2017); Chang et al. (2017). Recently, XGC has been generalized to 3-dimensional geometry for application to stellarator plasmas Moritaka and et al. (pted); Cole et al. (2019). In the present paper, XGC-S is verified for linear Ion Temperature Gradient-driven (ITG) microinstabilities with a delta-$f$ model global to the last closed flux surface in axisymmetric and stellarator geometries. In Section II, the model is described and simulation results are presented for global ITG calculations with a standard circular tokamak case. This case has been used for benchmarking between the tokamak XGC1 code and other gyrokinetic codes. XGC-S correctly reproduces results of this benchmark. In Section III, the analysis is extended to linear global ITG instabilities in the Wendelstein 7-X (W7-X) stellarator Nührenberg and Zille (1986); Wolf et al. (2017), in comparison with the EUTERPE code Helander et al. (2015); Riemann, Kleiber, and Borchardt (2016). We conclude in Section IV.

II Model verification

II.1 Model description

XGC is a whole volume global total-$f$ Ku et al. (2015, 2018) gyrokinetic Particle-in-Cell (PIC) code. Its nonlinear tokamak version for neoclassical and turbulent physics is XGC1. It discretises the distribution functions of plasma species by numerical marker particles advanced in cylindrical ($R$, $\varphi$, $Z$) coordinates according to the gyrocenter equations of motion. Charge and current densities are accumulated on nodes of unstructured meshes. Potentials are calculated on mesh nodes, and fields are then calculated by taking spatial derivatives of the potentials. These are then interpolated for the marker particle equations of motion. A fully nonlinear Fokker-Planck Landau operator Hager et al. (2016) is used to model Coulomb collisions and a Monte-Carlo method is used for neutral recycling and transport with atomic interactions Stotler et al. (2017).

XGC-S does not yet include all of these features. In this paper, the conventional collisionless delta-$f$ method is used in electrostatic mode, with a single thermal ion species and the adiabatic approxmation for electrons. In the delta-$f$ formulation, $f=f_{0}+\delta f$, where $f$ is the total distribution function, $f_{0}$ is a prescribed background distribution function, and $\delta f$ is a small perturbation to this background represented by the markers. In this paper the background distribution function is always Maxwellian.

Neglecting collisions and any background electric field, $\phi_{0}$, the perturbed distribution function evolution equation Brizard and Hahm (2007) solved by XGC-S is

[TABLE]

where $\vec{X}$ is the gyrocentre position,

[TABLE]

$v_{\|}$ is the gyrocentre velocity,

[TABLE]

and

[TABLE]

with $B$ being the local magnitude of the magnetic field, $q$ the ion charge, and $\mu$ the magnetic moment. $\langle\rangle$ represents the gyroaverage operation. The corresponding $\dot{\vec{X}}_{0}$ and $\dot{v}_{\|0}$ are as above but excluding terms containing $\phi$, while the corresponding $\dot{\vec{X}}_{1}$ and $\dot{v}_{\|1}$ contain only these terms.

The EUTERPE simulations with which XGC-S is compared in this paper have been run excluding nonlinear terms from the equations of motion. That is, $\dot{\vec{X}}_{0}$ and $\dot{v}_{\|0}$, without electrostatic perturbations, are used in place of $\dot{\vec{X}}$ and $\dot{v}_{\|}$ when acting on the derivatives of the perturbed distribution function. Since the initial field perturbations used with XGC are small, and the simulations are run for linear physics while the perturbations are small, the two models are very closely approximately equivalent in the simulations performed for this paper.

To solve these equations, the background magnetic field, $\vec{B}$, must be known at all positions. The components of the background magnetic field, and the flux label $s=\Psi/\Psi_{LCFS}$ where $\Psi$ is the poloidal flux, are obtained from the VMEC code. VMEC is an ideal MHD equilibrium code which assumes nested flux surfaces Hirshman and Whitson (1983). The VMEC output is mapped to cylindrical coordinates for use with XGC. To simulate linear ITG modes, the perturbed electrostatic potential must be calculated. The gyrokinetic Poisson equation in the Padé approximation is

[TABLE]

where $\rho_{i}$ is the ion gyroradius, $\lambda_{Di}$ is the ion Debye length, and $\tilde{n}_{s}$ are the perturbed species densities. In this paper, the perturbed electron density $\tilde{n}_{e}$ is obtained by use of the adiabatic electron approximation. The Padé approximation is used in this section, while the short wavelength term $\nabla_{\perp}\cdot\rho^{2}_{i}\nabla_{\perp}$ is set to zero (long wavelength approximation) is used for Wendelstein 7-X simulations in the following section to permit comparison with published results of the EUTERPE code, obtained with the long wavelength approximation.

In XGC1, the Poisson equation is solved on a series of identical unstructured meshes composed of nodes placed on flux surfaces and mapped to positions on the cylindrical coordinate grid. Nodes are placed such that, for the greatest practical proportion of each mesh, tracing a field line from a node on one plane to a neighbouring plane will arrive close to another node. In XGC-S, unstructured meshes must have varying geometry at different toroidal locations. A series of unstructured meshes is therefore produced based on VMEC output Moritaka and et al. (pted) with nodes lying equidistantly on surfaces in PEST poloidal angle, and following field lines with toroidal angle. In Figure 1, stylised (reduced radial and poloidal resolution) meshes produced for XGC-S are plotted at different toroidal angles to illustrate the twisting of nodes in toroidal angle. Similar stylised meshes for Wendelstein 7-X geometry are shown in Figure 2, illustrating also the change in mesh geometry with toroidal variation of the plasma boundary.

The use of unstructured meshes is a key distinguishing feature of the XGC code from the EUTERPE code, and other major global gyrokinetic codes. It facilitates extension of the simulation domain into the edge region, where the magnetic field can be irregular.

The Poisson equation is currently solved independently on each planar mesh, with the electric field calculated in 3D using a finite differences method along the equilibrium magnetic field lines between the planes. This is a valid approximation so long as the phenomena of interest are highly elongated along the magnetic field lines.

II.2 Tokamak electrostatic ITG verification

The 3D equilibria and mapping techniques for the stellarator code were previously verified Cole et al. (2019). In that work, particle tracing studies were performed and collisionless orbit loss fractions compared when calculated with the EUTERPE, BEAMS3D McMillan and Lazerson (2014) and XGC codes. They showed good agreement, and this capability was then used to investigate collisionless particle confinement in proposed stellarator reactors Moritaka and et al. (pted); Cole et al. (2019).

As an initial verification of the extensions to the stellarator code for solving the Poisson equation for ion-scale microinstabilities, we compare to a previously published, well-benchmarked calculation. This calculation has been performed using the tokamak codes XGC1, ORB5 and global GENE, finding good agreement Merlo et al. (2018). In this case, a circular tokamak geometry is chosen (Case V Burckel et al. (2010)). An electrostatic delta-$f$ gyrokinetic calculation to the last closed flux surface is performed, with analytic background profiles that localise a strong temperature gradient at the central flux surface. The analytical form of the background plasma profiles is given in detail in reference Merlo et al. (2018), Section III.

Although the tokamak XGC1 and new XGC-S versions should, in this case, produce the same result, there are significant implementation differences. XGC-S, for instance, uses a series of non-axisymmetric unstructured meshes to solve the Poisson equation, where XGC-1 solves the Poisson equation on the same axisymmetric mesh at each toroidal location Moritaka and et al. (pted). Another is the interpolation of the equilibrium magnetic field, which in the extended stellarator case is performed using a 3D rather than 2D spline interpolation. This equilibrium magnetic field model has been verified previously Moritaka and et al. (pted); Cole et al. (2019).

In Figure 3, the calculated mode structure and growth rate are compared for this case using the two code versions. In the stellarator version XGC-S the calculated growth rate is $\gamma_{XGC-S}=55.63~{}\mathrm{kHz}$, while in the tokamak XGC1 version the calculated growth rate is $\gamma_{XGC1}=55.75~{}\mathrm{kHz}$, giving a discrepancy of $0.22\%$. Perfect agreement is not expected due to the difference in meshing; in particular, the region in which nodes are not field aligned (unavoidable on irrational surfaces) is necessarily not the same due to the difference in mesh design. Nonetheless, the discrepancy is remarkably small, implying that the difference in meshing may not be important in calculating the ITG growth rate. In both cases, the dominant toroidal mode number $n=24$ is chosen by simulating a $1/24$th toroidal wedge of the tokamak, and the dominant poloidal mode is correctly calculated to be $m=33$, as can be predicted considering the local value of the safety factor, $q=1.375$. The number of toroidal slices (poloidal planes) within the wedge is $8$, as in the previous work, and the number of marker particles used is comparable within an order of magnitude ($10^{7}$).

III ITG microinstabilities in Wendelstein 7-X

The only published result for simulation of ion-scale microinstabilities in stellarator geometry with a radial global model was obtained with the EUTERPE code, in comparison with flux tube ensemble (surface global) GENE simulations Helander et al. (2015); Riemann, Kleiber, and Borchardt (2016). In this case, an artificial radially localised temperature gradient was used to approximate the local limit. This case is chosen as a robust test for the model and implementation of XGC-S for linear electrostatic microinstabilities.

The magnetic geometry is provided by a fixed boundary VMEC run corresponding to the Wendelstein 7-X high mirror configuration. The temperature gradient chosen has the piecewise form

[TABLE]

where $T$ is the ion temperature, $s$ is the normalised flux, $a$ is the minor radius (here, $0.505~{}\mathrm{m}$) and $L_{T}$ is the temperature gradient scale length. A temperature gradient of this form localises the ITG drive around the central flux surface at $s=0.5$. This allows a numerically clean simulation that can easily be compared with other codes, including surface global but radially local codes. In order to maximise the ITG growth rate for this initial run, we choose not to apply a density gradient. A density gradient is typically stabilising. Instead, we choose $n_{i0}=n_{e0}=3\times 10^{19}~{}\mathrm{m^{-3}}$ at all locations, a typical value for low density operation in W7-X experiments.

With the EUTERPE code, this case has been run linearly. We therefore perform an initial set of linear runs with the XGC-S code. For the field calculation XGC-S is run with 256 toroidal domains per field period of W7-X, each mesh with approximately $1.4\times 10^{5}$ mesh nodes. As in the EUTERPE case, only one field period is simulated. Mesh nodes lie on flux surfaces, chosen to be equidistant in $s$ with a radial resolution $N_{s}=64$ locations, as in the EUTERPE simulations. The number of nodes per flux surface varies throughout a cross-section. At $s=0.5$, there are approximately $6\times 10^{3}$ nodes distributed equidistantly in PEST poloidal angle. In total, approximately $2.5\times 10^{9}$ numerical marker particles are used. As in the published EUTERPE case, the long wavelength approximation is used in solving for the electrostatic potential. The simulation domain radially is $0<s<0.95$.

In figure 4, the linear growth rate dependence on the temperature gradient strength parameter $a/L_{T}$ is compared between the XGC and EUTERPE codes Riemann, Kleiber, and Borchardt (2016). Good qualitative and quantitative agreement is seen, with discrepancy in trend lines of less than $5\%$. The mode structure in terms of the perturbed electrostatic potential is shown in Figure 5, as calculated by XGC. The potential is plotted at two poloidal cross-sections within the simulated field period of W7-X, at $\varphi=0$ and $\pi/5$, which correspond to the bean and triangle shaped cross-sections. This case has a temperature gradient strength factor $a/L_{T}=1.41$. Unlike in a tokamak, strong poloidal, and toroidal, localisation is expected and observed in the perturbation of fields and density. This is seen with XGC.

In Figure 6, the electrostatic potential structure is plotted in terms of $s$ and $\theta$ at zero toroidal angle (bean-shaped cross-section, left) and in terms of $\theta$ and $\phi$ at $s=0.5$ (right). In both cases, $\theta$ is the PEST poloidal angle. The left-hand plot can be qualitatively compared with Figure 4 in Ref. Riemann, Kleiber, and Borchardt, 2016. However, exact comparison is not possible because here the real electrostatic potential is plotted without phase factor transformation, whereas in the reference the absolute value of a complex electrostatic potential is plotted with phase factor transformation. The phase factor transformation obscures the fine scale structure by renormalising $n=0$ and $m=0$ to some finite value for the toroidal and poloidal mode numbers. Nonetheless, the key feature of a perturbation localised close to the $s=0.5$ surface around $\theta=0$ is reproduced.

In Figure 7, the power spectrum is plotted: $|\phi|^{2}$ in terms of the toroidal and poloidal mode numbers, $n$ and $m$. This can be compared to Figure 6 (left) in Ref. Riemann, Kleiber, and Borchardt, 2016. Note that in the EUTERPE case, the phase factor transformation and Fourier filtering mean that only components of $\phi$ with values of $n$ and $m$ close to the expected solution are retained in the simulation. In XGC-S, no phase factor transformation or filtering is used. One can see that, as expected, the mode follows the line of resonance around $n+\iota m=0$.

IV Conclusion

In this paper XGC-S, the stellarator version of the full volume global total-$f$ gyrokinetic code XGC, has been used to simulate, and verify the code for, linear ITG microinstabilities in a circular tokamak benchmark case, and a Wendelstein 7-X case. XGC-S extends the existing XGC1 code Ku, Chang, and Diamond (2009); Ku et al. (2015, 2018) by the implementation of an interface for magnetic equilibria from the VMEC code and the ability to generate and utilise toroidally varying unstructured meshes for solving for the electrostatic potential with a toroidally varying plasma boundary and magnetic field Moritaka and et al. (pted). This code has previously been benchmarked and applied for orbit tracing and neoclassical studies Cole et al. (2019); Moritaka and et al. (pted).

As an initial verification study, XGC-S has been used to reproduce a circular tokamak benchmark previously performed with the XGC1, ORB5 and global GENE codes Merlo et al. (2018). Close agreement has been found between the two versions for ITG mode structure and growth rate. The linear growth rates calculated by XGC1 and XGC-S differed by less than $1\%$.

Having established the validity of the code improvements by direct comparison to previous XGC results, the verification has been extended to stellarator geometry. Linear simulations of ITG microinstabilities in the Wendelstein 7-X have been performed with a localised ion temperature gradient and flat density gradient, in comparison with published results from the EUTERPE code Helander et al. (2015); Riemann, Kleiber, and Borchardt (2016). Linear growth rates and mode numbers have been directly compared, finding good agreement. The calculated mode structure also shows good qualitative agreement.

In on-going work with the XGC-S code, neoclassical physics is being investigated in the Large Helical Device (LHD) heliotron. In connection with this work the XGC-S physics model including field solver is being extended to the first wall Moritaka and et al. (pted) and for a total-$f$ capability. Nonlinear turbulence simulations are envisioned, opening up the possibility of comparison with experimental measurements. In the outlook, XGC-S is an appropriate tool for modelling of neoclassical and turbulent physics simultaneously, and studying their nonlocal interaction.

Acknowledgements.

We thank P. Helander, J. Geiger and A. Dinklage for useful discussions and assistance. We gratefully acknowledge S. Seal, M. Cianciosa, J. Geiger and S. Hirshman for their work developing, maintaining and running the VMEC code. This manuscript has been authored by Princeton University under Contract Number DE-AC02-09CH11466. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.