Global electromagnetic turbulence simulations of W7-X-like plasmas with GENE-3D

The GENE-3D code, the global stellarator version of the established GENE framework, has been extended to an electromagnetic gyrokinetic code. This paper outlines the basic structure of the algorithm, highlighting the treatment of the electromagnetic terms. The numerical implementation is verified against the radially global GENE code in linear and nonlinear tokamak simulations, recovering excellent agreement between both codes. As a first application to stellarator plasmas, linear and nonlinear global simulations with kinetic electrons of ion temperature gradient (ITG) turbulence in Wendelstein 7-X were performed, showing a decrease of ITG activity through the introduction of electromagnetic effects via a finite plasma-$\beta$. The upgrade makes it possible to study a large variety of new physical scenarios, including kinetic electron and electromagnetic effects, reducing the gap between gyrokinetic models and physically realistic systems.


Introduction
Wendelstein 7-X (W7-X) is the first stellarator that has been optimised for low neoclassical transport (Wolf 2008), besides other criteria. As such, it has been shown (Klinger et al. 2019) that turbulence has become the limiting factor in the confinement for a broad range of its experiments. In order to understand the inherently nonlinear nature of plasma turbulence, numerical simulations based on gyrokinetic theory have become indispensable (Garbet et al. 2010).
Using some of the world's most powerful supercomputers, it is possible nowadays to simulate gyrokinetic turbulence in stellarators globally, making it possible to take into account the full variation of the magnetic field on a flux surface while simultaneously considering its radial variations as well as temperature and density profiles -all of which are inherently impossible in flux-tube simulations. As such, studies like those of Navarro et al. (2020), Cole et al. (2020), Sánchez et al. (2020) and Wang et al. (2020) have paved the way in investigating ion temperature gradient (ITG) turbulence globally in stellarators.
However, while these results already are of high importance, all of them relied on an adiabatic electron model, where only the ions were treated as a gyrokinetic particle species. This model fails to address things like trapped electron mode turbulence or the interaction between trapped electrons and ion turbulence, as well as electromagnetic effects, something that can become important for high-density operation of fusion devices.
All of these considerations motivated the upgrade of GENE-3D , the extension of the widely established GENE framework (Jenko et al. 2000) for three-dimensional magnetic field equilibria, to an electromagnetic gyrokinetic code, which is the focus of this work. In this paper, the underlying model of GENE-3D is presented, with special emphasis on the treatment of the electromagnetic terms in the gyrokinetic system of equations. Linear and nonlinear global verification studies using artificially heavy electrons to save computational cost against GENE are presented, showing excellent agreement between both codes, verifying the correctness and applicability of the electromagnetic upgrade. Finally, GENE-3D is used to perform global simulations of plasmas in W7-X with realistic electron mass and electromagnetic effects, something that, to the knowledge of the authors, has not been published before. It is found that, for the particular analytical background profiles used, the introduction of finite plasma-β leads to a reduction in linear as well as nonlinear ITG strength.
The rest of the paper is structured as follows. The GENE-3D framework is introduced in § 2, highlighting the gyrokinetic model it uses and the numerical methods to solve the corresponding set of equations. In § 3, results of verification studies of the code against the global version of GENE in linear and nonlinear, electromagnetic tokamak simulations are presented, testing the correctness of the numerical implementation. In § 4, the first global, nonlinear gyrokinetic simulations of W7-X incorporating kinetic electrons with realistic mass, as well as electromagnetic effects are discussed. Finally, the main results of this paper are summarised in § 5 and an outlook for future projects is given.

Extending GENE-3D to an electromagnetic code
This section focuses on the algorithmic details of the extension of GENE-3D to include effects of a perturbed parallel vector potential. First, a brief overview of the gyrokinetic system of equations is given in § 2.1, followed by details specific to the implementation in § 2.2. Finally, additional information about the numerics of the code is presented in Appendix A.

Gyrokinetic Vlasov-Maxwell system
The GENE-3D code is a Eulerian code that solves the gyrokinetic Vlasov equation, coupled to Maxwell's equations (e.g. Brizard & Hahm 2007) on a five-dimensional grid in phase space, using two velocity dimensions -the velocity parallel to the equilibrium magnetic field v || and the magnetic moment μ = mv 2 ⊥ /2B 0 -and three spatial dimensions (x, y, z), representing a field-aligned coordinate system (Beer, Cowley & Hammett 1995). Furthermore, it employs the so-called 'δf ' splitting (e.g. Garbet et al. (2010), and references therein). In this approach, one assumes that the full gyrocentre distribution function F σ of species σ can be split into a stationary background distribution function F 0,σ , describing the plasma at thermal equilibrium, and a time-dependent, first-order perturbation F 1,σ : where δ is a smallness parameter describing the scale relation between background and perturbed quantities. The GENE-3D code assumes that the background distribution is given by a local Maxwellian: Here, m σ , n 0,σ and T 0,σ are the mass, equilibrium background density and temperature of particle species σ at position x, respectively. Furthermore, B 0 gives the strength of the equilibrium magnetic field, v || and μ denote the velocity parallel to the magnetic field and the magnetic moment, respectively, and v th,σ = 2T 0,σ (x)/m σ is the thermal velocity of particle species σ . Under this assumption, the gyrokinetic equation for particle species σ , including a perturbation of the parallel vector potential, can be written as where χ = φ 1 − v || A 1,|| /c is called the gyrokinetic potential. The operator G{.} defines the forward-gyroaverage operator where X denotes the gyrocentre position and r is the gyroradius vector orthogonal to the local magnetic field. The drift velocities appearing in (2.3) are defined as follows: where Ω σ = (q σ B 0 )/(m σ c) is the gyrofrequency of species σ , b 0 = B 0 /B 0 is the unit vector in the direction of the equilibrium magnetic field B 0 , p 0 is the equilibrium thermodynamic pressure, β = 8πp 0 /B 2 0 is the ratio between plasma pressure and magnetic field strength and c is the speed of light.
At the current stage, GENE-3D uses a linearised Landau-Boltzmann collision operator for C [F 1,σ ]. Furthermore, the interplay between neoclassical and gyrokinetic physics is represented by the term coupling the curvature and ∇B drifts to the background Maxwellian (Oberparleiter et al. 2016), and the effects of a long-wavelength, radial electric field on the system (Helander & Simakov 2008;Mishchenko & Kleiber 2012) are represented by the term proportional to v E 0 in (2.3). Though they are incorporated in the main algorithm, these effects are neglected in the simulations presented in this paper for simplicity.
The system is closed by the gyrokinetic field equations. Poisson's equation, providing the electrostatic potential, reads where only the spatial derivatives perpendicular to the magnetic field are retained, which is in accordance with the gyrokinetic ordering. In order to calculate the density perturbation at the particle position x from the gyrocentre position X , GENE-3D employs a first-order pull-back transformation (e.g. Brizard & Hahm 2007): The operator K{.} denotes the backward-gyroaverage operator, which is defined here as Assuming a quasi-neutral limit (∇ 2 ⊥ φ 1 ≈ 0), Poisson's equation can be rewritten as At the current stage, magnetic compressibility, associated with a parallel perturbation of the magnetic field, is neglected, so that only the parallel vector potential A 1,|| is retained in GENE-3D, which is determined by the parallel component of Ampere's law: 2.2. Electromagnetic model Although (2.3), (2.9) and (2.10) are an analytically closed system, particular care has to be taken of the term containing the partial time derivative of A 1,|| in (2.3). The approach used in GENE for a long time was to introduce a new distribution function g σ = F 1,σ + (q σ v || F M,σ A 1,|| )/(T 0,σ c) and solving the corresponding gyrokinetic system for g σ instead of F 1,σ ). However, it was shown in Crandall (2019), that such a scheme becomes numerically unstable in nonlinear simulations at high plasma-β, although being stable linearly. In the same publication, an alternative scheme was proposed, similar to the one presented in Reynders (1993), which has been included in the standard version of GENE recently. This scheme is also used in GENE-3D and is presented in the following. By defining the parallel inductive electric field as where the term R σ contains all the contributions to the right-hand side of the gyrokinetic equation (2.3), except the one accounting for the parallel inductive electric field, and reads The next step is to take the partial time derivative of Ampere's law (2.10) and multiply by (−1/c): (2.14) Inserting definition (2.12) into (2.14) eliminates the explicit coupling of the equation to the particle distribution function F 1,σ . Therefore, one obtains an equation that treats the inductive electric field as an independent electromagnetic field: (2.15) Using (2.9), (2.10), (2.12), (2.13) and (2.15), the evolution of the distribution function F 1,σ can be calculated by treating time as a discrete variable and using a numerical scheme to solve (2.12) in the above system approximately. For this, GENE-3D provides several explicit Runge-Kutta schemes, with a fourth-order-accurate Runge-Kutta scheme (RK4) being the default option.

Code verification
Having introduced the major changes performed to the GENE-3D code, the following subsections focus on comparing its results of global, electromagnetic tokamak simulations with those obtained from GENE. As the implementation of the electromagnetic terms in the gyrokinetic equation is independent of the field geometry, these cases provide a simple scenario for code validation. Benchmarks with other stellarator codes capable of performing electromagnetic simulations are left for future work.

Linear simulations with varying plasma-β
This subsection presents a linear verification study between GENE-3D and GENE over a large range of plasma pressures. For this, linear growth rates and mode frequencies of the most unstable mode are calculated varying the plasma-β and are compared with the results presented in Görler et al. (2016).
The geometry is a tokamak with circular, concentric flux surfaces, with a magnetic field strength of B 0 = 2.0 T on axis, which is called B ref in the following, a major radius of R 0 = 1.67 m, an inverse aspect ratio of a/R 0 = 0.36 and a safety factor profile (3.1) Furthermore, density and temperature profiles are defined according to for both electrons and ions, where A = (n, T i , T e ). The corresponding normalised gradients are with T(x 0 ) = 2.14 keV, where x 0 /R 0 = 0.5 is the position at which the gradients have their maximum. The density at position x 0 is adjusted to obtain the desired value of the reference plasma β( The parameters κ A and w A set the amplitude and the width of the gradients, respectively, and are chosen to be κ T i = κ T e = 6.96, w T i = w T e = 0.3, κ n i = κ n e = 2.23 and w n i = w n e = 0.3. The box lengths are chosen to be (L and L y = 21.13 ρ s , which results in resolving multiples of the toroidal mode number n 0 = 19. Here, ρ s = c s /Ω s is the ion sound Larmor radius at reference position x 0 . The simulation is performed using deuterium as ionic species, as well as electrons that have twice their realistic mass. The resolution at which GENE-3D converged is (1344 × 16 × 16 × 64 × 32) in (x, y, z, v || , μ), respectively, and hyperdiffusion was set as η x = 2.0, η y = 0.05, η z = 2.0 and η v || = 0.2. As is mentioned in Appendix A, periodic boundary conditions are employed in the bi-normal direction, whereas zero-valued Dirichlet boundary conditions are used for the radial domain. In order to avoid numerical instabilities caused by the fixed boundaries, buffer zones with a Krook damping operator are used at the radial boundaries, each covering 5 % of the domain. Figure 1(a) shows the growth rate γ of the dominant linear instability as a function of β, whereas figure 1(b) shows the frequency ω of the respective mode. The red lines indicate the results obtained from GENE, whereas blue represents the results of GENE-3D. Excellent agreement is found for both frequencies and growth rates as the relative difference between the results of both codes is below 3 % in all cases. In particular, one can observe the expected stabilisation of the ITG mode by increasing β for both codes, as the growth rate decreases for β between 0.5 % and 1.4 %, while the frequency changes only slightly in the same interval. For values of β between 1.4 % and 1.75 % one can observe a transition of dominant modes from an ITG mode to a kinetic ballooning mode (KBM), which can be identified by the rapid increase in the mode frequency. The results for this interval are only shown here for the GENE code, since resolving this transition requires a much higher resolution than the one that was already used, which makes it impractical to investigate it with GENE-3D. Increasing β further results in an amplification of the KBM growth rate, while simultaneously causing a moderate decrease in its frequency, as observed with both codes. Finally, the radial and poloidal structures of the KBM at β = 2.5 % are compared in figure 2, again showing excellent agreement between both codes.

Nonlinear turbulence at finite plasma-β in a tokamak
While the geometry being used is kept the same, it is beneficial for a nonlinear comparison to choose a broader profile than the one used in the previous subsection, as there is less unwanted profile relaxation in a gradient-driven approach with GENE-3D. Therefore, new plasma profiles of the form are chosen, with x 0 /R 0 = 0.5, κ T i = κ T e = 6.66, κ n = 2.20, w T i = w T e = w n = 0.04 and Δ T i = Δ T e = Δ n = 0.8. Furthermore, the plasma pressure is such that β e (x 0 ) = 0.75 %.
In order to save computational time, the simulations were performed further increasing the electron-to-ion mass ratio to m e /m i = 0.01 and a finite-size parameter ρ * s = ρ s /L ref = ρ s /R 0 = 0.01. The simulation covers the radial domain 0.1 x/R 0 0.9, using buffer zones of 10 % at each side, with normalised box sizes being (L and a resolution of (n x , n k y , n z , n v || , n μ ) = (160 × 64 × 24 × 64 × 24) for GENE and (n x , n y , n z , n v || , n μ ) = (160 × 256 × 24 × 64 × 24) for GENE-3D. Finally, heat and particle sources, which are explained in Appendix A, were added with amplitudes being κ H = κ P = 0.1 and hyperdiffusion was set as η x = η y = 0.05, η z = 2.0 and η v || = 0.2.
In order to compare the results of the two nonlinear simulations, the heat fluxes Q σ for species σ are considered in the following. The heat flux itself can be split into an electrostatic and an electromagnetic component, Q σ = Q es,σ + Q em,σ , which are defined as (3.7) Here, F particle 1,σ is the perturbed distribution function of species σ , evaluated at particle position. Figure 3 shows the volume-averaged time traces of the ion and electron heat flux contributions. The fluxes are time-averaged over the interval t ∈ [100, 345]R 0 /c s , which roughly contains the same number of burst, as follows. The time traces are divided in disjoint intervals each approximately three autocorrelation times long. An average over each one of them is performed, the result of which is then used to compute an ensemble mean and variance of the heat flux. The results are presented in table 1. Both simulations are in good agreement with each other, as the relative difference between their respective ensemble mean values is below 10 % for all components.
The radial profiles of the heat flux contributions, averaged over a flux surface and in time, shown in figure 4, recover reasonable agreement between both codes as well. Here GENE and GENE-3D produce similar heat flux profiles, where ion and electron contributions are comparable in their magnitude. For both particle species, the  Both results differ by less than 6 %, giving additional confidence in the numerical implementation of GENE-3D. Finally, the spectra of the electrostatic heat fluxes at x/R 0 = 0.6 are compared in figure 5. There are some small deviations between the heat flux spectra of both codes at the smallest scale, which is the result of the different representations of the bi-normal coordinate. Nevertheless, there is excellent agreement of the spectra in the wavenumber interval where most of the transport is located. Therefore, overall the linear as well as nonlinear verification studies of GENE-3D can be considered successful.

Electromagnetic effects on ITG turbulence in W7-X
Motivated by the successful verification in axisymmetric scenarios, electromagnetic simulations of W7-X are presented in this section. As the aim is to address the effect of finite plasma-β on ITG turbulence (Snyder & Hammett 2001) in W7-X, kinetic electron simulations in the electrostatic limit of β e (x 0 ) = 10 −4 and electromagnetic simulations with β e (x 0 ) = 0.5 % were performed.  The profiles shown in figure 6 were used for all simulations discussed in this section. The variation of plasma-β was achieved by varying the reference density n ref . The specific choice of profiles localises turbulence in the centre of the plasma volume, such that one does not need to simulate the entire radial domain, with significant computational savings. The profiles are chosen to destabilise ITG modes by setting η i = L n /L T i > 1 and η e = L n /L T e = 1. The specific form of the profiles is given by (3.4), with the defining parameters set to κ T i = 4.0, κ T e = 1.0, w T i = w T e = 0.04, Δ T i = Δ T e = 0.17, κ n = 1.0, w n = 0.04 and Δ n = 0.17. For both cases, the standard configuration of W7-X was used, Linear simulations were performed to compare the growth rates of both scenarios. Since GENE-3D does not use a Fourier representation of the bi-normal direction, all modes are coupled and an initial value calculation will always only resolve the fastest growing mode present in the system. To overcome this issue, one can lower the bi-normal resolution to a point where said mode is no longer resolved, while still properly resolving the magnetic field geometry. While this method is not able to resolve local minima in the linear growth rate spectrum, it serves as a tool to calculate the maximum growth rate present in the system. An alternative is to use a numerical filter to single out the desired mode, as is done, for example, in the EUTERPE code. This is currently not implemented in GENE-3D, but it was shown in Sánchez et al. (2021) that both methods are valid approaches. Three simulations were performed for each case, using a resolution of (120 × 128 × 64 × 24) points in (x, z, v || , μ), with box lengths of size (L x , L y , L v || , L μ ) = (92.22 ρ s , 100.64 ρ s , 4.2 and hyperdiffusion parameters η x = η y = 0.05, η z = 2.0 and η v || = 0.2. The bi-normal resolution was set to (48,75,240), respectively, to resolve different linear modes. The corresponding growth rates are shown in figure 7, indicating a clear decrease induced by electromagnetic effects. The precise numerical values are listed in table 2, together with the ratio between the ion and electron heat fluxes, which is used later for comparison with the nonlinear results. One observes a reduction of the growth rates of approximately   20-25 %, with the maximum growth rate being decreased from 0.231c s /a to 0.179c s /a. For completeness, the radial and poloidal mode structures of the fastest-growing modes of both scenarios are compared in figure 8. Here one can observe that all fields peak radially close to x/a = 0.6, with only a slight shift between the electrostatic and the electromagnetic cases. Furthermore, while the electrostatic potential is highly localised around z = 0 in both cases, the parallel vector potential of the electromagnetic simulation extends over the entire poloidal domain. The stabilising properties of electromagnetic effects that were shown linearly can also be found in nonlinear simulations. For this, two simulations with the same parameters as for the linear cases were performed, except the resolution in y, which was set to n y = 120. The cost of each simulation was around 700 000 CPU hours on the Intel Xeon Gold 6148 processors of the MPCDF cluster Cobra. Figure 9 shows the volume-averaged time traces of the heat fluxes of both electrostatic and electromagnetic simulations. After approximately 200 time units, at which saturation is achieved, the electromagnetic case shows lower levels of turbulent heat fluxes than the electrostatic one. These differences in transport do not arise from nonlinear profile relaxation, but reflect the reduction of the linear growth rates, as can be seen in figure 10. In order to retain the background profiles, heat and particle sources of κ H = κ P = 0.03 were used for both simulations, resulting in nearly identical profiles for both cases, with only minor deviations from the initial ones.
In order to elaborate further on the quantitative differences in transport between both scenarios, one can compare the radial profiles of the total electrostatic and electromagnetic   peaking approximately at x/a = 0.46. Both ions and electrons undergo a reduction of their turbulent heat fluxes through electromagnetic effects by approximately 25 %, with the peak ion heat flux being reduced from 13.70Q GB to 10.34Q GB and the electron flux from 2.25Q GB to 1.65Q GB . One should mention, however, that the aforementioned reduction is observed in gyro-Bohm units only, as the increase in the plasma-β is created by a 50-fold increase in the reference plasma density. This stabilisation is in line with the decrease in linear growth rates shown in table 2. One can further compare the relative magnitudes of ion and electron heat fluxes for both cases to see that they share a common ratio of approximately Q es,e /Q es,i = 0.15. This is again similar to what was found in the linear simulations, although there seem to be slight deviations, especially for larger k y in the electromagnetic simulation.
To understand this, one has to take into account that not all modes contribute equally to turbulence, but that it is rather the low-k y modes that enter most, in particular for ITG turbulence, as was already shown, for example, in Navarro et al. (2020). To confirm this, the bi-normal wave spectra of the electrostatic heat fluxes at x/a = 0.46 are shown in figure 12. All the heat fluxes peak around k y ρ s ∼ 0.3, which is below the smallest k y that was resolved linearly, with only minor contributions from k y ρ s ∼ 1 and beyond. Overall, both linear and nonlinear simulations show a consistent reduction of ITG activity through electromagnetic effects for the given plasma-β. Nevertheless, increasing β will eventually result in the excitation of KBM turbulence, similar to the transition shown in § 3.1, which highlights the importance of electromagnetic effects for future investigations.

Conclusion
The present paper focuses on the upgrade of the GENE-3D code to capture electromagnetic effects by retaining a perturbation in the parallel vector potential. The main extension to the numerical scheme of GENE-3D with respect to Maurer et al. (2020) was the introduction of an additional field equation to compute the time derivative on the right-hand side of the gyrokinetic equation. Both linear and nonlinear verification studies show excellent agreement with the well-established tokamak code GENE, proving the correct implementation and viability of the electromagnetic upgrade.
The code was then applied to W7-X with the goal of addressing electromagnetic stabilisation of ITG turbulence. To this end, a set of analytical density and temperature profiles was considered in the standard configuration of W7-X. Linear and nonlinear simulations have been performed comparing the results between an electromagnetic case with β e (x/a) = 0.5 % and a case in the electrostatic limit. Electromagnetic effects cause a reduction by 25 % of the linear growth rates, which in turn results in a similar reduction of the saturated nonlinear fluxes. While these results were obtained under simplifying assumptions, e.g. using analytical profiles and artificial reference parameters, they demonstrate that GENE-3D is able to perform electromagnetic simulations and that finite-β stabilisation can be important for stellarators as well.
In order to approach such applications without these constraints, realistic simulations of W7-X discharges are planned for the future. Furthermore, it will be of great interest to study electromagnetic instabilities in medium-to high-plasma-β regimes in stellarators.
performed at the MARCONI-Fusion supercomputer at Cineca, Italy, as well as at the Cobra and Raven HPC systems at the Max Planck Computing and Data Facility (MPCDF), Germany.

Editor Per Helander thanks the referees for their advice in evaluating this article.
Funding 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.

Declaration of interests
The authors report no conflict of interest.

Appendix A. Additional numerical details
This appendix briefly mentions some additional numerical details of GENE-3D, a detailed discussion of which can be found in Maurer et al. (2020).
Phase-space variables are discretised in order to solve the gyrokinetic system numerically. For this, the velocity space parallel to the magnetic field lines is partitioned on an equidistant grid with zero Dirichlet boundaries, whereas the grid in the direction of the magnetic moment is distributed using a Gauss quadrature scheme with Gauss-Legendre weights and knots (Abramowitz & Stegun 1965). For the spatial discretisation, the three field-aligned coordinates x = ρ tor are introduced. Here, ρ tor = Φ tor /Φ edge is used as a radial coordinate, where Φ tor is the toroidal flux and Φ edge its value at the last closed flux surface. The bi-normal coordinate y is based on the field line label α = q(x)θ * − ϕ at a fixed flux surface, where q(x) is the safety factor, θ * is the poloidal PEST angle (Li, Breizman & Zheng 2016), ϕ is the geometrical toroidal angle and the constant C y is defined as C y = x 0 /|q 0 |, where q 0 is the safety factor at the reference point x 0 . Lastly, the parallel coordinate z describes the position along a magnetic field line. In addition, the sign of the poloidal magnetic field σ B p ensures that the parallel direction is always in the direction of the magnetic field. The magnetic field, as well as other geometric coefficients, such as the safety factor profile and metric coefficients, are provided via an interface by the Galerkin Variational Equilibrium Code (GVEC), an ideal magnetohydrodynamics equilibrium solver which follows the ideas of the well-established VMEC code (Hirshman & Whitson 1983;Hirshman & Betancourt 1991). The spatial grids in GENE-3D themselves are rectangular, equidistant meshes, which makes it particularly easy to calculate spatial derivatives using a fourth-order-accurate central-difference scheme. The boundaries in the radial direction have zero-valued Dirichlet conditions. The bi-normal coordinate is treated periodically, whereas the well-established twist-and-shift boundary condition (Beer et al. 1995) is used at the boundaries along a magnetic field line.
In order to avoid unphysical high-wavenumber modes that arise due to the central-difference approximations, one can introduce numerical hyperdiffusion terms to the right-hand side of the gyrokinetic equation to dampen these modes. The hyperdiffusion terms in the x, y, z and v || directions are all fourth-order-accurate terms with second-order stencils, meaning that they are of the form where u represents any of the previously mentioned coordinates and the damping parameter η u can be set for each coordinate individually as input by the user. By applying radial Dirichlet boundary conditions to the distribution functions, density and temperature profiles are fixed at the boundaries of the radial domain, which can cause strong turbulence and numerical instabilities at the boundaries. To avoid this, it is advisable to employ a Krook damping operator of the form to the right-hand side of (2.12), where the damping factor ν(x) is implemented as a fourth-order polynomial with compact support within a buffer region at the radial boundaries, typically 5 %-10 % at each side of the domain (Gï£¡ï£¡rler et al. 2011). As GENE-3D uses a gradient-driven approach Rath et al. 2016;Lanti et al. 2018) so far, numerical particle and heat sources have to be introduced to the gyrokinetic equation (  The symmetrisation of the distribution function with respect to v || in (A5) ensures the conservation of parallel momentum and the terms proportional to . . . FS / . . . FS avoids the unwanted numerical injection of particles (McMillan et al. 2008). The coefficients κ P and κ H are specified by the user and should be chosen to be around 5 %-10 % of the maximum linear growth rate of the system (Lapillonne et al. 2010).
The nonlinear term in (2.13), coupling v χ to F 1,σ and φ 1 , can be written in the form An Arakawa scheme (Arakawa 1997) is used to evaluate the Poisson brackets in (A8) numerically, as it ensures the conservation of free energy in nonlinear simulations (Bañón Navarro et al. 2011). Finally, in order to perform the gyroaverages efficiently, the distribution functions as well as the electromagnetic fields are discretised using a finite-element method in the directions perpendicular to the magnetic field, currently employing bicubic piecewise polynomials. Doing so, together with the finite-difference approximations of the derivatives, allows one to transform the field equations (2.9), (2.10) and (2.15) into a set of sparse, linear systems. This makes it possible to use the PETSc library (Balay et al. 1997(Balay et al. , 2019 in order to solve the equations for the electromagnetic fields. While providing an extensive set of iterative solvers already, PETSc can also serve as an interface to direct solvers such as MUMPS (Amestoy et al. 2001(Amestoy et al. , 2006 or SuperLU_ (Li & Demmel 2003).