Spectrally Accelerated Edge and Scrape-Off Layer Gyrokinetic Turbulence Simulations

Abstract

This paper presents the first gyrokinetic (GK) simulations of edge and scrape-off layer (SOL) turbulence accelerated by a velocity-space spectral approach in the full-$f$ GK code GENE-X. Building upon the original grid velocity-space discretization, we derive and implement a new spectral formulation and verify the numerical implementation using the method of manufactured solution. We conduct a series of spectral turbulence simulations focusing on the TCV-X21 reference case [Oliveira D. S. et al., Nucl. Fusion 62, 096001 (2022)] and compare these results with previously validated grid simulations [Ulbl P. et al., Phys. Plasmas 30, 107986 (2023)]. This shows that the spectral approach reproduces the outboard midplane (OMP) profiles (density, temperature, and radial electric field), dominated by trapped electron mode (TEM) turbulence, with excellent agreement and significantly lower velocity-space resolution. Thus, the spectral approach reduces the computational cost by at least an order of magnitude, achieving a speed-up of approximately $50$ for the TCV-X21 case. This enables high-fidelity GK simulations to be performed within a few days on modern CPU-based supercomputers for medium-sized devices and establishes GENE-X as a powerful tool for studying edge and SOL turbulence, moving towards reactor-relevant devices like ITER.

[ipp]organization=Max-Planck Institute for Plasma Physics, addressline=Boltzmannstr. 2, city=Garching, postcode=D-85748, country=Germany

[austin]organization=Institute for Fusion Studies, The University of Texas at Austin, city=Austin, postcode=TX 78712, country=USA

1 Introduction

Predicting turbulent transport in the edge and SOL regions is crucial to optimize fusion reactor performance, predict the divertor heat load, understand the L- to H-mode confinement transition, and design future magnetic confinement fusion devices, such as ITER [1] and DEMO [2]. Despite significant recent advancements in turbulent transport modeling [3], widely-used reduced turbulent transport (e.g., quasilinear) models and Braginskii-fluid simulations often fail to accurately capture edge and SOL turbulent transport due to, for instance, the importance of non-local effects and the fact that the edge is only marginally collisional. Therefore, because of the peculiar properties of the edge and SOL region, high-fidelity gyrokinetic (GK) turbulence simulations are necessary to overcome these difficulties and describe turbulent transport accurately.

While global [4, 5, 6] and local [7, 8, 9] GK codes for core turbulence are well-established, GK turbulence codes for the edge and SOL region remain less mature. One of the main reasons for this is (i) the lack of a clear separation between fluctuations and equilibrium quantities requiring a full-$f$ GK formalism and (ii) the complex magnetic geometry featuring open and closed field lines and X-points, which poses significant numerical challenges. In addressing this latter complexity, Braginskii-fluid turbulence codes such as GBS [10], TOKAM3X [11], and GRILLIX [12] have pioneered SOL turbulence, where the high-collisional assumption might be justified, in arbitrary magnetic configurations. Notably, the flux-coordinate independent (FCI) approach [13, 14], implemented in the GRILLIX fluid code [12, 15], has demonstrated promising performance in simulating turbulence in large devices [16], with flexible magnetic geometries [17].

Building on the experience gained from the GENE code [7, 4] and the flexibility of the FCI method, the GENE-X code [18] has been specifically designed to perform high-fidelity and high-performance GK turbulence simulations of the edge and SOL region with X-points. GENE-X is a full-$f$ GK code, i.e., it does not split the distribution function between an equilibrium and fluctuating parts. More precisely, GENE-X solves the full-$f$ electromagnetic and collisional GK Vlasov-Maxwell system and belongs to the continuum category of GK codes, where an Eulerian grid approach is utilized to discretize the velocity-space. Currently, GENE-X is one of the few GK codes able to perform edge and SOL turbulence with magnetic X-points. Among the other existing full-$f$ GK codes designed for edge and SOL applications, we can cite GKEYLL [19] and PICLS [20], which focus on either open and closed the open-field line region, and XGC [21, 22] and COGENT [23], which can also include magnetic X-points.

The GENE-X code has been validated against attached L-mode experiments in medium-sized devices such as in the ASDEX Upgrade [24] and TCV tokamak [25]. These validations demonstrate that GENE-X can provide predictions close to experimental measurements (e.g., OMP profiles, power balance, and divertor fall of length) in L-mode conditions. For instance, in the TCV-X21 reference case (L-mode discharge designed for code validation [26]), the GENE-X simulations have revealed the importance of the collisional cooling of trapped electrons in the edge to recover the correct electron temperature OMP profile within the experimental uncertainty, a kinetic mechanism that Braginskii-like fluid codes fail to capture due to the absence of trapped particles in these models. Despite these promising and encouraging results, the significant computational requirements of these first-principles grid GK simulations hinders the ability of GENE-X to simulate edge and SOL turbulence in reactor-relevant devices (such as ITER) and to explore high-performance and advanced experimental scenarios. Even for medium-sized devices, simulations frequently require several million CPU hours and span over several weeks to complete [24, 25]. Although GPUs and exascale HPC architectures offer potential, new numerical algorithms are needed to accelerate high-fidelity GK simulations.

This paper presents the first implementation of a spectral approach in velocity-space in a full-$f$ GK turbulence code such as GENE-X. The use of a spectral method is motivated by the fact that it can be particularly advantageous at high collisionality (e.g., in the SOL). The spectral method used in this work is based on a spectral expansion of the full-$f$ distribution function onto a Hermite and Laguerre polynomial basis in velocity-space. Using this basis, we derive and numerically implement the spectral formulation of the edge and SOL GK turbulence model solved in GENE-X. The numerical implementation is verified using the method of manufactured solution (MMS) [27]. We present the first spectrally accelerated GK edge and SOL turbulence simulations of the TCV-X21 reference case using GENE-X. It is noteworthy that the TCV-X21 scenario represents an ideal reference case for assessing the performance of the spectral approach, given that turbulence is dominated by TEMs [25], which can present a challenge for a global velocity-space spectral approach [28]. We compare our spectral results with the ones of Ref. [25] obtained from grid simulations with GENE-X. We find that the spectral approach can reproduce the OMP profiles in both the collisional and collisionless cases, dominated by TEMs, with excellent agreement and with a small number of spectral coefficients. In addition, by further increasing the spectral resolution in our simulations, we demonstrate that the agreement between the grid and spectral results is improved. Finally, the computational cost of the spectral simulations is assessed. We demonstrate that the spectral approach implemented in this work achieves a significant speed-up of CPU-based GK simulations with GENE-X. For the TCV-X21 reference case, a speed-up of nearly $50$ times is achieved compared to the previous grid simulations, allowing high-fidelity GK edge and SOL turbulence simulations to be conducted within a few days on current CPU-based supercomputers. This opens up new opportunities for studying edge and SOL turbulence through high-fidelity and high-performance GK simulations, which is crucial for the success of ITER and the design of future fusion power plants.

This paper is structured as follows. First, we introduce the GK turbulence model for the edge and SOL used in GENE-X in Section 2. We then present the spectral approach considered in this work and derive the spectral formulation of the GK model in Section 3. The numerical implementation of the spectral approach in GENE-X is detailed in Section 4, while its verification is carried out in Section 5. The first spectrally accelerated turbulence simulations of the TCV-X21 reference case are reported in Section 6, including a comparison with the grid simulations from Ref. [25]. Finally, the performance and computational cost of the spectral simulations are evaluated in Section 7. The conclusions and outlooks are presented in Section 8.

2 The GK turbulence model for the edge and SOL

GENE-X evolves the full-$f$ gyrocenter and gyroaveraged distribution function, $f_{\alpha}=f_{\alpha}(\bm{R},v_{\parallel},\mu,t)$ of particle species $\alpha$ in the gyrocenter phase-space described by the coordinates $(\bm{R},v_{\parallel},\mu,\theta)$. Here, $\bm{R}$ denotes the gyrocenter position related to the particle position $\bm{r}$ by the transformation $\bm{r}\simeq\bm{R}+\bm{\rho}_{\alpha}(\mu,\theta)$, with $\bm{\rho}_{\alpha}=\bm{v}_{\perp}(\mu,\theta)/\Omega_{\alpha}$ the particle Larmor radius and $\Omega_{\alpha}=q_{\alpha}B/(cm_{\alpha})$ the gyrofrequency, where $m_{\alpha}$ and $q_{\alpha}$ are species mass and charge, respectively. The velocity dependence of the gyrocenter phase-space is described by the parallel and perpendicular components of the gyrocenter velocity $\bm{v}$ relative to the equilibrium magnetic field, i.e., $v_{\parallel}=\bm{v}\cdot\bm{b}$ with background magnetic field unit vector $\bm{b}=\bm{B}/B$ and $\bm{v}_{\perp}=\bm{v}-\bm{b}v_{\parallel}$, respectively. Finally, $\mu=m_{\alpha}\norm{\bm{v}_{\perp}}^{2}/(2B)$ is the gyrocenter magnetic moment, and $\theta$ is the gyrocenter gyroangle, such that $\partial_{\theta}f_{\alpha}=0$. To obtain the self-consistent evolution of $f_{\alpha}$, the long-wavelength electromagnetic and collisional GK Vlasov equation is used, which is given in conservative form,

$$\frac{\partial\left(B_{\|}^{*}f_{\alpha}\right)}{\partial t}+\gradient\cdot\left(B_{\|}^{*}\dot{\bm{R}}f_{\alpha}\right)+\frac{\partial}{\partial v_{\|}}\left(B_{\|}^{*}\dot{v}_{\parallel}f_{\alpha}\right)=\sum_{\beta}B_{\|}^{*}C_{\alpha\beta}(f_{\alpha},f_{\beta}),$$ (1)

where the conservation of the magnetic moment, $\dot{\mu}=0$, is used. In Eq. (1), the equations of motion are

$$\dot{\bm{R}} =v_{\|}\frac{\bm{B}^{*}}{B_{\|}^{*}}+\frac{c}{q_{\alpha}B_{\|}^{*}}\bm{b}\times\left(\mu\gradient B+q_{\alpha}\gradient\psi_{1\alpha}\right),$$ (2a)

$$\dot{v}_{\parallel} =-\frac{\bm{B}^{*}}{m_{\alpha}B_{\|}^{*}}\cdot\left(\mu\gradient B+q_{\alpha}\gradient\psi_{1\alpha}\right)-\frac{q_{\alpha}}{m_{\alpha}c}\frac{\partial A_{1\parallel}}{\partial t},$$ (2b)

with $\bm{B}^{*}=\bm{B}_{T}+m_{\alpha}v_{\|}c\gradient\times\bm{b}/q_{\alpha}$. Here, $\bm{B}_{T}=\bm{B}+\gradient A_{1\|}\times\bm{b}$ is the sum of the constant (in time and space) magnetic field, $\bm{B}$, and the small perpendicular magnetic fluctuation, $\delta\bm{B}_{\perp}\simeq\gradient A_{1\|}\times\bm{b}$, with $A_{1\|}$ the parallel component of the perturbed magnetic vector potential. We remark that the gyrocenter Jacobian, $B_{\parallel}^{*}/m_{\alpha}$, is proportional to $B_{\parallel}^{*}=\bm{b}\cdot\bm{B}^{*}=B+m_{\alpha}v_{\|}c\bm{b}\cdot\gradient\times\bm{b}/q_{\alpha}$. The second term in $B_{\parallel}^{*}$ is referred to as the guiding-center correction and is proportional to the ratio of the particle Larmor radius $\rho_{\alpha}$ to the scale length of the equilibrium magnetic field, $L_{B}\sim R$ (with $R$ being the major radius of the fusion devices), i.e. it is a small correction proportional to $\rho_{\alpha}/L_{B}\ll 1$. It is worth noticing that the grid implementation of GENE-X (see Section 4) solves the GK Vlasov equation in the advection form rather than the conservative form given in Eq. (1). However, the conservative form is more convenient to derive the spectral approach in Section 3.

In Eq. (1), we introduce the generalized potential, $\psi_{1\alpha}$, defined by

$$\psi_{1\alpha}=\phi_{1}-\frac{m_{\alpha}c^{2}}{2q_{\alpha}B^{2}}\left|\gradient_{\perp}\phi_{1}\right|^{2},$$ (3)