Effects of plasma resistivity in FELTOR simulations of three-dimensional full-F gyro-fluid turbulence

A full-F, isothermal, electromagnetic, gyro-fluid model is used to simulate plasma turbulence in a COMPASS-sized, diverted tokamak. A parameter scan covering three orders of magnitude of plasma resistivity and two values for the ion to electron temperature ratio with otherwise fixed parameters is setup and analysed. Two transport regimes for high and low plasma resistivities are revealed. Beyond a critical resistivity the mass and energy confinement reduces with increasing resistivity. Further, for high plasma resistivity the direction of parallel acceleration is swapped compared to low resistivity. Three-dimensional visualisations using ray tracing techniques are displayed and discussed. The field-alignment of turbulent fluctuations in density and parallel current becomes evident. Relative density fluctuation amplitudes increase from below 1% in the core to 15% in the edge and up to 40% in the scrape-off layer. Finally, the integration of exact conservation laws over the closed field line region allows for an identification of numerical errors within the simulations. The electron force balance and energy conservation show relative errors on the order of 10−3 while the particle conservation and ion momentum balance show errors on the order of 10−2. All simulations are performed with a new version of the FELTOR code, which is fully parallelized on GPUs. Each simulation covers a couple of milliseconds of turbulence.


Introduction
Turbulence in the edge and scrape-off layer (SOL) regions of magnetically confined plasmas displays very efficient (and mostly unwelcome) transport properties [1,2].In fact, the observed levels of transport of particles and thermal energy out of the confined region by far exceed the ones predicted by collisional transport theory [3,4] even if neoclassical effects from the magnetic field geometry are taken into account.This has led to the alternative denomination of turbulent transport as 'anomalous' transport.Since particle and energy confinement are the goal of magnetic confinement fusion devices, plasma turbulence is subject to intensive research.
Numerous challenges exist when modelling plasma turbulence.For example, it is observed that relative density fluctuation levels increase from the edge into the SOL and may approach and even exceed order unity [5][6][7][8][9].This was recently also found close to the X-point region [10].This means that a linearisation of equations around a background profile is inadmissible in modelling.Avoiding such a separation between stationary profile and dynamic fluctuations in models has the additional advantage that a profile can interact with turbulence and evolve self-consistently in time.The profile is then an output of the model rather than a given input.
Furthermore, it is observed that the ratio of ion-temperature relative to electron temperature is above one in the edge and scrape-off layer regions [11][12][13].Turbulent eddies in the edge and in blobs in the scrape-off layer are of the size ρ s = √ T e m i /(eB 0 ) where T e and m i are electron temperature and ion mass respectively, e is unit charge and B 0 is the reference magnetic field strength.With ρ i = √ T i m i /(eB 0 ) ≈ ρ s (with T i the ion temperature) this leads to finite Larmor radius and polarization effects being important for the dynamics of turbulent eddies and blobs [14][15][16].
Full-F gyro-fluid models are able to evolve large fluctuation amplitudes, steep background profiles and include finite Larmor radius effects [14,[16][17][18][19]. Gyro-fluid models in general result from taking velocity space moments over an underlying gyro-kinetic model and share many of its advantages: finite Larmor radius corrections, consistent particle drifts, an energy and momentum theorem based on variational methods in the underlying gyro-kinetic model and an inherent symmetry in moment equations with regards to multiple ion species.These advantages are absent from so-called drift-fluid models that result from a drift-expansion of the Braginskii equations [20][21][22][23][24].A downside of gyro-fluid models is that closed expressions for scattering collisions and plasma neutral interactions remain an open issue, despite recent formulations in the long wavelength limit [25].Compared to gyrokinetic models, gyro-fluid models invoke a closure scheme that can be tailored to specific physical regimes of interest, e.g. the collisional regime.Such closures can be adopted at the chosen number of moments, which emerge typically from a Hermite-Laguerre expansion in velocity space of the gyro-averaged gyro-centre distribution function [17,19].The number of moment equations is usually small (2 in the present work) and the associated reduced velocity space resolution translates to a corresponding saving in computational cost over gyro-kinetic models.This implies that gyro-fluid models are more computationally efficient for parameter scans or for resolving larger plasma volumes than gyro-kinetic models.
Further challenges arise in numerical approaches to plasma turbulence.The dynamics of a magnetized plasma is highly anisotropic with respect to b, the magnetic unit vector.Fluctuations along b typically have a much larger extension L ∥ than fluctuations perpendicular to it L ⊥ ≪ L ∥ .In a numerical simulation the use of field-aligned coordinates, in particular flux-tube coordinate systems thus seems appropriate.The field alignment translates to a low spatial resolution requirement along the field line following coordinate [26][27][28].However, field aligned coordinate systems cannot include X-points in the geometry.This is a major downside as one or more X-points in the magnetic field equilibrium are a crucial ingredient to current tokamak design and in particular ITER [29].The X-point is connected to the construction of a divertor, which separates the plasma-wall interactions from the confined plasma region [1].Further, it plays a crucial role in and at least facilitates the transition to the high confinement mode [30][31][32].Correct modelling of magnetic field equilibria that include one or even several X-points is thus critical.
Two solutions to the problem exist to date.With the increase in computational resources it is possible to directly discretize and simulate model equations on non field-aligned coordinate systems [33,34].This allows simulations including Xpoints as exemplified by the GBS [35], STORM [36] or TOKAM-3X [37] codes.However, such an approach does not exploit the field-aligned character of turbulence and can thus only be used for small to medium sized tokamaks due to both strong numerical diffusion and extreme computational cost [38,39].An alternative approach is the so-called fluxcoordinate independent approach [38,[40][41][42].Here, the grid is not field-aligned while at the same time the toroidal direction is resolved by only a few points.Turbulence simulations of AUG were successfully performed with the GRILLIX code [43,44].
For the verification of codes the method of manufactured solution is often used [33,35,37,45].This method tests if the numerical solution produced by an implementation converges with the expected order to a manufactured solution.In this way (coding) mistakes of the type that prevent convergence can be found.Of course, this does not necessarily guarantee that any given simulation is indeed converged or free from any other systematic errors.An instructive counter-example is applying an explicit Euler method to the Kepler problem as outlined in [46].While the method is convergent, the numerical energy theorem is not fulfilled and the solution is qualitatively wrong.Another problem is convergence in a turbulent system.Consider that even in two-dimensional turbulence simulations numerical errors on the order of machine precision exponentially increase to order one within a short period of time [47].This is fundamentally due to the turbulent nature of the underlying model and not an issue of the numerical implementation.Thus, turbulence simulations due to their very nature cannot reach pointwise convergence after sufficiently long simulation time.We thus advocate to not only consider convergence of a numerical scheme.Properties like structure, volume or energy preservation may be equally important and help ascertain that an obtained numerical solution is indeed physical.
In this contribution we address the above challenges in a new version of the simulation code FELTOR [47,48].As opposed to the drift-fluid models discretized in the mentioned GRILLIX, TOKAM-3X, GBS and STORM codes FELTOR discretizes a full-F gyro-fluid model and thus benefits from finite Larmor radius effects, an exact energy conservation law and consistent particle drifts.Polarization effects are taken in the long wavelength limit in order to avoid inversion of an operator function [16,49].Similar to the GRILLIX code FELTOR uses an FCI scheme for its parallel dynamics but in a recently upgraded finite volume FCI formulation [42] that has significantly improved conservation properties compared to previous versions [38,40,41].For the perpendicular dynamics FELTOR chooses discontinuous Galerkin methods [50,51] in contrast to the above codes, which rely on finite difference methods.FELTOR is the only code among the mentioned ones that is fully ported to GPUs using a platform independent MPI+X implementation.Recently, all the above codes including FELTOR were part of a validation effort in TORPEX and TCV [52,53].
FELTOR allows stable simulations encompassing several milliseconds of turbulent dynamics.In a parameter scan with 12 simulations we vary the plasma resistivity and the ion to electron temperature ratio.We present techniques for three-dimensional visualisations using ray-tracing in order to gain visual intuition of the magnetic field, the density and the parallel current.In particular the field-alignment of turbulent fluctuations with L ∥ ≪ L ⊥ is visible.In order to quantitatively analyse the simulation data we introduce the flux-surface averages and integration.Numerically, these are accurately computed by transforming on a flux-aligned grid [42].We discuss flux-surface averaged density and fluctuation profiles.Afterwards, we focus on an error analysis of the implementation.Specifically, we quantify errors in exact analytical conservation laws.These include mass, parallel momentum and energy conservation as well as the electron force balance.We suggest to use volume and time integration to define a numerical error of simulation results.At the same time we are able to identify the largest and most important terms in each of the mentioned conservation equations and further in the total parallel momentum balance.Applied to the mass and energy conservation, we can compute and discuss the mass and energy confinement times.The latter relate back to our initial statement of confinement being an important goal for the magnetic fusion devices.This work is structured as follows.In section 2 we present the gyro-fluid model including resistivity and diffusive terms, the density source and boundary conditions as well as a description of the magnetic field.A parameter scan over plasma resistivity and ion temperature is setup for model simulations of the COMPASS tokamak in section 3 discussing the COMPASS magnetic field and the exact physical parameters in use.Further, we present computational performance measurements of our GPU implementation.In section 4 we present the results.We discuss three-dimensional visualisations and density and fluctuation profiles.In particular, here we show the numerical errors of mass, energy, ion momentum and parallel force balance.Finally, we discuss particle and energy confinement times computed from previously analysed terms in the mass and energy conservation equations.We conclude in section 5.

The model
In the following we denote ϕ as the electric potential, A ∥ the parallel magnetic potential, m the species mass, q the species charge, T the constant species temperature, N the gyro-centre density, U ∥ the gyro-centre parallel velocity, b the magnetic unit vector field and B the magnetic field strength.Note that all species dependent quantities m, q, N, U ∥ , T have an implied species index s that we omit in the notation.We define two magnetic field curvature vectors as well as perpendicular and parallel derivatives Notice the formulary in appendix A and the supplementary material in appendix B, which includes code documentation showing further details of the model.

Gyro-fluid moment equations
The gyro-centre continuity and parallel momentum conservation equations read for each species [17,19,54,55] (omitting the species label) ∂ ∂t The model is isothermal with constant temperature T for each species.The system is closed by the parallel Ampere law and the polarisation equation where we sum over all species.We have the density current momentum current (10) and the electric and mirror force terms The definition of the diffusive terms Λ N and Λ mNU and the resistivity R ∥ are shown in section 2.2 while the gyro-centre density source term S N is defined in section 2.3.No source is added in the parallel momentum equation.We use These are the Padé approximated gyro-average operator Γ 1 with thermal gyro-radius ρ 0 , the perpendicular magnetic field perturbation b ⊥ , the gyro-centre potential ψ and temperature T.
We keep a 2nd order accurate gyro-averaging operator Γ 1 independent of particle position that closely mimics an exponential to arbitrary order [19].The polarisation in the second term in equation ( 8) is taken in a long wavelength limit.Finite Larmor radius effects in the parallel magnetic potential A ∥ can be neglected [56].
In equation (9) we can identify the density flux parallel to the magnetic field b perturbed by magnetic fluctuations b ⊥ , followed by the E × B, the curvature and the grad-B drifts.
The first term in the momentum current equation (10) consists of the parallel momentum current quadratic in the parallel velocity U ∥ .This term is an expression of Burger's term and can lead to shocks if no parallel viscosity was added to the system.The term ∇ • (NT( b + b ⊥ )) stemming from ∇ • J mNU with J mNU from equation ( 10) can be combined with the mirror force NT( b + b ⊥ ) • ∇ ln B in equation (12) to yield the familiar pressure gradient Further, in equation (10) we have the E × B and curvature drift transport of parallel momentum.In the parallel electric force equation (11) we have the parallel and perturbed gradients of the gyro-centre electric potential ψ together with a correction due to the magnetic curvature.Even though the latter term is small it must be kept to guarantee energetic consistency.The equivalent correction also appears in the mirror force term equation (12).

Electron-ion plasma.
Even though the model is formulated inherently as a multi-species model we here only treat an electron-ion plasma, specifically with Deuteron ions (q i = e, m i ≈ 2m p with m p the proton mass).The model can also be used to simulate electron-positron plasmas [57].Multi-species gyro-fluid simulations were presented in [58,59].
We take the electron gyro-radius to be zero ρ 0,e = 0 and thus have [14,15] This is combined with neglecting the electron mass in the polarisation equation, which thus reads Note here that we denote the electron gyro-centre density as n e and gyro-centre parallel velocity as u ∥,e (as opposed to N e and U ∥,e ) to signify that these quantities coincide with the actual fluid particle density and parallel particle velocity.

Toroidal field line approximation.
First, we employ cylindrical coordinates (R, Z, φ), with φ anti directed to the geometric toroidal angle (clockwise if viewed from above) to obtain a right handed system.We denote êφ as the toroidal unit vector.We then approximate b ≈ ±ê φ in all terms that relate to the perpendicular dynamics, (i.e.all terms involving a cross product e.g.terms of the form ∇ • ( f b × v), ∆ ⊥ f, curvature operators, etc), which we denote as the toroidal field approximation [54,55].We retain b in all terms relating to the parallel dynamics (∇ • ( bf ), −T∇ ∥ N, ∆ ∥ f, etc).Note that we allow the negative sign b ≈ −ê φ for the case b φ < 0. This approximation effectively renders the drift-planes equal to planes of constant toroidal angle φ.
Technically, the approximation is justified by using the flute-mode ordering k ∥ ≪ k ⊥ [2] as well as neglecting terms that are second order in the magnetic field pitch angle ϑ b • êφ =: ± cos (ϑ) .(18) where the negative sign applies if b φ < 0. We have therefore Note that . In this way one can determine derivatives in a locally aligned coordinate system êR , êZ , b and proceed to neglecting terms of order O(ϑ 2 ) and O(k ∥ /k ⊥ ).For example we have terms of the form for functions f and g and Starting with the exact formula equation (A.9) we can insert the above expressions (19) and (20) and obtain in our ordering (21) with higher order corrections in ϑ 2 and k ∥ /k ⊥ .Similarly in the same ordering we have The curl of b reduces to ∇ × b ≈ − ±1 R êZ .The curvature operators are simplified to: and which maintains a vanishing divergence of the total curvature ∇ • K = 0.As is pointed out in [54,55] this quality is important to maintain energetic consistency in the model.The toroidal field approximation is motivated computationally.The unapproximated terms (see appendix A) contain derivatives in the φ direction.From equation (20) we see that in flute mode ordering we have k φ ∼ k R , k Z .This means that in order to accurately compute the φ derivative, the φ direction must have a similar resolution as the R and Z direction, which is highly undesirable because it is computationally expensive.In the given ordering all derivatives are in R and Z only, independent of the resolution in φ.Importantly, the φ direction can then be trivially parallelized in an implementation of the model.In fact, this is very similar to why a direct discretization of is inferior to using an aligned scheme like FCI.Computing ∂ φ requires a very high resolution in φ, which the FCI scheme avoids [38,39,42].

Resistivity and diffusive terms
Here, we discuss the terms Λ N in equation ( 5) and Λ mNU , R ∥ in equation (6).These terms take the form and We first notice that the diffusion terms have the form of total divergences ).Under volume integration these terms vanish modulo surface terms, which is important for mass and momentum conservation.Second, we notice the term −∇ • (mUj ν,N ) in the momentum diffusion (26) has the form of a velocity convection.This is a correction term that prevents energy from being generated by mass diffusion as we will see explicitly in section 4.3.2 and was suggested by for example [42,60].The consistent treatment of the diffusive terms is particularly important for the parallel ion momentum equation.The alternative variant Λ mNU,∥ := µ ∥ ∆ ∥ U ∥ + µ N,∥ mU ∥ ∆ ∥ N has the advantage that in velocity formulation Λ U,∥ = µ ∥ ∆ ∥ U ∥ /(mN) simplifies [43].However, in this formulation the term µ N,∥ mU ∥ ∆ ∥ N unphysically generates momentum, leading to artificial toroidal rotation after a long enough simulation time.Other works on drift-fluid models completely neglect the parallel ion and electron viscosities [35][36][37].
In equations ( 25) and ( 26), µ N,⊥ and µ U,⊥ are ad-hoc artificial numerical diffusion coefficients that are added to stabilize perpendicular advection and are thought to be small.In the same sense µ N,∥ represents artificial parallel diffusion necessary to stabilize the parallel advection [42].
The parallel velocity difference u ∥,i − u ∥,e := (N i U ∥,i − n e u ∥,e )/n e determines the parallel resistive term R ∥ in equation (27).The term is positive for electrons with q e = −e and negative for ions with q i = e.This form both conserves parallel momentum and vanishes for zero current but leads to a quadratic energy dissipation term only in the long-wavelength limit as we discuss in section 4.3.2.
For the parallel viscosity µ ∥ and the parallel resistivity η ∥ we use the parallel resistive and viscous terms from the Braginskii fluid equations [20].The electron-ion and ion-ion collision frequencies are given by ν ei = √ 2z 2 e 4 ln Λn e /(12π ) ( T e eV with ν ei,0 := ν ei (n 0 , T e ) as well as with ln λ ≈ 10, Ω i0 = eB 0 /m i the ion gyro-frequency and Ω e0 = eB 0 /m e the electron gyro-frequency.Finally, in order to prevent unreasonably small simulation timestep we need to impose a maximum and minimum on ν ∥,e and ν ∥,i : We emphasize that this restriction is numerically motivated.The physical implications of equation (31) are discussed in section 4.

Sources
We provide a constant (in time), toroidally symmetric influx of particles.For electrons the source term S N in equation ( 5) becomes where ω s is the source strength parameter (with unit 1/s) and n s (R, Z) is an in principle arbitrary toroidally symmetric profile, which we discuss further in section 3.2.In order to not generate potential with the source term the ion gyro-centre source needs to fulfill ) for given particle source S ne and potential ϕ, which follows from a time derivative of equation (8).We were unable to invert this equation numerically.Only in the long wavelength limit can it be inverted to yield the approximation [25] S Ni ≈ The long wavelength limit should be well-fulfilled for the source strengths we use in this work since the amplitude ω s is typically quite small.Note that the additional terms besides S ne in equation ( 33) are total divergences, which means they do not contribute to the volume integrated total particle number created by the source ´dVS Ni = ´dVS ne , where we integrate over the entire simulation volume.The surface terms are zero assuming that S ne vanishes on the surface.
A second task of the source S N is to globally ensure a minimum density.This is required since through sheath dissipation the density can in principle become arbitrarily close to zero.This is, however, both detrimental to the stability of the simulation as well as the CFL condition (and thus the allowed time step) of the simulation and in reality also never happens due to e.g.wall-recycling.For both electrons and ions we choose the additional source term (34) where H α (x) is a continuously differentiable approximation to the Heaviside function with width α.The Heaviside function ensures that this source term only acts when the density is below the lower limit.In our simulations we choose ω min = Ω i0 , n min = 0.2n 0 , α = 0.05n 0 .

Boundary conditions
Following [43] we setup boundary conditions with the immersed boundary method using volume penalization [61].In order to do this we first formally define a wall function where Ω w is the wall domain.Analogously, a sheath function χ s can be defined using a sheath domain Ω s .Both χ w and χ s are further specified in section 3.1.We have Ω s ∩ Ω w = ∅.We can then enforce boundary conditions on the wall and sheath by ∂ ∂t where (6).We choose ω s = 5 and ω w = 0.01.The polarization equation is penalized according to the immersed boundary method We do not penalize the parallel Ampere law due to numerical stability.
We choose the wall conditions N w = 0.2n 0 and U ∥,w = 0. Further, we have ϕ w = 0 and ∇ ⊥ A ∥,w = 0 for the electric and magnetic potential.The latter two are however only enforced at the domain boundaries rather than through a penalization method.We have the insulating sheath boundary conditions N sh is chosen such that ∇ ∥ N| sh = 0 for both electrons and ions.
We remind the reader that T e and T i are constants in our model.

The magnetic field
This section discusses FELTOR's general capabilities to represent toroidally symmetric magnetic fields.The specific magnetic field used for the main physical discussion in section 4 is presented in section 3.1.
In cylindrical coordinates the general axisymmetric magnetic field obeying an MHD equilibrium (µ 0 j = ∇ × B, ∇p = j × B) can be written as [62] Here, ψ p is the poloidal flux function and I(ψ p ) is the current stream function.For the sake of clarity we define the poloidal magnetic field ) and the toroidal magnetic field B t = I R êφ .Note that with a typically convex function ψ p (second derivative is positive), I(ψ p ) > 0 and the previously defined coordinate system the field line winding is a left handed screw in the positive êφ -direction.Also note that then B × ∇B points down, which for a lower single null configuration is towards the magnetic X-point, and we have the favourable drift direction (in experiments H-mode is reached more easily in this configuration [32,63,64]).
We have the contravariant components of and the covariant components In FELTOR we have various ways to represent the flux function ψ p and its derivatives.In this work we use a general solution to the Grad-Shafranov equation using Solov'ev pressure and current profiles [65,66] with A, P ψ free constants, P I = ±P ψ for A ̸ = 0 and P I arbitrary for A = 0 (purely toroidal equilibrium current).We introduce R ≡ R/R 0 and Z ≡ Z/R 0 where R 0 is the major radius and B 0 is a reference magnetic field strength.The dimensionless base functions ψpi are listed in [65].
Since equations (42) are given in terms of analytical base functions we can numerically evaluate ψ p (R, Z) and I(ψ p ) and all their derivatives at arbitrary points to machine precision, which is simple to implement and fast to execute.This translates to an exact representation of the magnetic field and related quantities, for example curvature (23), in code.In particular, the X-point(s) and O-point can be determined to machine precision via a few Newton iterations.
The choice of the coefficients c i and A determines the actual form of the magnetic field.We can for example represent single and asymmetric double X-point configurations, forcefree states, field reversed configurations and low and high beta tokamak equilibria [65,66].The scaling factors P ψ and P I are mainly introduced to maximize the flexibility e.g. to adapt the solution to experimental equilibria or to reverse the sign of the magnetic field.If one or more X-points are present, we choose c 1 such that ψ p (R X , Z X ) = 0 for the X-point closest to the Opoint that is the separatrix is given by ψ p (R, Z) = 0.

The magnetic flux, the wall and the sheath
The first step in setting up a simulation with FELTOR is to choose an appropriate magnetic field.In this work we choose to model the COMPASS tokamak and fit the magnetic flux function described in [67] with a Solov'ev equilibrium described in equation ( 42).One X-point is situated at found with a few iterations of a Newton solver).In figure 1(a) we plot the normalized poloidal flux In figure 1(b) we plot the chosen wall and sheath functions χ w and χ s , which signify the penalization regions for the immersed boundary conditions in equation (36) and equation (37).The wall region is given simply as a flux aligned region Here we choose ρ w = 1.15 for the scrape-off layer and the private flux region at ρ F = 0.97.For the sheath region we first define an angular distance φ w of each point (R, Z) to the bounding box via the integration of with initial condition (R, Z) until R((φ w ), Z(φ w )) intersects the bounding box.The intersection can be found with a bisection algorithm.The sheath is then given by where we choose φ s = 14π/32.Note that for numerical reasons we implement a continuously differentiable transition at the boundary of the regions Ω w and Ω s .Both plots in figure 1 show the numerical simulation domain in the R-Z region as

Initial profiles and sources
To initialize our simulation we choose equal for both electrons and ions such that the profile given in [67] is approximately reproduced with a peak density of n peak = 8.5 • 10 19 m −3 and a separatrix density of n sep = 10 19 m −3 .In the SOL the profile exponentially decreases to the background density of n min = 0.2 • 10 19 m −3 .The initial parallel velocity for both electrons and ions is zero everywhere except in the scrape-off layer where it varies linearly between ± √ (T e + T i )/m i with the sheath angle coordinate φ w defined in equation (45).This is to conform to the sheath boundary conditions in equation (38).
The velocity profile is initially symmetric in φ while the toroidally symmetric density profile is perturbed by small fluctuations in order to trigger turbulence.
We define the source profile in equation (32) as We choose ρ p,b = 0.55 for the source region, which is depicted as a dashed line figure 1.This form for the source is mainly intended to keep a fixed profile in the core region of our simulations and keep the total particle number inside the confined region constant.

The q-profile
We follow the methods presented in [68] and define the geometric poloidal angle Θ as the field-line following parameter around the O-point We then have with B given by equation ( 39) We can then directly integrate any field-line as from Θ = 0 to Θ = 2π.The safety factor results via With our chosen orientation of the Θ and φ angles the safety factor is positive for a right-handed field winding and negative for left-handed winding [62].Figure 2 shows the (absolute) q-profile of the chosen equilibrium.As expected the q-profile diverges at the separatrix situated at ρ p = 1.This is because B Θ = 0 at the X-point and thus the integration in equation ( 49) diverges.At the O-point around ρ p = 0 the q-profile converges to a finite value q ≈ 1.9.In the domain between ρ p = 0.4 and ρ p = 0.9 the value of q lies between 2 and 3.

A parameter scan
A dimensional analysis of our equations outlined in section 2 shows that (disregarding numerical diffusion parameters) there are 6 free physical parameters: major radius R 0 , electron temperature T e , ion temperature T i , reference density n 0 , reference magnetic field strength B 0 , and ion mass m i that map to only 5 dimensionless parameters: τ = T e /T i , η (as defined in equation ( 28)), µ = m e /m i , β = n 0 T e /µ 0 B 2 0 , and R0 = R 0 /ρ s .A priori, each simulation thus represents a one-dimensional set of physical parameters, however, the map can be made invertible by fixing the major radius (R 0 = 0.545 m in our work).
We setup parameters for in total 12 simulations as two sets of 6 simulations each.The first set uses T i = 0 while the second set uses T i = T e .The 6 simulations within each set vary the dimensionless plasma resistivity η equation (28), while keeping the plasma density n 0 = 10 19 m −3 and ρ s = 1 mm constant (as well as major radius R 0 = 0.545 m and ion mass m i = 2m p , the mass of the Deuteron).This is achieved by changing the electron temperature T e (to set η) and the The (absolute) q-profile as a function of the normalized poloidal flux ρp (43).q diverges at it approaches ρp = 1 (the separatrix) but converges to a finite value q ≈ 1.9 at ρp = 0 (the O-point).magnetic field strength B 0 (to keep ρ s ∝ √ T e /B 0 constant) as shown in table 1.This results in a constant value for the plasma beta β := n 0 T e /(B 2 /(2µ 0 )) = 10 −4 .Note that with a magnetic field strength between 0.8 T and 2.1 T and typical temperature profiles between 20 eV in the SOL and 100 eV in the edge, the COMPASS experimental values for resistivity lie somewhere in the range 4 • 10 −7 < η exp < 1.2 • 10 −5 [67].The source strength parameter ω s in equation ( 32) is constant for the duration of each simulation and chosen (differently for each simulation) such that the volume integrated source roughly matches the total density flux out of the last closed flux-surface.
We set the dimensionless parallel density diffusion necessary for numerical stability of the FCI scheme to a constant value ν ∥,N = 500.The perpendicular hyper-diffusion coefficients are set to ν ⊥,N = ν ⊥,U = 10 −3 .
The simulation domain is a rectangle in the R-Z plane chosen such that the closed field line region as well as the SOL, wall and sheath regions are captured as displayed in figure 1.It is important for the stability of the FCI scheme that the boundary of the wall region does not intersect the boundary of the simulation domain except at the sheath region.This is because boundaries that are neither perpendicular nor parallel to b are difficult to resolve for the FCI scheme and a primary source of instability in the code such that we would like to avoid them as much as possible.The domain is symmetric in The value 100 000 is chosen as a compromise between a short simulation wall time (of about a week) and a long enough, i.e. statistically significant, time series for our analysis in the following section 4. The end-time in units of ms is however different for each simulation and depends on the magnetic field strength corresponding to the chosen resistivity as depicted in table 1.Since we keep ρ s ∝ √ T e /B 0 constant, changing the electron temperature T e yields a corresponding change in B 0 and thus Ω i0 .

Performance observations
The given resolution of 2 • 10 7 grid points corresponds to an array size of 150 MB for each of the density, velocity and potential variables.With simulation data written to file at every 150 Ω −1 i0 the total file size of one simulation is about 500 GB.The grid size is of similar order of magnitude as other fluid-type simulation runs [35,44] but is a factor 10 larger than the spatial grid size of some of the largest (computationally) gyro-kinetic runs [69,70].We assume that a gyro-fluid model does not differ much in the spatial resolution requirement from its underlying gyro-kinetic model.At the same time the velocity space resolution of a gyro-kinetic code is about 10 3 grid points [69].A hypothetical gyro-kinetic code containing our model as a limit would thus require a resolution of 2 • 10 10 grid points to repeat simulations presented here.The expense in computational resources would likewise be about 1000 times higher compared to what is used in this work underlining the computational advantage of the gyro-fluid approach.
Our simulations were run on 16 NVidia V100 GPUs (equivalent to 4 nodes on the M100 GPU cluster).In table 3 we present the average runtime in seconds per Ω −1 i0 for each simulation with the error being the standard deviation.These timings include the times for input/output and diagnostics but exclude the times for initialization and restarting of the code.Typically, we achieve a computation time of 6 ± 1 s per Ω −1 i0 but outliers at 4.6 ± 0.6 s and 8.3 ± 0.2 s exist.This corresponds to a total simulation time of 7 ± 1 days for 100 000 Ω −1 i0 .The differences may be due to slightly different viscosity parameters that we chose to stabilize some simulations and subsequent smaller or larger simulation time steps.The evaluation of a single right hand side of equations ( 5) and ( 6) including solutions of all elliptic equations and evaluation of the parallel advection-diffusion terms takes about 0.20-0.25 s in all simulations.The polarization equation ( 8) is solved in typically 0.05 s and less than 0.1 s.The right hand side on average has to be evaluated 24 ± 4 times per Ω −1 i0 in all simulations corresponding to about (2.4 ± 0.4) • 10 6 evaluations per simulation.
As pointed out in our performance study [47] the observed code performance is bound by memory bandwidth and memory latencies.We emphasize that due to our structured grid approach our matrix-vector multiplications are almost as fast as vector additions since the matrix elements can be kept in cache.This and the major reduction in memory requirements that comes with it are the major benefits over unstructured grids.Of the total peak performance of 14 400 GB s −1 our implementation (of vector addition, matrixvector multiplication, scalar products) reaches on average 70%.We can compare this to the conventional Skylake partition on the Marconi cluster where one node has a theoretical peak bandwidth of 256 GB s −1 of which our implementation on average (vector addition, matrix-vector multiplications, scalar products) achieves 183 GB s −1 .With 16 nodes we thus have an available throughput of 4096 GB s −1 , which is a factor 3.5 less than what is available on 4 nodes on the M100 cluster.We see about a factor 3 in practice, i.e. a runtime of 15 s per Ω −1 i0 for the η = 10 −4 simulations and approximately 0.7 s per right hand side evaluation.

A study of resistivity and temperature
In this section we analyse the simulations previously setup in section 3.In section 4.1 we show selected three-dimensional renderings of the magnetic field, plasma density and parallel current.Following this we establish the flux surface average in section 4.2 as a diagnostics tool for a numerical error analysis of the simulations in section 4.3.We focus on the parallel acceleration in section 4.4 and mass and energy confinement in section 4.5.

Three-dimensional visualisations
Here, we present three-dimensional renderings of the magnetic field and the density and parallel current of the η = 10 −4 , T i = T e simulation.The ParaView visualisation toolkit [71] is used and all results are rendered on a NVidia RTX3090 card.Specifically, we use ray-tracing based on the OptiX path tracer using 1000 progressive passes for each image.
In order to render the φ direction smoothly we face a challenge due to the low resolution of only 32 toroidal planes.To solve the issue we temporarily extend the available data to 384 toroidal planes by interpolating along the magnetic field lines with the methods presented in [68].This allows for a smooth visualisation of the field-line following structures but the reader should keep in mind that the displayed data has been extended in this way.After the visualisation is complete we can free the used memory such that the memory consumption remains negligible towards the total size of the simulation data.We begin by showing a threedimensional rendering of the magnetic streamlines in figure 3. We use the streamline tracer module in ParaView [71] to integrate magnetic field lines of equation (39).A low-opacity isocontour of ρ p = 1.10 is plotted in order to remove spurious visualisation artifacts.The colour scale is chosen from [72,73] and is used to represent the magnetic field strength following the 'dark-is-more' bias [74] for easier interpretation.
The streamlines follow a left handed winding with the positive B direction clockwise if viewed from the top.Only magnetic streamlines in the scrape-off layer are visible, which originate at the numerical divertor at the bottom.The magnetic field strength is clearly higher on the interior side (high field side) than on the outside (low field side) following the general 1/R dependence of equation (41).As mentioned in section 2.5 the B × ∇B direction points towards the magnetic X-point and we have a favourable drift direction.

Electron density.
The electron density is depicted in figure 4. Here, we create an iso-volume for n e /n 0 ⩾ 0.22 between the angles 0 and 250 • in the φ direction.This enables the viewer to see both the field-alignment in the scrape-off layer as well as a cross-section of the perpendicular turbulence in the edge and core regions.
As a colour-scale we create a three-colour map with the help of the ColorMoves tool [72] with transitions at 0.8 and between 6 and 7.The three colours can be interpreted as visualisations of scrape-off layer (grey-blue), edge (redyellow) and core (brown-grey).Here, the core is the region where our particle source is active (cf the dashed line in figure 1).The motivation for choosing such a colour scale for the density is the large data volume spanning almost two orders of magnitude with relatively small fluctuations on top.We follow the colour-name variation metric as promoted by [75] as opposed to a colour scale that varies purely in luminance say.The latter would help to visually establish order, that is darker regions correspond to higher density values.However, we found no single-colour scale that could span the large data volume while maintaining good colour discriminability.
The scene itself shows the largest turbulent fluctuations in the core and edge regions on the low field side, especially at the outboard midplane.Fluctuations on the high field side are smaller in perpendicular extension.This points towards a ballooning mode.Further, we notice that fluctuations are highly elongated vertically at the top of the domain as well as at the bottom around the X-point both in the edge as well as the scrape-off layer.The scrapeoff layer fluctuations appear field aligned judging from the form of the contours in between the two poloidal planes.The next visualisation is the parallel current j ∥ = e(N i U ∥,i − en e u ∥,e ) in figure 5. We create two separate iso-volumes for j ∥ : one for j ∥ /(en 0 c s ) ⩾ 0.5 and one for j ∥ /(en 0 c s ) ⩽ 0.5.Here, we use c s = √ T e /m i = 4.64 • 10 4 m s −1 .Two separate colourmaps are chosen for each region; a blue one for the negative and a red one for the positive values.Both colourmaps begin at ±0.5 and are truncated at ±1 (actual values lie between ±4.7).In order to guide the viewer we plot a low-opacity iso-contour of ρ p = 1 (the separatrix).
The resulting image highlights the localized 'field-aligned tubes' in which current flows in the simulation.These tubes have a typical extension of about 5 mm and thus carry a current of approximately 25en 0 c s ρ 2 s ≈ 2.5 A. It is further visible that the current is positive (flow direction clockwise viewed from above) mainly on the high-field side and negative mainly on the low-field side.However, a couple of individual current lines of the opposite signs are discernible in either region.Few current lines exist in the scrape-off layer and only close to the separatrix.
The large scale regions of positive and negative current on the high-and low-field side could be a form of Pfirsch-Schlüter current.A Pfirsch-Schlüter current is typically associated to the parallel part of an equilibrium current.Now, the background magnetic field in equation ( 39) is externally given and any equilibrium current that generates it is absent from our model.At the same time the observed large-scale currents in figure 5 lead to an associated magnetic potential A ∥ that adds a small poloidal magnetic field b ⊥ equation ( 14) to the background field.With the toroidal field approximation 2.1.2we have A comparison with equation (39) reveals that ψ p + RA ∥ can be viewed as a modified magnetic equilibrium flux function.However, in a contour-plot RA ∥ barely changes the contours of ψ p (see supplementary material in appendix B).We conclude that the perturbation RA ∥ of the flux-function is rather small.

The flux surface average-profiles and fluctuations
Before we can turn to an error analysis of our simulations we first need to establish appropriate integration routines.More specifically we here want to compute so called flux-surface averages and integrals.The flux-surface average is defined as a differential volume average according to [62]: where H(x) is the Heaviside function.In order to accurately integrate equations ( 51) and ( 52) we use the methods described in [68].The first step is to construct a flux aligned coordinate system as we show in figure 6.
There are several numerical pitfalls that should be considered when numerically constructing a flux-aligned grid.As pointed out in [68,76] the volume element in flux-aligned grids diverges and care must be taken when constructing such grids close to the X-point.This is especially true if the fluxsurface average of the separatrix (or a surface close to it) is to be computed.We follow [76,77] for the construction of our grid.
In flux-aligned coordinates η, ζ, φ the flux-surface average simplifies to where √ g is the volume element in the flux aligned coordinate system.
The numerical integration in the φ direction is straightforward.The resulting toroidal average can be interpolated onto the flux-aligned grid displayed in figure 6.Then, equation ( 53) can be used to compute the flux surface average.Since the grid in figure 6 exists also outside the last closed flux surface, we can use equation ( 53) to compute flux-surface averages in the scrape-off layer as well.
In figure 7 in the top half we show the flux-surface averages of the density ⟨n e ⟩ as a function of ρ p defined in equation (43).In fact, we show a time averaged ⟨n e ⟩ profile for all simulations.The average profiles for T i = 0 and T i = T e are visibly very similar.For the high resistivity simulations η = 3 • 10 −4 and η = 10 −4 (both T i = 0 and T i = T e ) the average profile appears linear in ρ p up to the separatrix at ρ p = 1.The remaining simulations have accumulated density in the core at about ρ p < 0.4.This is the region where the source S ne is active and continuously increases the density, which also translates to a large variation amplitude in the core.The edge and scrape-off layer at 0.95 < ρ p < 1.075 are shown enlarged.The density on the separatrix increases with resistivity from 0.5 • 10 19 m −3 to about 1.5 • 10 19 m −3 for both T i = 0 and T i = T e simulations.Afterwards, in the scrape-off layer at ρ p > 1 the density sharply drops.Notice that the black dashed line in the enlarged region signifies the minimum density n e,min = 0.2 • 10 19 m −3 in equation (34).The average densities thus cannot reach below the artificially enforced lower boundary.It may be preferable to run simulations with lower n e,min to study if the lower resistivity simulations converge at a different value, however then also the parallel viscosities ν ∥ must be adapted in equation (31) in order to not also lower the CFL condition.
We define the relative fluctuation amplitude as In the lower part of figure 7 we show the time averaged σ ne for our simulations.Again, both the T i = 0 and T i = T e simulations exhibit similar behaviour.The fluctuation levels in the core region lie between 10 −3 and 10 −2 at the smallest ρ p where higher resistivity corresponds to higher fluctuation levels.The relative fluctuation amplitudes increase for all simulations to about 15% at the separatrix.There is sharp increase in fluctuations for ρ p > 1 to a maximum of 35% for T i = 0 and 40% for T i = T e , visible in the enlarged regions of figure 7.
Furthermore, between about 1 < ρ p < 1.01 the amplitudes for all simulations overlap before they decrease again at about ρ p = 1.02.The small resistivity simulations decrease furthest in fluctuation amplitudes.
The decrease in fluctuation amplitude again is a numerical effect of the enforced minimum density in equation (34).From experimental measurements we expect the fluctuation amplitudes to increase in the SOL [6][7][8][9].
Note that setting cold ions T i = 0 in the model removes finite Larmor radius effects with Γ 1i = 1 equation ( 13) as well as the parallel ion pressure gradient T i ∇ ∥ N i and part of the ion curvature fluxes in equations ( 9) and (10).We conclude that these terms have only minor influence on the density and fluctuation profiles, which appear similar for both T i = 0 and T i = T e .With respect to finite Larmor radius effects we can further estimate from figure 5(a) radial fluctuation extension of l ⊥ ≈ 5ρ s and thus have (k ⊥ ρ s ) 2 ≃ (ρ s σ ne /l ⊥ ) 2 ≃ (0.4/5) 2 = 0.0064 ≪ 1.Small scales on which finite Larmor radius effects become important are not present in our simulations.However, we may not yet conclude that finite Larmor radius effects are unimportant in general.Higher order polarization terms currently not present in equation ( 8) may lead to the so-called gyro-amplification strongly exciting small scales with k ⊥ ρ s ≃ 1 in the E × B vorticity, as is shown for isolated blob simulations in [14,16].
The observed radial profiles for density and its fluctuations can be tentatively compared with [78] where a nonisothermal drift-fluid model is used to simulate the turbulent dynamics in a limiter configuration using buffer regions to exclude the core region from the simulation domain.There, the fluctuation level at the separatrix peaks only for small resistivities.Furthermore the separatrix densities are highest for smallest resistivities instead of largest resistivities as in our case.This is likely a consequence of how the source term S N depends on η.In the present case the source strength is adapted (see table 1) such that the density profiles across simulations remain similar, while [78] keeps an absolute source strength.

Error analysis of conservation laws
With a reliable way to compute flux-surface averages and volume integration we can now turn to defining a suitable error norm for a numerical error analysis.First, we again emphasize that due to the turbulence nature of our simulations, we cannot show pointwise convergence.In fact, in [47] it is shown that even computational errors on the order of machine precision in two-dimensional simulation exponentially increase to order one within a short period of time.We here therefore follow a different strategy where we compute the volume and time integrated error of conservation laws.

Assume that our model equations in section 2 allow for a local analytical balance equation of the form
that is a sum of individual terms t i balances to zero.First, we define a time average via The time interval [t 0 , t 1 ] in equation ( 56) will in the following section 4.3.1 be manually defined for each simulation by identifying a saturated turbulence state.Under a further volume integration we can convert the t i to The spatial integration region in equation ( 57) is chosen as the closed field line region Ω := {(R, Z, φ) : Z > Z X ∧ ρ p (R, Z) < 1} and shown in colour in figure 6.Note that once we have the flux-surface average ⟨t i ⟩ on a sufficiently fine grid in ψ p we can integrate We then have ∑ i T i = 0 analytically, however, numerically due to discretization errors we usually have where E is the total numerical error and T num i is the numerical result given by the discrete version of equation (57).We would consider the conservation law well fulfilled numerically, if E is small compared to the which motivates the definition of the relative global error The error E consists of the contributions E i of the errors of each individual term We are interested in the error for each term, however, given E we a priori cannot deduce E i .In order to get an error estimate nevertheless, we here assume that the error contribution E i of each term is determined by its magnitude |T num i |, i.e. we define It is important to state that in our approach the terms t num i and T num i are numerically computed from the data available in post-processing independently from the numerical methods used during the simulation.For example, for the term ∇ • (n e u E ) in the mass conservation an upwind scheme is used in the simulation code while simple centred differences are used in post-processing.A part of the error E is thus the discretization error due to using different methods for computing the same term.
The main advantage of our method is that it applies to the same simulation data that is used for data-analysis, i.e. it provides numerical errors E i for terms T i that are used for physical interpretation.For example, we can quantify the numerical errors of the terms we use for the mass and energy confinement times in section 4.5.At the same time the method makes a statement of how large or important these errors are relative to the quantity of interest.Our approach may also catch coding mistakes provided the caused error is much larger than the discretization error.This is because the simulation code and the diagnostics code are independent and serve as consistency checks for each other.Finally and importantly, our method catches certain systematic errors in a conservation law.For example, terms of the form ∇ • ( f b) vanish under flux-aligned volume integration, which is not necessarily true numerically.Any numerical diffusion out of the confined region will show up as an error contribution.
A disadvantage of the method is that it cannot determine which contribution exactly causes a certain error, only the size of the total error relative to a certain quantity of interest.It is up for discussion whether this size is sufficiently small or cause for concern.Overall, we consider any error below 1% as excellent, while anything above merits further discussion.
We now analyse the mass conservation in section 4.   with where  [62] and thus Equation ( 64) signifies that the total flux out of the last closed flux surface is given by either the flux-surface average of the radial current ⟨ j • ∇v⟩ or the volume integral of the divergence ´⟨∇ • j⟩ dv.Assuming that we can numerically compute j then equation (64) gives two methods to numerically compute the total flux out of a flux surface: (i) we can numerically compute ∇ • j, compute ⟨∇ • j⟩ and finally numerically integrate ´⟨∇ • j⟩ dv or (ii) we can form j • ∇ψ p and then compute ⟨ j • ∇v⟩ = ⟨ j • ∇ψ p ⟩ (dv/dψ p ), where dv/dψ p can be computed highly accurately from the Jacobian of the flux aligned grid [68].Both methods converge because the involved numerical derivative and integration formulas converge and equation (64) guarantees that they converge to the same solution.We tried both methods in our analysis and it turns out that the former method is slightly superior to the latter.The global errors E are systematically lower using volume integrals of divergences than using the average radial flux.Furthermore, we see in figure 9 that ⟨j • ∇v⟩ exhibits slightly more fluctuations than the direct volume integral.
In figure 10 we plot the volume integrated terms of the mass conservation (61) as a function of time for the T i = T e and η = 10 −4 simulation.We immediately see that the two largest actors in this figure are the E × B flux ⟨ j E • ∇v⟩ on the last closed flux surface and the density source ´Sne dV, which is constant throughout the simulation.The time derivative of the total mass fluctuates around zero.Note that the remaining terms including the error given by the sum of all terms ∑ i t i are too small to be visibly different from zero in the plot.
Further, notice that the flux surface average ⟨∇ • (j 0 b)⟩ = d dv ⟨ j 0 b • ∇v⟩ = 0 vanishes for any parallel current j 0 b.Any deviation from zero is thus purely numerical.This applies in particular to the terms ∇ • j ∥ and Λ ne,∥ in equation (61).In our recent work in [42] we individually study the deviations from  zero in those terms and find them to be negligibly small.We will thus here and in the following ignore parallel terms accepting that they may contribute to the errors visible in figure 8.
From the E × B flux in figure 10 we manually identify a time interval where fluctuations appear around a constant average.We do this for all 12 simulations.This allows us to identify suitable t 0 and t 1 = t end in equation ( 57) and thus we can compute the relative global error in equation (59).We plot terms together with error bar from equation ( 60) in figure 11.The left plot shows simulations with T i = 0 for the various resistivities η and the right plot shows corresponding simulations with T i = T e .We can immediately confirm that the E × B flux as well as the source term are the largest terms for all simulations while the time derivative follows with lesser importance.Note here that the density source strength ω s in S ne in equation ( 32) was chosen differently for each simulation.The magnetic flutter term as well as the curvature flux and the perpendicular diffusion terms have negligible importance on the evolution of the global mass balance.We emphasize that this does not necessarily imply negligible importance on the local dynamics just that the volume integrated mass balance is unaffected.
The relative errors in the terms are invisible in this plot, which is why we separately plot these in figure 8.There we see that the relative error of the terms in the mass conservation is at an excellent maximal 3% for all simulations and below 1% for simulations with η > 10 −5 .

Energy theorem.
The terms of the energy theorem are with ) where in the energy flux j E we neglect terms containing time derivatives of the electric and magnetic potentials and we sum over all species.The energy density E consists of the Helmholtz free energy density for electrons and ions, the E × B energy density, the parallel energy densities for electrons and ions and the perturbed magnetic field energy density.In Λ we insert the dissipative terms of section 2.2 and use The dissipation term can be further simplified to where we use The dissipation term thus consists of a diffusive energy current under a total divergence and a dissipation contribution.Focusing on the parallel diffusion terms we find for the dissipative contribution: The first two terms are always negative and thus always dissipate energy.The last term containing the potential vanishes under species summation at least to zeroth order with n e ≈ N i and ψ i ≈ ϕ.
The term R E is approximately quadratic in the sense that R E ≈ −η ∥ j 2 ∥ , which is the familiar Joule heating term.Since we have an isothermal model this term appears as an energy dissipation term.The source term S E dissipates parallel kinetic energy −0.5mU 2 ∥ S N < 0 but generates free energy ln NS N > 0. The integration region in time remains unchanged and we can compute the time and volume integrated terms equation (57) with error bar equation ( 60) in figure 12.The relative errors of the terms must again be read from figure 8 and are below 1% for all simulations.The global relative error in energy is generally a factor 2-5 smaller than the error in mass.
In figure 12 we see that the energy source S E is the largest (and only) negative contributor in the equation.From equation (69) we see that it is in fact the density source S ne that translates to a source of energy.The magnitude of the energy source decreases by approximately a factor 10 from smallest to highest resistivity.Since the density source does not vary much in figure 11, this is likely a simple consequence of the decreasing electron temperature in our parameter scan in table 1.The energy source is balanced by the energy flux out of the last closed flux surface j E , the parallel energy dissipation Λ E,∥ , the Joule heat R E , the perpendicular energy dissipation Λ E,⊥ and the energy gain ∂ t E.
Few clear trends with resistivity can be inferred from the plot.The parallel energy dissipation is systematically larger than the perpendicular energy dissipation.The resistivity term R E becomes relatively less important for smaller resistivities η than for higher resistivities.For T i = 0 the energy gain ∂ t E is most important for small resistivities η < 10 −4 but least important else.The energy flux term j E is most important for η ⩾ 10 −4 but small compared to the other terms for η < 10 −5 .
The fact that the energy gain ∂ t E is the largest contributor especially for small resistivities indicates that the system constantly gains energy during our simulations.It thus appears that these simulations have not yet reached a steady state where ∂ t E = 0 on average.Further analysis reveals that the energy at the end of the small resistivity simulations is up to 10% higher than at the beginning.Longer simulation runs need to confirm the trends shown in figure 12  In the parallel momentum equation ( 6) for ions we insert the mirror force equation ( 12) and use with ion momentum current as well as resistivity term and the parallel electric force Note that the total divergences ∇ • j mNU,∥ and Λ mNU,∥ , parallel flux and viscosity terms, again vanish exactly under the flux-surface average.We plot the terms of the ion momentum equation in the top half of figure 13.Again, the error bars are invisible and are separately plotted in figure 8.There we find relative errors for T i = T e between 10 −3 and 3 • 10 −2 .Each term of the ion momentum equation thus has a relative error of maximal 3%.This is true also for T e = 0 and η > 10 −5 simulations.However, for T i = 0 and η ⩽ 10 −5 the relative error climbs to about 10%.This can be reasoned in the smallness of the terms in figure 13, i.e. the absolute error of the equation remains the same across simulations but the term 59) is small for T i = 0 and small η.
In figure 13 the largest positive term is the parallel electric force eN i ∇ ∥ ψ.To this add negative contributions from the gauge term eN i ∂ t A ∥ and the magnetic flutter eN i b ⊥ • ∇ψ.The resistivity term R ∥,e , as expected, makes a significant contribution only for large η > 10 −5 for both T i = 0 as well as T i = T e .The E × B flux is the final significant term and decreases in magnitude with η.The absolute value is however larger for T i = T e than for T i = 0.
For T i = T e and small resistivities the term m i ∂ t N i U ∥,i is the largest positive term.This indicates positive acceleration, while for large resistivites η > 3 • 10 −5 there is acceleration in the opposite direction.For T i = 0 the same trend can be observed, however, the magnitude of the term is about a factor 10 smaller than for the T i = T e simulations.We will discuss this further in section 4.4.In the bottom half of figure 13 we plot the terms of the parallel force balance.The relative global error of this equation is generally the smallest among all the equations that we test.
In figure 8 we see that the error is of excellent orders 10 −4 and 10 −3 , which lies in the range of the value for m e /m i = 2.7 • 10 −4 .This confirms that at least under volume integration equation ( 77) is very well fulfilled even if it is not analytically exact.
Analogous to the ion momentum equation the largest term in the electron force balance is the parallel electric force en e ∇ ∥ ϕ.Notice here that the colours of figure 13(top) and (bottom) coincide for analogous terms.In fact, visually the terms en e ∇ ∥ ϕ, en e b ⊥ • ∇ϕ and en e ∂ t A ∥ , i.e. all terms of the electric field are indistinguishable from eN i ∂ t A ∥ , eN i b ⊥ • ∇ψ and eN i ∇ ∥ ψ.We will use this to further study the total momentum equation in section 4.4.

Parallel acceleration
Figure 13 is visually overburdened due to the number of displayed terms and thus hard to physically interpret further.Thus, we here simplify the discussion by focusing on the total momentum balance.First, we see in figure 13 that the electron and ion components of the electric field and the resistivity are visually equal.Neglecting those terms we sum the ion and electron momentum equations to get We further neglect the term ∇ • j mNU,A and Λ mNU,⊥ and approximate The result is shown in figure 14.
The error bars in figure 14 are visible in particular in the T i = 0 plot, however the plot is easier to interpret than figure 13.We now clearly see the positive acceleration in the T i = T e plot for η ⩽ 10 −4 .For η ⩾ 10 −4 the parallel acceleration is negative.The T i = 0 plot shows the same trends but the acceleration is more than a factor 10 smaller than for T i = T e .
Four candidates explain the observed accelerations.The E × B flux of parallel momentum is negative signifying that positive momentum is lost to the plasma (or negative momentum enters the plasma) via the radial transport.The E × B flux decreases in magnitude with η for both T i = 0 and T i = T e but is about a factor 2 − 4 larger for T i = T e than for T i = 0.For T i = 0 the two terms ∇ • j mNU,C and m i N i U ∥,i K ∇× b • ∇ψ are close to zero for all η.The only remaining term for T i = 0 is thus the parallel gradient T e ( b + b ⊥ ) • ∇n e , which remains roughly constant in η.
For T e = T i the term (T e + T i )( b + b ⊥ ) • ∇n e is positive but much smaller than the curvature contribution.The second curvature term ∇ • m i N i U ∥,i K ∇× b • ∇ψ is strongly negative for η < 10 −4 but jumps to a positive contribution at η = 10 −4 thus facilitating the associated negative acceleration.The term ∇ • j mNU,C in figure 14 represents the total flux of ion momentum through the last closed flux surface by curvature drifts, while the term m i N i U ∥,i K ∇× b • ∇ψ appears as a drift correction to the parallel electric force term (11).In our previous theoretical analysis both curvature terms were neglected as small [79] but for T i = T e each term has similar contribution in magnitude to the radial E × B momentum flux.
Again, similar to the energy balance it can be argued that the appearance of a net parallel acceleration indicates that our simulations have not all reached a steady state.Our analysis thus gives only indications of the situation where the parallel momentum saturates.We would expect that the positive/negative parallel acceleration then translates to a net positive/negative parallel momentum.

Mass and energy confinement times
From our analysis of the mass conservation equation in figure 11 and the energy conservation equation in figure 12 it is straightforward to extract confinement times.As we explained before in equation ( 64) the volume integral of ∇ • j yields the total flux out of the closed fieldline region ´j • dA.We thus start with the definition of the total particle number and energy within the confined region We can then compare these with the total loss of particles and energy.The particle loss is simply the total flux j ne (equation ( 62)) integrated over the last closed flux surface.We can neglect the diffusive transport from figure 11 as close to zero.The losses of energy consist of the energy flux out of the last closed flux surface, but also of the energy dissipation through diffusion and the resistivity.We thus define ) In figures 15 and 16 we present the resulting values for our simulations.Note that the total particle number ⟨M⟩ t = (2.3 ± 0.1) • 10 19 is roughly constant for all simulations.From figure 12 we should keep in mind that the total energy has not saturated for small resistivities.Longer simulations are needed to confirm the trends in the following discussion.The error bars are computed from the fluctuation amplitudes of all quantities in equations ( 81) and (82).The relative numerical errors are negligible at 1% as established in section 4.3.Two regimes are visible in both plots with a transition at η crit ≈ 5 • 10 −5 for both T i = 0 as well as T i = T e .The mass confinement times in figure 15 reach roughly constant values for η < 3 • 10 −5 while for η > 10 −5 there is a decrease of confinement with increasing resistivity.The drop in mass confinement above the critical η could be related to the discussion of the density limit [80,81] in the operational space of tokamaks.The constant regime should be regarded tentatively as the fluctuations are particularly large in this regime, especially for T i = T e .The values for T i = 0 are a factor √ 1 + T i /T e larger than the ones for T i = T e within the error bars.We can tentatively fit a power law of for η < 5 • 10 −5 η −1/3 for η > 5 • 10 −5 (83) where c M (n 0 , m i , R 0 , ρ s ) signifies the unknown dependency on the parameters n 0 , m i , R 0 and ρ s that we kept constant during our parameter scan.We remind the reader here that the values for both T e and B 0 decrease for increasing η in our parameter scan as seen in table 1 as well as our dimensional analysis in section 3.4.
For the energy we see a clear maximum in the confinement time at η = 3 • 10 −5 .The fluctuations are systematically smaller for the energy confinement times than for the particle confinement times.However, the energy confinement times are also approximately a factor 100 smaller than the mass confinement times.This may be due to the fact that we have an isothermal model where Joule heat is not converted to an increase in temperature and is instead lost to the system.A tentative fit reveals τ E = c E (n 0 , m i , R 0 , ρ s ) √ 1 + T i /T e { η +1/4 for η < 3.5 • 10 −5 η −1/3 for η > 3.5 • 10 −5 (84) where similar to equation (83) the factor c E (n 0 , m i , R 0 , ρ s ) encapsulates a yet unknown dependence on the parameters n 0 , m i , R 0 and ρ s .The existence of a critical value for the plasma resistivity η crit ≈ 5 • 10 −5 for both mass and energy confinement points towards two different turbulent regimes above and below the critical value.Various candidates are discussed in the literature with the most likely ones being drift-wave turbulence for small η and resistive ballooning type turbulence for high η [81][82][83].
According to [83] the transition between the two regimes happens at the resistive ballooning threshold at α t,crit = 1 with turbulence parameter α t := ηq 2 R 0 /ρ s ≈ 5 • 10 3 η.With η crit = 5 • 10 −5 we obtain α t,crit, num ≈ 0.25, which is only a factor 4 away from the theoretical prediction.The difference may be explained by geometrical factors like the presence of the X-point.
There is an apparent discrepancy in this explanation however, insofar the transport in drift-wave turbulence reduces for small η (converging to the adiabatic case) and thus the confinement time should increase for decreasing η instead of remaining constant.An explanation for the observed plateau in the mass confinement time could be so-called reactive instabilities, which are independent of η and are due to a finite electron inertia [2].Reactive instabilities are unphysical insofar they are an artefact of an isothermal gyro-fluid model and have no gyro-kinetic counterpart where Landau damping counteracts the effect of electron inertia.Note that this does not contradict figure 13 where the electron inertia effect vanishes under volume integration.Locally, the electron inertia may still be important.

Conclusion
We present a new version of the three-dimensional gyrofluid turbulence code FELTOR.12 simulations covering several milliseconds with different values for plasma resistivity and ion temperature and fixed values for plasma density and gyro-radius are setup, analysed and discussed.An efficient implementation on GPUs allows for simulation runtimes of about 1 week per simulation covering a couple of milliseconds of turbulence.Simulation results are analysed using volume and time integrated conservation laws, mass, energy, momentum and force balance.Relative errors are generally below 1% for energy conservation and force balance while for mass and momentum conservation the errors climb to about 3% as seen in figure 8.Only in the ion momentum balance and for vanishing ion temperature and small resistivity do we see relative errors of about 10%, which is reasoned in the smallness of the parallel acceleration compared to T i = T e simulations with at the same time equal absolute errors.
We systematically investigate the importance of the terms in the parallel momentum generation where we find that for increasing resistivity the direction of acceleration is swapped.This is caused mainly by an interplay of decreasing E × B momentum transport and curvature drifts across the separatrix.The analysis of the momentum density m i N i U ∥,i is related to intrinsic toroidal rotation in tokamaks and the angular momentum density m i N i U ∥,i R [79,84].At the same time the appearance of net acceleration indicates that the parallel momentum has not yet saturated during the simulation time and longer simulations are needed to confirm the trends shown here.A detailed analysis of rotation profiles and the angular momentum balance is here postponed to future analysis.
Similar transitions from a low resistivity regime to a high resistivity regime happen for the mass and energy confinement times.Beyond the critical resistivity the mass and energy confinement decrease with increasing resistivity.Below it, the mass confinement remains roughly constant, while the energy confinement decreases with decreasing resistivity.This behaviour could be explained by so-called reactive instabilities, which are an artefact of electron inertia in isothermal gyrofluid models and have no gyro-kinetic counterpart.A dynamic electron temperature should help counteract this effect in future works.The transition from drift-wave turbulence to resistive ballooning roughly coincides with the value predicted by the literature.Further parameter studies in ρ s and n 0 need to clarify the unknown dependence factors c M (n 0 , ρ s ) and c E (n 0 , ρ s ) in the observed scaling laws for τ M (η) (83) and τ E (η) (84).
We did not find qualitative differences in our analysis of mass, energy or momentum confinement between the T i = 0 and T i = T e simulations.Also, the density and fluctuation profiles are similar in both scenarios.Quantitative differences exist, in particular a factor √ 1 + T i /T e in the mass and energy confinement times.Small scales with ρ s k ⊥ ≃ 1 on which finite Larmor radius effects become important are absent in our simulations.This may be due to the absence of arbitrary wavelength polarization terms in equation (8), which are known to strongly excite small scales ρ s k ⊥ ≃ 1 in the E × B vorticity, an effect called gyro-amplification [14,16].
The capability of running numerically stable simulations for a set of different parameters with FELTOR is an important milestone.We offer a first high level analysis of the run simulations and quantify numerical errors, leaving many questions open for future work as outlined above.Furthermore, various physical model improvements can be added fairly straightforwardly within the FELTOR framework.These include for example, dynamic temperature equations [15], plasma-neutral collisions [25], arbitrary order polarisation terms [16,19] and more.
opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission.Neither the European Union nor the European Commission can be held responsible for them.This work was supported by the UiT Aurora Centre Program, UiT The Arctic University of Norway (2020).This research was funded in whole or in part by the Austrian Science Fund (FWF) [P 34241-N].For the purpose of Open Access, the author has applied a CC BY public copyright license to any Author Accepted Manuscript (AAM) version arising from this submission.This work was supported by a research Grant (15483) from VILLUM Fonden, Denmark.

Figure 1 .
Figure 1.Calibration of the simulation box.The normalized magnetic flux ρp = √ (ψ p,O − ψp)/ψ p,O on the left and the wall and sheath regions on the right.The magnetic flux ψp is modified to a constant inside the wall region.On the right plot colours range linearly from 0 to 1. Two contour lines indicating the wall at ρp = 0.97 in the private flux region and ρp = 1.15 in the scrape-off layer region are plotted in solid black lines.The separatrix ρp = 1 and the boundary of the source region at ρp = 0.55 in the core are plotted in black dashed lines.

Figure 2 .
Figure 2.The (absolute) q-profile as a function of the normalized poloidal flux ρp(43).q diverges at it approaches ρp = 1 (the separatrix) but converges to a finite value q ≈ 1.9 at ρp = 0 (the O-point).

Figure 3 .
Figure 3. Streamlines of the magnetic field vector B integrated and visualised in ParaView [71].One low-opacity iso-contour of ρp = 1.10 is plotted (corresponding to ψp = 4).The positive B direction is clockwise if viewed from above and the field-line winding is left-handed.B × ∇B points towards the magnetic X-point and we have a favourable drift direction.

Figure 4 .
Figure 4.The electron density ne at 5.18 ms for η = 10 −4 , T i = Te = 7.8 eV and B 0 = 0.4 T. We show an iso-volume of ne/n 0 ⩾ 0.22 and choose a wave colourmap constructed with the ColorMoves tool from [72] mapped to logarithmic density values.The three colour regions (blue-grey, red-yellow and brown-grey) roughly coincide with the three regions scrape-off layer, edge and core/source region (cf figure 1(b)).(A video for this figure is available).

Figure 5 .
Figure 5.The parallel electric current j ∥ /(en 0 cs) at 5.18 ms for η = 10 −4 , T i = Te = 7.8 eV and B 0 = 0.4 T. We plot two isovolumes j ∥ ⩽ −0.5en 0 cs and j ∥ ⩾ 0.5en 0 cs.The colour-scale is cut at −1 and 1 respectively.A translucent contour of the separatrix ψp = 0 is shown.Current mainly flows in field-aligned tubes.Each tube has a typical extension of 5 mm and thus carries approximately 25en 0 csρ 2 s ≈ 2.5 A.

Figure 6 .
Figure 6.The flux-aligned grid (with 20 × reduced resolution to see the grid points) used for the computation of flux-surface averages and flux-volume integration.The closed field line region Ω for our analysis is shown in blue and contains a volume of 0.5 m 3 .The grid allows for a definition of a flux-surface average outside the separatrix.

Figure 7 .
Figure 7.The time averaged density profiles (top) and the relative fluctuation amplitudes (bottom) for T i = 0 (left) and T i = Te (right) as a function of ρp equation (43).Colours correspond to different values of plasma resistivity η.The separatrix corresponds to ρp = 1.The edge and scrape-off layer regions 0.95 < ρp < 1.075 are shown enlarged.

Figure 8 .
Figure 8.The relative global errors as defined by equation (59) of the terms in the mass conservation in section 4.3.1, the energy theorem in section 4.3.2, the parallel momentum balance in section 4.3.3 and the electron force balance in section 4.3.4 for T i = 0 (left) and T i = Te (right).
3.1, the energy theorem in section 4.3.2, the parallel momentum balance in section 4.3.3 and the electron force balance in section 4.3.4.The resulting relative global errors are presented in figure 8 .

4. 3 . 1 .
Mass conservation.The electron density equation (5) directly yields the particle conservation ∂ ∂t n e + ∇ • j ne − Λ ne − S ne = 0 (61) we split the density flux into the E × B flux j ne,E := n e b × ∇ϕ /B, the curvature flux j ne,C = −n e TK/e − m e n e u 2 ∥,e K ∇× b, parallel flux j ne,∥ = n e u ∥,e b and magnetic flutter flux j ne,A = n e u e,∥ b ⊥ .The diffusive part consists of Λ ne,⊥ = −µ ne,⊥ ∆ 2 ⊥ n e and Λ ne,∥ = µ ne,∥ ∆ ∥ n e .Analytically, the volume integral of divergences like ∇ • j ne yields the total flux out of a flux surface.It holds ⟨∇ • j⟩ = d dv ⟨ j • ∇v⟩

Figure 9 .
Figure 9.A comparison of two ways to numerically compute the flux ´dA • j ne,E equation (64).Shown is the radial E × B flux as a function of ρp equation (43) for the simulation T i = Te, η = 10 −4 at t = 5.18 ms.

Figure 10 .
Figure 10.The time evolution of volume integrated terms in the mass conservation equation for T i = Te and η = 10 −4 .The length of the shaded regions signifies the time interval which we consider for our statistics while the widths signify the standard deviations within that region.

Figure 11 .
Figure 11.The mass conservation equation (61): volume integrated and time averaged terms equation (57) with error bar equation (60) for T i = 0 (left) and T i = Te (right).The error bars are too small to be visible in the plot and are separately shown in figure 8.

Figure 12 .
Figure 12.Energy conservation equation (66): the terms equation (57) with error bar equation (60) for T i = 0 (left) and T i = Te (right).The error bars are too small to be visible in the plot and are separately shown in figure 8. .

Figure 13 .
Figure 13.The parallel momentum balance (top) equation (73) and the parallel electron force balance (bottom) equation (77): the terms equation (57) with error bar equation (60) for T i = 0 (left) and T i = Te (right).The error bars are too small to be visible in the plot and are separately shown in figure 8.

4. 3 . 4 .
Parallel electron force balance.The parallel electron momentum equation is given by equation (73) with electron instead of ion labels.In a plot of the terms analogous to the ion momentum plot figure 13(top) it turns out that most of the terms are very close to zero.We thus gather only the dominant terms in the electron momentum equation neglecting all terms proportional to the electron mass with m e = 0.This leaves the parallel force balance − T e ( b + b ⊥ ) • ∇n e + en e ( ( b + b

Figure 14 .
Figure 14.The sum of electron force balance and the parallel ion momentum equation (figure 13) neglecting small terms.The summed electric force is close to zero and drops out as does the resistivity.The error bars in the T i = 0 (left) plot become visible for η ⩽ 3 • 10 −5 while staying invisible for T i = Te (right).

Table 2 .
Simulation end times in units of Ω −1 i0 and in physical units reached after an equal amount of simulation time for all parameters.Simulations are run on 16 NVidia V100 GPUs.
the φ direction.The resolution is chosen as 192 cells in R and 336 cells in Z direction with 3 polynomial coefficients in each cell in both R and Z.The number ratio N R /N Z corresponds approximately to the aspect ratio of the simulation domain such that the grid cells are square in R-Z and the distance between neighbouring points is h ≈ 0.8ρ s .In φ we choose 32 planes.In total we thus have 576 • 1008 • 32 ≈ 2 • 10 7 grid points.Each simulation is run to roughly the same end time of 100 000 Ω −1 i0 with exact values displayed in table 2.

Table 3 .
Average computational time per Ω −1 i0 in seconds.A runtime of 6 s per Ω −1 i0 corresponds to a total simulation time of 7 days for 100 000 Ω −1 i0 .