Frédéric Schwander

Eric Serre

Hugo Bufferand

Guido Ciraolo

Philippe Ghendrih

Global fluid simulations of edge plasma turbulence in tokamaks: a review

Keywords: magnetic fusion, fluid modelling, numerical discretization, turbulence simulation

Introduction

Research in magnetic confinement fusion plasmas explores the possibility of producing power by using fusion in deuterium-tritium plasmas heated to temperatures of up to 10 7 -10 8 K, and confined by magnetic field in machines of toroidal shape known as tokamaks. ITER will be the largest machine of this type to date, and many challenging issues remain to solve before its operation at full performance. Power exhaust is certainly one of the most challenging among these, as fusion at ignition requires very high power in the core while the technological constraints of the material on the plasma-facing components impose to maintain manageable heat power fluxes to guarantee reactor safety [START_REF] Loarte | Power and particle control[END_REF][START_REF] Pitts | Physics basis for the first iter tungsten divertor[END_REF]. Controlling and predicting these heat fluxes remain uncertain, largely owing to an incomplete understanding of the mechanisms at play that result from the complex interplay of turbulence and transport processes in the plasma, losses at the wall, and the wealth of atomic and molecular interactions.

Heat exhaust properties are determined largely by the "edge plasma", which covers the region affected globally by the contact of the plasma with the wall. It encompasses both the outer boundary layer, in which the plasma is magnetically connected to the wall by open magnetic flux surfaces, referred in the literature as Scrape-Off Layer (SOL), and the outermost fraction of closed flux surfaces into the core, although its extension in the core is not rigorously defined and depends as much on the physics of interest as on authors (Fig. 1). Closed and open flux surfaces are separated by the Last Closed Flux Surface (LCFS) or separatrix. The thickness of the SOL, i.e. the decay length of density and temperature outside the LCFS, and over which the thermal power flows from the core to the wall, is a particular point of interest which depends to a large extent on the ratio of the turbulent transport across magnetic field lines to the very rapid transport along them [START_REF] Loarte | Power and particle control[END_REF].

The difficulty to get global experimental measurements in tokamaks makes reliable numerical simulations of the edge plasma desirable to complete the interpretation and understanding of tokamak discharges. The international fusion community has therefore been called upon to make a major modelling effort to construct efficient and reliable simulation tools, which are now at the heart of fusion research.

Despite the exponential growth of computer speed along with significant improvements in computer technology, ab initio simulation of high-performance plasma in realistic tokamak configurations (see some characteristics in Tab. 1) is still unreachable today. ITER, with its larger size and more intense magnetic field, makes the challenge of fullscale simulations even more demanding. Kinetic models, which require the knowledge of the particle distribution function lead to very costly simulations in a 6-dimensional phase space. Turbulent plasma fluctuations have long been recognized to have long wavelengths along the magnetic field, of the order of the parallel connection length ∝ q s R 0 (where q s is the safety factor, and R 0 the tokamak major radius defined thereafter), whilst turbulent fluctuations in the perpendicular direction reach scales of order on the ion or electron gyroradius ρ i,e depending on the type of turbulence [START_REF] Hennequin | Scaling laws of density fluctuations at high-k on Tore Supra[END_REF]. In extreme ranges, scales range from the electron gyro-radius (∼ 10 -5 m) and time scale (∼ 10 -11 s), to the length of the magnetic field line (∼ 10 2 m) and the diffusion time of the electrical current through the plasma (∼ 10 1 s). The difficulty is even larger at the edge of the tokamak where geometrical effects and plasma-wall interactions add complexity, even if some pioneering simulations start appearing in the fusion community (see for example in Ref. [START_REF] Chang | Gyrokinetic projection of the divertor heat-flux width from present tokamaks to ITER[END_REF][START_REF] Hager | Gyrokinetic study of collisional resonant magnetic perturbation RMP-driven plasma density and heat transport in tokamak edge plasma using a magnetohydrodynamic screened RMP field[END_REF][START_REF] Dif-Pradalier | Transport barrier onset and edge turbulence shortfall in fusion plasmas[END_REF][START_REF] Michels | Full-f electromagnetic gyrokinetic turbulence simulations of the edge and scrape-off layer of ASDEX Upgrade with GENE-X[END_REF]), but mainly addressing physical phenomena of fundamental interest for magnetic fusion.

The reduced dimensionality of fluid models compared to kinetic descriptions makes their computational cost significantly lower, and makes codes based on these models very popular to study transport and turbulence properties at the plasma edge. Therefore they remain a standard in the international community when considering tokamak relevant configurations with plasma-wall interactions in realistic geometries and parameter ranges. Various 2D/3D codes exist in the literature that are widely used in several international research programs to increase fundamental knowledge of tokamak physics, and also to support experimental interpretation. Depending on the range of scales resolved, these codes are usually classified as transport or turbulence codes, which could correspond quite well in CFD to Reynolds-Averaged-Navier-Stokes and Large-Eddy-Simulation codes, respectively. Transport codes are based on reduced models in which the turbulence has been smoothed by averaging. The resulting reduction in the computational cost makes it possible to account more precisely for plasma-wall interactions, which often induces coupling of codes because they are multi-physics. Thus, they are the current workhorse of the physicists closely connected to the operation and used to investigate and design optimal scenarios for reactors. In a non-exhaustive way, state-of-the-art transport codes are EDGE2D-EIRENE [START_REF] Simonini | Models and numerics in the multi-fluid 2-D edge plasma code EDGE2D/U. Contribut[END_REF][START_REF] Wiesen | EDGE2D/EIRENE code interface report[END_REF], EMC3-EIRENE [START_REF] Feng | Recent improvements in the emc3-eirene code[END_REF], SOLEDGE3X-EIRENE [START_REF] Bufferand | Implementation of multi-component Zhdanov closure in SOLEDGE3X[END_REF], SOLPS-ITER [START_REF] Wiesen | The new SOLPS-ITER code package[END_REF] or UEDGE-DEGAS2 [START_REF] Rognlien | A fully implicit, time dependent 2-D fluid code for modeling tokamak edge plasmas[END_REF]. Most of these simulations are 2D assuming the axisymmetry of the flow, restricted to the poloidal cross-section of the machine. 3D simulations can be however performed by EMC3-EIRENE (3D Edge Monte Carlo) and since very recently by SOLEDGE3X-EIRENE [START_REF] Genova | Modelling of tungsten contamination and screening in WEST plasma discharges[END_REF]. However, one of the main uncertainty with these codes is the transport across the magnetic field, which cannot rigorously be described by diffusion [START_REF] Ghendrih | Thermodynamical and microscopic properties of turbulent transport in the edge plasma[END_REF], and is thought to be turbulent [START_REF] D'ippolito | Convective transport by intermittent blob-filaments: Comparison of theory and experiment[END_REF]. Thus, since the fluctuations can be of similar spatial scales and magnitude to average profiles, significant efforts are focused on flux-driven fluid turbulence simulation codes. There have been numerous local studies of turbulence in the Scrape-Off Layer, mostly in two-dimensional reductions and in rather academic configurations with simplified geometries (see as examples HESEL [START_REF] Juul Rasmussen | Numerical modeling of the transition from low to high confinement in magnetically confined plasma[END_REF], TOKAM2D [START_REF] Sarazin | Intermittent particle transport in two-dimensional edge turbulence[END_REF]). However, the realization that turbulence in the edge and Scrape-off Layer is of the interchange type, which is ballooned and therefore asymmetric between Low Field Side (LFS) and High Field side (HFS), has pushed towards global, three-dimensional simulations of the edge. Thus, 3D codes were developed for the simulation of tokamak configurations. In a non-exhaustive way, we mention GBS [START_REF] Ricci | Simulation of plasma turbulence in scrape-off layer conditions: the GBS code, simulation results and code validation[END_REF], GDB [START_REF] Zhu | GDB: A global 3D two-fluid model of plasma turbulence and transport in the tokamak edge[END_REF], FELTOR [START_REF] Wiesenberger | Reproducibility, accuracy and performance of the Feltor code and library on parallel computer architectures[END_REF][START_REF] Wiesenberger | Effects of plasma resistivity in three-dimensional full-f gyro-fluid turbulence simulations[END_REF], GRILLIX [START_REF] Stegmeir | Global turbulence simulations of the tokamak edge region with GRILLIX[END_REF], STORM [START_REF] Easy | Three dimensional simulations of plasma filaments in the scrape off layer: A comparison with models of reduced dimensionality[END_REF] and HERMES [START_REF] B D Dudson | Hermes: global plasma edge fluid turbulence simulations[END_REF] (both built on the BOUT++ framework [START_REF] Dudson | BOUT++: A framework for parallel plasma fluid simulations[END_REF]), or SOLEDGE3X [START_REF] Tamain | The TOKAM3X code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF][START_REF] Bufferand | Implementation of multi-component Zhdanov closure in SOLEDGE3X[END_REF]. Here, all scales of the flow larger than the grid spacing are resolved in a self-consistent manner. These codes currently tackle realistic configurations but they can still only provide a rather limited description of plasma-wall interactions when compared to the transport codes previously cited.

As mentioned above, plasma turbulence involves a large range of space and time scales to resolve. In a fluid context, long parallel wavelengths are the result of the large electric conductivity of a fusion plasma, which favours slow variations of the electric potential in the parallel direction, and of large thermal conduction which smooths out temperature fluctuations along a magnetic field line. As a result, plasma turbulence in a tokamak is quasi-2D, requiring fine resolution in the plane across the magnetic field but only coarse resolution along it.

Even though some recent developments in these codes have proved the attractiveness of high-order finite-element methods like in SOLEDGE3X-HDG and FELTOR [START_REF] Giorgiani | A hybrid discontinuous Galerkin method for tokamak edge plasma simulations in global realistic geometry[END_REF][START_REF] Giorgiani | A new high-order fluid solver for tokamak edge plasma transport simulations based on a magnetic-field independent discretization[END_REF][START_REF] Scotto D'abusco | Core-edge 2D fluid modeling of full tokamak discharge with varying magnetic equilibrium: from WEST start-up to ramp-down[END_REF][START_REF] Wiesenberger | Reproducibility, accuracy and performance of the Feltor code and library on parallel computer architectures[END_REF][START_REF] Wiesenberger | Effects of plasma resistivity in three-dimensional full-f gyro-fluid turbulence simulations[END_REF], all the other codes are classically based on finite-difference / finite-volume numerical schemes (first and second-order accurate excepted in GBS which is fourthorder accurate) associated to an explicit or (semi)-implicit time discretization. Their major difference lies in the shape of the grid and the way the equations are discretized on it. They can consider Cartesian poloidal planes, assuming a rectangular wall cross-section like in GBS, or using an immersed boundary technique like in SOLEDGE3X, FELTOR and GRILLIX to handle general wall contours embedded in the rectangular computational domain, or they can use non-Cartesian coordinate systems with curvilinear coordinates like in BOUT++. A major complexity in the discretization is related to the strong anisotropy of the flow (k ⊥ k where k ⊥ is the wavenumber in the direction perpendicular to the magnetic field and k the wavenumber along the magnetic field) in edge plasma which can induce spurious numerical diffusion when dealing with the parallel operators [START_REF] Soler | A new conservative finite-difference scheme for anisotropic elliptic problems in bounded domain[END_REF][START_REF] Giorgiani | A high-order non field-aligned approach for the discretization of strongly anisotropic diffusion operators in magnetic fusion[END_REF]. The problem is made even more difficult because discretization schemes must simultaneously handle the complex magnetic geometry of the edge with the X-point(s), the vessel wall, and the divertor targets. The analysis recently provided in Stegmeir et al. [START_REF] Stegmeir | Analysis of locally-aligned and non-aligned discretisation schemes for reactor-scale tokamak edge turbulence simulations[END_REF] defines two large classes of methods used in the edge plasma community depending on the alignment or not of the discretization of parallel operators. The non-aligned methods allow an unequalled flexibility in discretizing geometries of any shape, but introduce spurious numerical diffusion. The capability of high-order schemes to limit the latter, as for example proposed in GBS, remains an open question, as shown by studies with contradictory conclusions on this subject [START_REF] Stegmeir | Analysis of locally-aligned and non-aligned discretisation schemes for reactor-scale tokamak edge turbulence simulations[END_REF][START_REF] Giorgiani | A high-order non field-aligned approach for the discretization of strongly anisotropic diffusion operators in magnetic fusion[END_REF]. With the locally-aligned methods used to discretize parallel direction [START_REF] Hariri | A flux-coordinate independent field-aligned approach to plasma turbulence simulations[END_REF], the anisotropy is computationally exploited by sparsifying the numerical grid along the resulting parallel direction, with a sparse toroidal mesh and a field-line map acting between few poloidal planes as for example in GRILLIX. However, with these methods a precise description of the geometry, and in a general way of the plasma-wall interactions, can be cumbersome.

All this suggests that edge plasma simulation in relevant tokamak configurations remains a computationally challenging problem. This paper proposes to focus the review on the numerical methods developed and implemented in the state-of-the-art turbulence codes, which are more demanding in terms of numerical performances than the transport codes due to the much larger number of degrees of freedom involved in the simulations. The paper is organized as follows: the geometrical configuration is introduced in Sec.2 ; the main assumptions and characteristics of the fluid models are detailed in Sec.3. Sec.4 discusses the main numerical techniques implemented in the codes. Sec.5 summarizes the work carried out around the verification and validation of these codes. Sec. 6 shows a selection of recent results illustrating the stat-of-the-art in edge plasma fluid simulations. Conclusion and discussions are provided in Sec. 7, and some perspectives in developments are outlined in Sec. 8 to complete the paper.

Geometrical configuration

The design of the tokamak wall and magnetic equilibrium essentially aims at satisfying constraints imposed by heat and particle exhaust issues.

All modern tokamaks use a rather complex magnetic equilibrium configuration, called divertor configuration. The objective is to separate as much as possible the plasma-solid interaction surface from the core where fusion reactions take place. The limiter configuration used in former tokamaks will not be discussed here (see for example a sketch in Ref. [START_REF] Giorgiani | A magnetic-field independent approach for strongly anisotropic equations arising plasma-edge transport simulations[END_REF] and simulation results in Refs. [START_REF] Tamain | TOKAM-3D: A 3D fluid code for transport and turbulence in the edge plasma of tokamaks[END_REF][START_REF] B D Dudson | Hermes: global plasma edge fluid turbulence simulations[END_REF]).

A sketch of the resulting geometry is shown on Fig. 1. The intersections between the flux surfaces and an arbitrary poloidal plane have the shape of the phase portrait of a harmonic oscillator, with a point of zero poloidal field appearing called the X-point. The volume of the plasma on the opposite side of the X-point with respect to the magnetic axis is called the divertor. It is composed of two distinct regions: a fraction of the so-called main SOL, which encompasses the field lines connected to the vessel wall and radially connected to the main confined plasma, and the so-called private-flux region, which is a micro-SOL radially connected to the main SOL but not to the confined plasma.

Let R 0 be the major radius, i.e. distance from the magnetic axis, and (R, Z, ϕ) a fixed polar coordinates system, Fig. 1. The magnetic equilibrium is assumed to be toroidally symmetric and to encompass both closed and open flux surfaces, even if some very recent codes capabilities allow to deal with 3D magnetic equilibrium as for example in Ref. [START_REF] Shanahan | Fluid simulations of plasma filaments in stellarator geometries with BST-ING[END_REF][START_REF] Coelho | Validation of GBS plasma turbulence simulation of the TJ-K stellarator[END_REF]. The 2D equilibrium magnetic field is defined as B = B b, ( b = B/B is the unit vector along the magnetic field) that can be expressed in a general way as [START_REF] Wesson | Tokamaks. International Series of Monogr[END_REF]:

B = F ∇ϕ + ∇ψ × ∇ϕ (1)
where ϕ denotes the toroidal angle, F is a toroidal flux and ψ(R, Z) the poloidal flux function. According to Eq. 1, the iso-ψ surfaces are tangent to the magnetic field and ψ labels flux surfaces (one value for each flux surface). It is thus natural to define a curvilinear system of coordinates denoted (ψ, θ , ϕ), where θ denotes a curvilinear abscissa along the poloidal direction in the (R, Z) plan that defines the poloidal plane, i.e. along iso-ψ surfaces and orthogonal to ∇ϕ (Fig. 1). In the base ( e ψ , e θ , e ϕ ) associated to (ψ, θ , ϕ) (Fig. 1), the magnetic field writes:

B = B p e θ | e θ | + B t e ϕ | e ϕ | (2)
where (B p , B t ) denote the poloidal and the toroidal components of the magnetic field, respectively. In a tokamak, the magnetic field can be of strong intensity (Table 1) and with |B p | |B t |. This leads to define a privileged flow direction In the core and edge the magnetic field lines are either closed, fill a poloidal flux-surface or volume of space (extracted from [START_REF] Held | Full-f gyro-fluid modelling of the tokamak edge and scrape-off layer[END_REF]). (right)

Coordinates system and geometrical parameters R 0 be the major radius, a the minor radius on the magnetic axis and (R, Z, ϕ) a fixed cylindrical coordinates system. (ψ, θ , ϕ) curvilinear system of coordinates (extracted and modified from [START_REF] Lütjens | The CHEASE code for toroidal MHD equilibria[END_REF]).

along which the governing equations are projected using differential operators ∇ = b • ∇ and ∇ ⊥ = ∇b • ∇ to define gradients along the parallel and perpendicular direction, respectively (See Sec. 3).

To end this section, the characteristic parameters of some of the major tokamaks in operation, as well as those of ITER currently under construction, are summarized in Table 1. The Larmor radius ρ s = √ 2T i m i /(eB 0 ) is given for each machine at a temperature characteristically found at the separatrix taken as 100eV. It fixes the order of the smallest theoretical scale to be resolved in the flow at the edge of a tokamak. The ratio a/ρ s , evaluated at the same temperature, is given as a measure of scale separation between turbulence scale at the separatrix and machine size. This ratio, which grows with magnetic field intensity and minor radius illustrates the challenge of full-scale simulation of ITER.

Fluid-based model

2D/3D fluid conservation equations for electrons and ions remain the standard in all edge plasma codes of the literature. We attract here the reader's attention, on the fact that depending on the assumptions made and the degree of simplification introduced in the model, the terms in the equations can slightly change from a code to another one. It is therefore impossible to be exact and exhaustive in the formulation when introducing fluid-based model.

This section aims to introduce the drift-reduced Braginskii framework (Sec. 3.1) in which a generic set of fluid equations is introduced (Sec. 3.2) together with boundary conditions (Sec. 3.3). The reader will have to refer to the various works in the literature in order to know the fluid model exactly resolved. In the plasma edge modelling, neutrals (Sec. 3.4) are mandatory to estimate the sources of particles in the plasma equations due to charge exchange, ionization and recombination reactions between the plasma and one atomic neutral species. Indeed, neutral particles play a major role in establishing the density and temperature profiles through recycling, i.e. neutralisation of the plasma at the wall and re-ionization of the neutral gas.

While beyond the scope of this review, it is however relevant to close this Section by briefly mentioning gyrofluid models as well as specific closure for multi-species plasma when taking into account impurities (Sec. 3.5).

The drift-reduced Braginskii framework

Derived as the successive moments of the kinetic equations, most of the models are based on Braginskii closures, which use the Chapman-Enskog approximation with two Laguerre-Sonine polynomials [START_REF] Braginskii | Transport processes in a plasma[END_REF].

5.2×10

Table 1: Main characteristics of some of the major tokamaks in operation together with ITER currently under construction. The Larmor radius ρ s = √ 2T i m i /(eB 0 ) is given for each machine at a temperature characteristically found at the separatrix taken as 100eV.

The assumption of drift-ordering, introduced for fluid plasmas in [START_REF] Zeiler | Nonlinear reduced Braginskii equations with ion thermal dynamics in toroidal plasma[END_REF] for example, is frequently used in the Braginskii fluid models to make the computations more accessible (see for example in [START_REF] Madsen | Collisional transport across the magnetic field in driftfluid models[END_REF][START_REF] Tamain | The TOKAM3X code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF][START_REF] B D Dudson | Hermes: global plasma edge fluid turbulence simulations[END_REF][START_REF] Stegmeir | GRILLIX: a 3D turbulence code based on the flux-coordinate independent approach[END_REF][START_REF] Giacomin | The GBS code for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary[END_REF]), although its validity outside the SOL is widely discussed in the literature [START_REF] Gath | Consistency in drift-ordered fluid equations[END_REF][START_REF] Frei | Local gyrokinetic collisional theory of the ion-temperature gradient mode[END_REF]. This corresponds to a model reduction assuming the turbulence is characteristically low frequency and long wavelength in nature leading to a strong scale separation between the characteristic plasma frequency ω of the turbulence compared to the ion cyclotronic frequency ω c with which the particles spiral along magnetic field lines (gyromotion). This scale separation allows to treat the fluid-momentum equation in a singular perturbation expansion around the small parameter = ω/ω c 1 that makes it possible to solve explicitly for the perpendicular components to the magnetic field of the fluid velocity, the so-called drift velocities. Then the total velocity writes,

v = v b + v ⊥ (3)
where the perpendicular components of the velocity can be explicitly described in terms of drifts:

v ⊥ = v E + v + v p (4)
Owing to the large magnetic field present in a tokamak, the fluid velocity in the plane orthogonal to the magnetic field is given by a quasi-static balance between Lorentz force, pressure gradient and electromotive force. The following expressions for the perpendicular drift velocities can be derived from this balance (see for example in Refs. [START_REF] Stegmeir | GRILLIX: a 3D turbulence code based on the flux-coordinate independent approach[END_REF][START_REF] Bufferand | Progress in edge plasma turbulence modelling-hierarchy of models from 2D transport application to 3D fluid simulations in realistic tokamak geometry[END_REF]):

• v E := ( B× ∇φ)/B 2 is the perpendicular electric drift with φ the electric potential (identical for ions and electrons),

• v := ( B × ∇p)/(nqB 2 ) is the diamagnetic velocity with q and n are the charge and the particle density, respectively,

• v p = -1 n ∂ t ω + ∇ • v ⊗ ω is the second-order polarisation velocity where ω = m qB 2 n ∇ ⊥ φ + 1 q ∇ ⊥ p (with the vorticity Ω = ∇ • ω).
It is linked to inertia effects and proportional to the mass and therefore much larger for the ions compared with the electrons.

If the first-order drift velocity v E is taken into account in all models, the two others can be considered or not, and possibly treated differently depending on the assumptions made. For example, the diamagnetic velocity can be sometimes expressed as a magnetic drift as v ∇B = ± (+ for hydrogen ions,for electrons), because it is more easily implemented in terms of fluxes through cell faces than the standard form, and hence it is more suitable the conservative numerical schemes [START_REF] Tamain | The TOKAM3X code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF][START_REF] Dudson | BOUT++: A framework for parallel plasma fluid simulations[END_REF]. The second order polarisation velocity linked to inertia effects, can be totally neglected, or for the electrons only, when electronic inertia is neglected [START_REF] Madsen | Collisional transport across the magnetic field in driftfluid models[END_REF][START_REF] Tamain | The TOKAM3X code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF][START_REF] Stegmeir | GRILLIX: a 3D turbulence code based on the flux-coordinate independent approach[END_REF].

Thus, the perpendicular component of the velocity being no longer an independent variable of the problem the resolution of each equation is simplified.

A generic system of fluid equations

The paper focuses on the drift-reduced Braginskii equations, which are resolved by all codes (except FELTOR which solves a set of gyrofluid equations). For all charge species (typically ions and electrons), they describe the time evolution of particle density n, velocity v , and temperature T .

We only present here a generic system of conservation equations for a simple quasi-neutral (n i = n e ) hydrogenic plasma (species charge number Z = 1), assuming both gases as ideal gases (p = nT ) and derived from Ref. [START_REF] Bufferand | Progress in edge plasma turbulence modelling-hierarchy of models from 2D transport application to 3D fluid simulations in realistic tokamak geometry[END_REF]:

∂ t n + ∇ • n v = ∇.(D ⊥ ∇ ⊥ n/n) + S n (5)
∂ t (mnv ) + ∇ • mnv v = -∇ p + nqE + R (6)
+ ∇ • ν ∇ v b + mnν ⊥ ∇ ⊥ v + S v ∂ t E + ∇ • E v = -∇ • pv b + nqv E + R v (7)
+ ∇ • ν v ∇ v b + mnν ⊥ v ∇ ⊥ v + ∇ • κ ∇ T b + χ ⊥ n ∇ ⊥ T + S E
where v is defined by Eq.3, v D = -D ⊥ ∇ ⊥ n/n stands for an anomalous diffusive transport in the density equation, m denotes the species mass. E = 3 2 nT + 1 2 mnv 2 is the species total energy. Parallel viscosity ν , parallel heat conductivity κ take values computed by Braginskii closure. S n , S v and S E are the source terms, and R a friction term between species. Finally, cross-field diffusivity coefficients emulating turbulent transport (D ⊥ , ν ⊥ and χ ⊥ ) are constant, taking classical values of the order of ∼ 10 -2 m 2 • s -1 . As a comparison, for transport simulations classical values are of the order of ∼ 1m 2 •s -1 and can be non constant, determined from experimental profiles at mid-plane or be determined more self-consistently as in RANS eddy viscosity models in CFD [START_REF] Baschetti | the EUROfusion MST1 team, and the TCV team. Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF][START_REF] Baschetti | Self-consistent cross-field transport model for core and edge plasma transport[END_REF][START_REF] Coosemans | A self-consistent mean-field model for turbulent particle and heat transport in 2D interchange-dominated electrostatic E×B turbulence in a sheath-limited scrape-off layer[END_REF][START_REF] Dekeyser | A self-consistent κ-model for anomalous transport due to electrostatic, interchange-dominated E×B drift turbulence in the scrape-off layer and implementation in SOLPS-ITER[END_REF]. For electrons, only the energy equation is solved, the electron density being computed assuming quasi-neutrality and electron velocity being computed assuming ambipolarity.

The drift approximation has the specificity of making the treatment of perpendicular dynamics very similar to that of two-dimensional Navier-Stokes, where the electric potential and ion or electron pressure play a role close to that of the flux function. Not unlike the Navier-Stokes equation, the electric potential needs to be determined through its relation to a vorticity-like quantity. The equation determining the time-evolution of this vorticity is usually obtained through the quasi-neutrality condition.

∇ • j = 0 ( 8
)
Assuming a single ion population with single positive electric charge results in:

∇ • [ne( v p,i + v * ,i ) -ne v * ,e ] + ∇ • j = 0 ( 9
)
where e is the electron charge. Note that the E × B drift does not appear in this equation, since it does not carry any current, and that the electron polarization drift is omitted here owing to the smallness of the electron/ion mass ratio m e /m i 1. Moreover, a diffusion current j diff can be also added in the codes (as for example in SOLEDGE3X) for numerical stability that expresses as j diff = ζ ∇ ⊥ Ω. This equation becomes an equation for the electric potential using

v p,i = B eB 2 × m i d( v E + v * ,i ) dt ( 10
)
where the total time derivative is that of the ion fluid, which allows to rewrite an equation for the vorticity

Ω = ∇ • B eB 2 × nm i ( v E + v * ,i ) = -∇ • m i [n∇ ⊥ φ + ∇ ⊥ p i /(e)] eB 2 (11)
where ∇ ⊥ is the projection of the gradient in the plane orthogonal to B. This equation might found in the literature written as:

∂Ω ∂t + ( v E + v * ,i + v ,i ) • ∇Ω + ne( v * ,i -v * ,e ) + ∇ • ( j b) = 0 ( 12
)
Several forms of the equation exist, especially in the treatment of the vorticity flux, in which modelling is often liberal in the commutation of spatial differential operators and projection along the magnetic field. This equation also illustrates the main difference between electrostatic and electromagnetic models described in the recent literature.

In electrostatic framework, parallel current is related to electric potential and electron pressure (neglecting electron inertia), by projecting the resistive version of Ohm's law along the direction of the electric field, to obtain:

η j = -∇ φ + ∇ p e ne + 0.71∇ T e e ( 13
)
where η is the plasma resistivity. This relation can then be used to substitute j in Eq. 12. One must recall that this equation is used to determine the electric potential through the relation between Ω and φ. Let's mention that considering electron inertia as presented in the electrostatic model of GBS in [START_REF] Paruta | Simulation of plasma turbulence in the periphery of diverted tokamak by using the GBS code[END_REF] for example, makes j an independent variable through nv ,i and nv ,e .

In electromagnetic framework [START_REF] Zholobenko | Electric field and turbulence in global Braginskii simulations across the ASDEX Upgrade edge and scrape-off layer[END_REF][START_REF] Düll | Introducing electromagnetic effect in Soledge3X. Contrib. to[END_REF], the parallel current is defined differently since it becomes linked to the potential A through Ampère law which yields:

µ 0 j = -∇ 2 ⊥ A (14)
An evolution equation for the parallel potential is obtained by using the relevant definition of the parallel electric field E = -∇ φ -∂ t A , which in turn yields the electromagnetic form of parallel Ohm's law Eq.13, still ignoring electron inertia:

η j = - ∂A ∂t -∇ φ + ∇ p e ne + 0.71∇ T e e (15)
Thus, the electric and magnetic potential are coupled and must be resolved together (see for example in [START_REF] Zholobenko | Electric field and turbulence in global Braginskii simulations across the ASDEX Upgrade edge and scrape-off layer[END_REF][START_REF] Giacomin | The GBS code for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary[END_REF]) by using Eqs. 12 and 15, which yields a system for φ and A which describes shear-Alfvén waves. Fluctuations of the potential result in fluctuations in the poloidal magnetic field:

B = ∇ × Ã b (16)
where the unit vector b is obtained from the equilibrium magnetic field.

Neutrals equations

Neutral particles and their interaction with the plasma are central to issues of heat and particle exhaust for magnetic fusion. The momentum and energy volumetric losses they are responsible for when recycling back into the plasma are the main drive for accessing detached conditions which will be mandatory in ITER and DEMO to keep heat fluxes at the target plates below the engineering limit. Increasing experimental evidence also demonstrates an impact of density regimes on transverse plasma transport, which might in turn influence the access and stability of detachment.

There is a wide variety of models based on either fluids or kinetic equations, depending on the expected accuracy and computational cost. However, while edge transport codes have for long included neutrals dynamics, edge turbulence modelling has only recently started addressing the topic and many open questions subsist in terms of numerical implementation and physics [START_REF] Wersal | A first-principles self-consistent model of plasma turbulence and kinetic neutral dynamics in the tokamak scrape-off layer[END_REF][START_REF] Fan | Self-consistent coupling of the three-dimensional fluid turbulence code TOKAM3X and the kinetic neutrals code EIRENE[END_REF][START_REF] Thrysøe | Plasma particle sources due to interactions with neutrals in a turbulent scrape-off layer of a toroidally confined plasma[END_REF].

The kinetic approach is well justified because the neutral mean-free path in the tokamak edge can be larger than the machine size, but it remains computationally expensive. It is therefore only widely used in transport simulations [START_REF] Bufferand | Near wall plasma simulation using penalization technique with the transport code SOLEDG2D-EIRENE[END_REF][START_REF] Wiesen | The new SOLPS-ITER code package[END_REF][START_REF] A V Chankin | EDGE2D-EIRENE simulations of the influence of isotope effects and anomalous transport coefficients on near scrape-off layer radial electric field[END_REF] through code coupling with EIRENE [START_REF] Reiter | The EIRENE and B2-EIRENE codes[END_REF] or DEGAS2 [START_REF] Russell | A reduced model of neutral-plasma interactions in the edge and scrape-off-layer: Verification comparisons with kinetic Monte Carlo simulations[END_REF] based on a Monte Carlo algorithm. In turbulence simulations, the first attempt to couple a 3D turbulence code with EIRENE was carried out in Ref. [START_REF] Fan | Self-consistent coupling of the three-dimensional fluid turbulence code TOKAM3X and the kinetic neutrals code EIRENE[END_REF]. The model was not extensively exploited since due to its high computational cost. In GBS, another approach has been chosen based on the implementation of a kinetic equation directly embedded into the plasma solver avoiding code coupling [START_REF] Wersal | A first-principles self-consistent model of plasma turbulence and kinetic neutral dynamics in the tokamak scrape-off layer[END_REF]. This allows 3D turbulence simulations in realistic configurations as in TCV [START_REF] Giacomin | The GBS code for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary[END_REF] including kinetic neutrals.

The fluid approach remains however the standard one in 3D turbulence simulations with varying degrees of refinement. Generally, a simple fluid-diffusive model, the neutrals diffusing homogeneously, is used based on the neutrals density conservation. Despite its simplicity, this modelling provides a good approximation for the source of plasma generated by recycling as well as power losses by radiation in the divertor as shown in [START_REF] Zholobenko | Electric field and turbulence in global Braginskii simulations across the ASDEX Upgrade edge and scrape-off layer[END_REF][START_REF] Bufferand | Progress in edge plasma turbulence modelling-hierarchy of models from 2D transport application to 3D fluid simulations in realistic tokamak geometry[END_REF]. In Ref. [START_REF] Leddy | Simulation of the interaction between plasma turbulence and neutrals in linear devices[END_REF], the model is extended to consider parallel momentum and neutrals pressure evolution.

Boundary conditions

Boundary conditions for fluid models are required at the tokamak wall, but generally also at the core edge boundary, as most codes do not extend to the centre of the tokamak. It is necessary to distinguish the parallel and perpendicular directions to the magnetic field lines which are not governed by the same physics as mentioned above. They are generally implemented as for the fluid equations, with varying degree of sophistication. The reader is referred to go in the detail of each paper in order to get the right expressions which have been implemented. As for the fluid equations we can however provide a set of generic boundary conditions, here inspired from [START_REF] Bufferand | Progress in edge plasma turbulence modelling-hierarchy of models from 2D transport application to 3D fluid simulations in realistic tokamak geometry[END_REF] and as for the equations neglecting electron inertia.

-In the perpendicular direction, zero Neumann boundary conditions for all plasma variables, i.e. ∂ ⊥ (.) = 0 is imposed both at the wall and the core edge boundary, when the models do not include the entire plasma core. Let's mention that this core edge boundary condition does not have a clear physical analogue.

-In the parallel direction, boundary conditions are derived from the generalized Bohm-Chodura sheath boundary conditions [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF]. They model the physics of the sheath located next to the limiter wall, and in which many assumptions used to derive the fluid models (quasi-neutrality, drift-ordering) are no longer valid.

• Outgoing velocity normal to the wall larger that parallel sound speed normal to the wall:

| v • n wall | ≥ |c s b • n wall | (17)
This property guarantees that the total plasma velocity is oriented outward

• Sheath transmission factor between energy and particle fluxes:

φ E,se = γT φ n,se (18)
Where for each species, φ E,se is the total energy flux at sheath entrance, φ n,se is the particle flux at sheath entrance and γ is the sheath transmission factor about 2.5 for ions and 4.5 for electrons.

• Total plasma current on the wall is given by

j wall = 1 -exp Λ - φ T e φ n,se (19)
Where the ions saturation current is computed from ions particle fluxes φ n,se and where Λ denotes the normalised potential drop in the sheath Λ ∼ 3.

In electromagnetic framework the boundary conditions above do not change. An additional condition on A is provided, which generally corresponds to A = 0 at the magnetic pre-sheath entrance [START_REF] Giacomin | The GBS code for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary[END_REF][START_REF] Düll | Introducing electromagnetic effect in Soledge3X. Contrib. to[END_REF].

Beyond the Braginskii closure

To conclude this Section, we briefly mention two closures beyond the standard Braginskii closure, i.e. gyrofluid models as well as specific closure for multi-species plasma which are beginning to emerge in the literature.

The drift-fluid models described above invoke a collisional closure and thus miss reactor relevant kinetic effects, such as small amplitude fluctuations at the gyro-radius scale and low collisionality effects. For decades indeed, gyrofluid models have been developed by considering a finite number of velocity space moments of the GK equation to extend the validity of Braginskii-like fluid by introducing some collisionless kinetic effects. This gyrofluid models appear as the most balanced compromise between fidelity and computing cost to be used as the high-end of the multi-fidelity hierarchy of models. Although large similarities exist in the structure of equations, such a model is only exploited in FELTOR [START_REF] Wiesenberger | Reproducibility, accuracy and performance of the Feltor code and library on parallel computer architectures[END_REF] and BOUT++ [START_REF] Ma | Global kinetic ballooning mode simulations in bout++[END_REF] to our knowledge due to the added theoretical and numerical difficulties. The interested reader will find very relevant references and analysis in the recent works in Refs. [START_REF] Wiesenberger | Long-wavelength closures for collisional and neutral interaction terms in gyrofluid models[END_REF][START_REF] Frei | Local gyrokinetic collisional theory of the ion-temperature gradient mode[END_REF].

A fusion plasma is inherently composed of multiple ion species. Due to the fusion reaction itself (Deuterium, Tritium and Helium) but also from the plasma-wall interactions. Indeed, intrinsic impurity ions, like tungsten, carbon or oxygen, are present, originating from the wall-material through erosion or resulting from an incomplete vacuum. Also, other ion species are brought deliberately into the plasma, for example through neutral beam injection or gas puffing for heating, fuelling or diagnostic purposes like Nitrogen, Neon or Argon. Thus, an accurate description of a fusion edge plasma therefore must include capabilities to treat multi-species physics. However in this case the Braginskii closure fails and such limitation is a blocking point.

From a modeling perspective there is then no difference between a main ion species and an impurity species. In gyro-fluid models mentioned above, this extension is straightforward since the equations for these models are by construction entirely symmetric for each species. However, the inter-species collisional terms in a gyro-fluid model are intricate to work out in practice. In contrast, in drift-fluid models the collisional terms are well-understood but the multi-species extension leads to two problematic issues [START_REF] Poulsen | Collisional multispecies drift fluid model[END_REF][START_REF] Raghunathan | Multi-temperature generalized Zhdanov closure for scrape-off layer/edge applications[END_REF]: (i) the equations for electric potential and all ion pressures are coupled in the time-derivative such that numerically a non-linear mass matrix has to be solved in each time-step. In an implicit time-stepper this leads to a nested inversion problem which is non-trivial and computationally expensive to solve. (ii) the equations are not symmetric in the sense that multiple equal ion species equations do not describe the same as the single ion species model. Even though a high order effect, in turbulence simulations this may be problematic. The Zhdanov closure for multi-component plasma proposes an alternative to Braginskii closure and has been recently implemented in transport codes like SOLPS-ITER [START_REF] Sytova | Derivation of the friction and thermal force for SOLPS-ITER multicomponent plasma modeling[END_REF] or in turbulence codes as SOLEDGE3X [START_REF] Bufferand | Implementation of multi-component Zhdanov closure in SOLEDGE3X[END_REF] or GBS [START_REF] Mancini | Self-consistent multi-component simulation of plasma turbulence and neutrals in detached conditions[END_REF]. Results have show significant improvements in the prediction for turbulent transport in different machines.

Numerical schemes

Fluid modelling of edge plasmas is primarily based on conservation equations (Sec.3.2), and thus shares many features of fluid modelling in other contexts. It carries however specific features which make these equations specially challenging to resolve and stressful for the numerical methods. Among them, the strong magnetization of the plasma that leads to a marked anisotropy, the complex geometry of both tokamak wall and magnetic equilibrium with singularities, the high temperatures and low density, the presence of charged species that make them sensitive to electromagnetic forces, and finally the interaction of the plasma with the wall which require the development of dedicated numerical algorithms. This section reviews the key features of the latter, implemented in state-of-art simulation codes of the literature, for both the spatial and time discretizations.

Spatial discretizations

The specific features of the fluid equations mentioned above influences the spatial discretization strategies, both in the meshing and in the development of the discretization schemes for the parallel differential operators. In particular, turbulence in a tokamak being quasi-2D it requires fine resolution in the plane across the magnetic field but only coarse resolution along it. Moreover, in diverted geometry (Sec. 2), the poloidal angle is no longer a valid concept owing to the singularity resulting from the presence of the X-point (saddle-point for the poloidal flux function) which makes the LCFS non-differentiable.

A recent analysis of locally-aligned and non-aligned discretisation schemes for reactor-scale tokamak edge turbulence simulations can be found in Ref. [START_REF] Stegmeir | Analysis of locally-aligned and non-aligned discretisation schemes for reactor-scale tokamak edge turbulence simulations[END_REF]. It appears to us that the most recent edge modelling codes ought to be differentiated by their mesh construction as well by the manner in which they handle parallel gradients (aligned or not). We must comment that the classification is not standard, as the issues of meshing and spatial discretization are often entangled, but is seems to us that this classification gives a clear overview of choices made in existing codes. Summary of this classification for some of the edge modelling codes is given in Table 2, which illustrates that no consensus on the optimal combination, if it exists, has been found.

Code

Mesh design Parallel discretization FCIM FCDM Aligned Non-aligned BOUT FELTOR GBS GDB GRILLIX SOLEDGE3X SOLPS Table 2: Summary of mesh designs and discretization strategies of parallel gradients for fluid codes. FCIM stands for Flux-Coordinate Independent Mesh and FCDM for Flux-Coordinate Dependent Mesh.

Meshing strategies

We separate codes in two classes defined as using Flux-Coordinate Dependent Mesh (FCDM) or Flux-Coordinate Independent Mesh (FCIM), illustrated in the following.

-Flux-Coordinate Dependent (FCD). Here, in order to maintain flux-aligment, local spatial variables are constructed so that all but one of the contravariant basis vectors are tangent to the magnetic surface, following strategies already adopted in early 2D transport codes [START_REF] Rognlien | A fully implicit, time dependent 2-D fluid code for modeling tokamak edge plasmas[END_REF][START_REF] Knoll | NEWEDGE: a 2D fully implicit edge plasma fluid code for advanced physics and complex geometries[END_REF]. This leads to a multi-block structured approach for the meshing of the poloidal cross-section due to the presence of an X-point in the computational domain, which makes the radial and poloidal coordinates singular (see for example codes as SOLEDGE3X or SOLPS). Typically, the mesh is decomposed in six blocks, covering the core, the SOL and the private region, as illustrated in Figure 2. This strategy can be coined as Flux-Coordinate Dependent Mesh (FCDM), as one of the spatial coordinates can be seen as a function of the poloidal flux function.

This strategy presents the advantage that the parallel gradient operates only on a 1D-lattice in 2D, or 2D lattice in 3D, on each of the blocks constituting the global mesh. It carries the disadvantage that the structure and conforming nature of meshes used result in computational cells that are coarser close to the X-point than away from it, especially at the separatrix, Fig. 3 (left). [START_REF] Wiesenberger | Streamline integration as a method for structured grid generation in X-point geometry[END_REF]). The X-point has a different connectivity with respect to the other cells and a coarser resolution around it. (FCIM) (right) (extracted from Ref. [START_REF] Stegmeir | GRILLIX: A 3D turbulence code for magnetic fusion devices based on a field line map[END_REF]). The Cartesian mesh allows a homogeneous resolution and the X-point is no longer a singularity for the mesh.

Generating the mesh in this framework is a task specific to fusion plasmas, especially for diverted plasmas with one or several X-points. Most codes rely on dedicated grid generation strategies in the 2D poloidal cross-section [START_REF] Wiesenberger | Streamline integration as a method for structured grid generation in X-point geometry[END_REF] and tools [START_REF] Marchand | CARRE: a quasi-orthogonal mesh generator for 2D edge plasma modelling[END_REF][START_REF] Garcia | INGRID: An interactive grid generator for 2D edge plasma modeling[END_REF], and this 2D mesh is simply extruded in the toroidal dimension to obtain the final 3D mesh. Though advantageous, the latter meshing strategies are not capable of taking advantage of the quasi-2D nature of the tokamak plasma turbulence : typically, in a conventional 3D curvilinear mesh, the quasi-aligned character of turbulent fluctuations imposes that number of grid points N ϕ is linearly related to the number of grid points in the poloidal direction N θ by

N ϕ ≈ N θ q (20)
where q is the safety factor (typically q ∼ 1 -10). Combined with the properties that the smallest turbulent structures are isotropic in the poloidal plane, and that the same structures are of size comparable to the ion gyroradius, this induces the unfavourable scaling of total number of points N grid with normalized gyroradius:

N grid ∝ ρ 3 * q (21)
-Field-Aligned Mesh (FAM). Early attempts had been made to circumvent this difficulty in local, gyrokinetic turbulence simulations [START_REF] Beer | Field-aligned coordinates for nonlinear simulations of tokamak turbulence[END_REF], in which a curvilinear mesh is constructed so as to be aligned with the magnetic field, and the parallel gradient thus operates on a 1D-lattice in three-dimensional simulations. This indeed permits the use of a very coarse mesh in the parallel direction. This kind of approach does not translate easily to global simulations since magnetic field lines do not systematically close: a strictly aligned mesh therefore has to be non-conforming, rendering a conventional finite-difference treatment more involved. To our knowledge, BOUT is the only code that has pursued this strategy [START_REF] Umansky | Status and verification of edge plasma turbulence code BOUT[END_REF], in which the 2D mesh in the poloidal section is constructed by extrusion following magnetic field lines. To our knowledge, no current fluid turbulence code uses this method.

-Flux-Coordinate-Independent Mesh (FCIM). Finally, and perhaps surprisingly, the last strategy appearing is one where the mesh is constructed in the poloidal cross-section without any regard for the poloidal flux function of magnetic field. The prototypes for this strategy are FELTOR, GBS and GRILLIX on the one hand for turbulence codes, which use a straightforward 3D Cartesian mesh using the (R, ϕ, Z) coordinates (Fig. 3 (right)), and the high-order, finiteelement code SOLEDGE3X-HDG [START_REF] Giorgiani | A hybrid discontinuous Galerkin method for tokamak edge plasma simulations in global realistic geometry[END_REF] for transport codes, where the volume inside the tokamak walls is meshed using a classical triangulation.

4.1.2. Field-aligned discretizations for parallel gradients ∇ = b. ∇ All of the above codes, except SOLEDGE3X-HDG and FELTOR, have been constructed in a finite-difference setting and, although some care has been put in the discretization of advection terms and diffusion terms, standard finite-difference discretizations for spatial operators are used in FCDM, FCIM and FAM codes.

A change in paradigm on these discretizations was incited in [START_REF] Scott | Shifted metric procedure for flux tube treatments of toroidal geometry: Avoiding grid deformation[END_REF], in which it was proposed that parallel gradients be computed by finite-differences, not using node values, but using values along the field line obtained via a perpendicular interpolation. This forms the basis of the shifted-metric procedure, in which the resolution in one direction (typically the toroidal one) can be coarsened without loss of information owing to slow variations in the parallel direction as is the case in FAM, and without the constraint that the mesh follow the magnetic field. This technique effectively exploits efficiently the quasi-2D nature of turbulence of tokamaks and avoids systematic over-resolution in the parallel direction. It has also led to the development of an even more aggressive technique, in which the mesh is constructed without any consideration for the magnetic field. The coarsening of the mesh in one direction is there made possible by a specific evaluation of parallel gradients, which is performed by differentiating the function along the field line with the help of local interpolation. This is the field-aligned discretization proposed in [START_REF] Hariri | A flux-coordinate independent field-aligned approach to plasma turbulence simulations[END_REF], which asserted the fact that the quasi-aligned nature of turbulent fluctuations could be exploited to coarsen resolution in one direction even in FCIM. The use of this possibility is one of the defining features of GRILLIX and is also used by BOUT++ [START_REF] Shanahan | Fluid simulations of plasma filaments in stellarator geometries with BST-ING[END_REF]. A schematic illustration of this method is given in Fig. 4 where the parallel operators are discretised via field line map, i.e. field line tracing and interpolation. We note that recent generalization to finite-volumes has been proposed in [START_REF] Wiesenberger | A finite volume flux coordinate independent approach[END_REF].

As in the original paper of Hariri and Ottaviani [START_REF] Hariri | A flux-coordinate independent field-aligned approach to plasma turbulence simulations[END_REF], Flux-Coordinate-Independent is related to the construction of the mesh, while aligned or non-aligned qualifies the discretization of parallel gradients. This view is further reinforced by the fact that, in its most recent version, GBS uses the FCIM approach, with a rectangular Cartesian mesh in the poloidal cross-section, but without aligned discretization [START_REF] Boinnard | Plasma turbulence simulations in a diverted tokamak with applied resonant magnetic perturbations[END_REF].

Discretization of the wall -boundary conditions

Edge modelling is concerned with the interaction between the plasma and the wall. At the lowest modelling level, the wall enters as a boundary condition. The wall itself is crucial in determining the particle and energy circulation pattern inside the plasma, and control of the heat power deposition on the first wall is one of critical technological issues on the way towards power production through magnetically confined fusion. Though they rely on first-principle (Right) Scheme for toroidally staggered FCIM approach. Grid is indicated with crosses and dual grid with circles. Extracted from [START_REF] Stegmeir | GRILLIX: a 3D turbulence code based on the flux-coordinate independent approach[END_REF].

approaches, plasma edge modelling codes have therefore moved towards an increasingly detailed account of the geometry of the wall, and of the physical phenomena occurring at the interface between the plasma and the first wall. In order of priority, the first issue that an edge modelling ought to treat is the particle and heat flux into surfaces impinged by the magnetic field, such as limiters or divertor plates, which are meant to contribute the most to power exhaust. As presented in Sec. 3.4, these wall elements are modelled through use of outflow boundary conditions on the plasma flow, and must be included in the computational domain for any relevant SOL modelling. Turbulence codes classically focus more on turbulence around the LCFS and flux impacting divertor plates or limiters. This justifies the use of numerical boundary conditions on a flux surface outside and sufficiently far of the LCFS, so that the wall is often left out of the computational domain by most turbulence codes.

A more precise representation of the wall, and thus a more precise account of wall plasma/interaction, is mostly the concern of transport codes as SOLPS-ITER and SOLEDGE3X-HDG, where heat deposition and particle recycling are crucial elements of particle and energy circulation. It has been shown however that this precise representation could also be performed in turbulence codes, as examplified by SOLEDGE3X. The latter code, still retaining its block-structured approach, has gained the capacity of a precise wall representation through the application of penalization methods firstly for the Bohm condition on Mach number [START_REF] Isoardi | Penalization modeling of a limiter in the tokamak edge plasma[END_REF] for isothermal plasmas, and then for Bohm boundary conditions for plasmas with non-uniform temperature, including boundary conditions for both ion and electron temperatures [START_REF] Paredes | A penalization technique to model plasma facing components in a tokamak with temperature variations[END_REF]. The penalization methods are simple in their principle: they rely on the use of a mask function, which typically takes the value 0 inside the plasma domain and 1 outside of it, and appears as a factor of extra terms added to the model equations. This allows the use of a global solver on a structured mesh, taking into account the complex geometry of the wall, as illustrated in Figure 5. A protoype equation is given in [START_REF] Stegmeir | GRILLIX: a 3D turbulence code based on the flux-coordinate independent approach[END_REF] in the form:

∂ ∂ f f = (1 -χ)F f + χ ( f p -f ) ( 22
)
where f represents a dynamical variable of the model, F f defines the time-variation of this variable in the plasma domain. Boundary conditions on f can then be enforced by suitable choices of and f p . These choices are however outside of the scope of this contribution. Suitable choices for Bohm boundary conditions have been first devised in the transport framework, but have been taken further by SOLEDGE3X to solve its anisothermal drift model in realistic geometries [START_REF] Bufferand | Progress in edge plasma turbulence modelling-hierarchy of models from 2D transport application to 3D fluid simulations in realistic tokamak geometry[END_REF]. GDB also has used the penalization method in order to enforce Bohm boundary conditions on limiters [START_REF] Zhu | GDB: A global 3D two-fluid model of plasma turbulence and transport in the tokamak edge[END_REF]. Finally, GRILLIX has used a generalization of penalization methods cited above to solve for the electric potential with Dirichlet conditions on wall boundaries [START_REF] Stegmeir | GRILLIX: a 3D turbulence code based on the flux-coordinate independent approach[END_REF] in diverted geometries.

We note here that [START_REF] Pamela | A wall-aligned grid generator for non-linear simulations of MHD instabilities in tokamak plasmas[END_REF] introduced a different block-structured meshing technique for the MHD code JOREK, using a Flux-Coordinate Dependent mesh for the core, and using wall-fitted structured mesh for the region between the LCFS and the wall. This meshing technique would be applicable by all turbulence codes using block-structured, curvilinear meshes and provides an alternative to penalization techniques for the description of the tokamak wall in the poloidal cross-section.

Discretization of fluxes associated to drifts

The crucial competence of 3D edge modelling lies in the computation of turbulence and turbulence transport. Turbulence is associated with motion of the plasma in the plane across the magnetic field, through the drift-velocities introduced in Sec. 3.1. We only discuss here the case of the particles flux but it is similar for all fluxes associated to drifts.

A first flux of interest, which is present in cold and hot plasmas such as fusion plasmas, is that carried by the E× B drift. It reads

∇ • (n v E ) = ∇ • n B×∇φ B 2 (23)
= ( B×∇φ)•∇n B 2 + n∇φ • ∇ × B B 2
The first term is the leading term. It is present in a uniform magnetic field, and has the structure of a Poisson bracket acting on φ and n, denoted by

[φ, n] = B B 2 • (∇φ × ∇n) (24)
It lends itself to a direct application of the celebrated Arakawa scheme. This scheme, initially designed for geophysical flows, has many useful properties described in its initial presentation [START_REF] Arakawa | Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow[END_REF] (conservativity, preservation of total kinetic energy and total enstrophy for the 2D Navier-Stokes equation) and has been present in fluid plasma modelling from the early days of 2D modelling until today in 3D modelling codes. Its presence has been further reinforced by studies that highlighted its superiority to conservative finite-differences using WENO reconstruction of the flux introduced in [START_REF] Liu | Weighted Essentially Non-Oscillatory schemes[END_REF] for example, in that the Arakawa scheme, formally of second order in space, is less dissipative than an equivalent discretization using 3 rd order WENO reconstruction [START_REF] Volker Naulin | Accuracy of spectral and finite difference schemes in 2d advection problems[END_REF]. It should be noted that the Arakawa scheme, because of its numerical properties, demands high quality meshes with adequate spatial resolution to avoid marked numerical instabilities, which explains that more standard schemes can be preferred. Turning back to Equation 24, its second term, associated with spatial variations of the magnetic fields is often referred to as a curvature term acting on φ, using the notation:

C(φ) = ∇ × B B 2 • ∇φ (25)
It must be noted that this curvature term is sometimes truncated:

C(φ) = ∇ × B B 2 + B × ∇B B 3 • ∇φ ≈ B × ∇B B 3 • ∇φ (26)
This truncation does omit a part of the divergence of the particle flux related to the plasma equilibrium current µ -1 0 ∇× B, which is therefore related to the poloidal magnetic field and hence of lower magnitude than the retained term. The latter is largely due to the 1/R dependence of the magnetic field, and is systematically kept as it is instrumental in the development of the interchange instability for example. In practice, the truncation is also convenient since it allows to rewrite the curvature term as a Poisson bracket and therefore to use the same discretization:

C(φ) ≈ B × ∇B B 3 • ∇φ = B B 2 • ∇B B × ∇φ = [ln B, φ] (27)
Turning now to the diamagnetic velocity, it carries a particle flux that reads:

∇ • (n v * ) = ∇p q • ∇ × B B 2 (28)
This relation illustrates that particle flux carried by the diamagnetic drift is divergence-free in a uniform magnetic field. Only the curvature term is present. Several authors have noted that it is absolutely essential to make sure that the discretization of this divergence must be constructed so as to vanish for uniform magnetic field to limit systematic errors that might occur using a naive finite-volume discretization.

Treatment of the vorticity equation and parallel Ohm's law

The vorticity equation is one of the main difficulties of fluid modelling of edge plasmas.

In electrostatic models, it requires the inversion of an elliptic, Laplacian-like operator which relates the electric potential to the vorticity Ω, Eq. 12. In a fusion plasma, resistivity is small and the term ∇ • (∇ φ/η) has to be treated implicitly to avoid too sharp a limitation of the time-step. This in turn implies that the determination of the electric potential must done through the inversion of a three-dimensional, anisotropic, Laplacian-like operator. SOLEDGE3X relies on time-scheme where perpendicular advection terms as well as the parallel plasma acoustic waves, are treated explicitly, whilst the electric potential and diffusion terms are treated implicitly. In the model used, the time-step is limited by a CFL condition fixed by parallel acoustic waves with phase velocity √ (T i + T e )/m i . If considering electron inertia, another approach can be followed as in GBS [START_REF] Paruta | Simulation of plasma turbulence in the periphery of diverted tokamak by using the GBS code[END_REF]. This approach presents the advantage that j , through nv ,i and nv ,e is then an independent variable. The smallness of resistivity is then no longer an issue, but this model comes at the cost to computational performance of having to handle electron sound waves that propagate along the magnetic field. In an explicit time-marching strategy, these waves limit the time-step through a CFL condition involving the electron sound speed c e = √ T e /m e , which is about 40 times faster than plasma sound speed if T i ≈ T e . In contrast in electromagnetic models, the parallel current is related to the potential A Eq. 14 and by a different form of parallel Ohm's law Eq. 15. Coupling Equations 12 and 15 yields a system for φ and A which describes shear-Alfvén waves. These waves are transverse and propagate along the equilibrium magnetic field with a phase velocity:

c A = B √ µ 0 nm i (29)
Again, an explicit time-marching scheme results in a CFL condition for the time-step where the velocity is c A . A thorough discussion of choices for Ohm's law focusing on the dispersion relation of linear waves is given by [START_REF] Dudson | On ohm's law in reduced plasma fluid models[END_REF]. It focuses on a minimal reduced drift model, investigating electrostatic and electromagnetic models, with or without electron inertia and solves this model for an isolated plasma filament in typical edge conditions in the cold-ion limit. It identifies that modes in the electrostatic model without electron inertia diffuse along the magnetic field with a parallel diffusivity:

D = T e ηn 0 e 1 k 2 ⊥ ρ 2 s + 1 ( 30
)
The addition of electron inertia alters the dispersion of linear modes, which become travelling waves with phase velocity equal to the electron sound speed in the large-perpendicular-wavenumber limit. It also illustrates that the inclusion of electromagnetic effects yields similar behaviour for small waves numbers with or without electron inertia, but that the inclusion of electron inertia yields pure electron acoustic waves at larger k ⊥ . Strikingly, it identifies the most complete model with electromagnetic effects and electron inertia as yielding linear modes with the longest characteristic time.

Discretization of thermal conduction

Thermal conduction is a key mechanism in the dynamics of fluctuations because of the magnitude of the thermal conductivity: although fusion plasmas exhibit large temperature variations between the core and the plasma, plasma acoustic waves behave as isothermal, rather than adiabatic, waves. This fact is particularly important in the SOL, where it sets the outflow velocity on plasma-facing components such as limiters or divertors, and the magnitude of thermal conductivity, through the boundary condition on the heat flux on the same components, determines the regime of the SOL, between conduction-limited and sheath-limited regimes [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF]. These arguments illustrate the necessity to consider thermal conduction in fluid modelling of the edge, but we illustrate here why it also presents a significant difficulty. The first issue with thermal conduction is its anisotropy in fusion plasmas: whilst particles travel freely along the magnetic field, their movement is hindered by the Lorentz force in the orthogonal plane. This fact is reflected in a significant anisotropy of the thermal conductivity. Expressions for both parallel and perpendicular conductivities are summarized in [START_REF] Huba | NRL Plasma Formulary[END_REF]. The parallel conductivity in a strongly increasing function of temperature, roughly increasing as T 5/2 α , whereas the perpendicular conductivity decreases as T -1/2 . More precisely, κ α /κ α ⊥ ∝ ω 2 cα τ 2 α , where τ α ∝ T 3/2 α is the collision time and ω cα is the cyclofrequency. The difficulty of anisotropy is augmented by the fact that the dominant conduction occurs along the magnetic field which, because of its winding, typically does not follow one the directions of the mesh. In this situation, naive discretizations are notoriously troublesome. This difficulty is now regularly solved in codes using structured meshes by implementing the scheme proposed in [START_REF] Günter | Modelling of heat transport in magnetised plasmas using nonaligned coordinates[END_REF]: the anisotropic Laplacian is treated using a conservative finite-difference method where fluxes are evaluated on a staggered grid. This methods leads typically to a 9-point stencil, classical for anisotropic diffusion, but with the advantages of being conservative, as well as reproducing at the discrete level the property that the anisotropic Laplacian is semi-definite negative. It is extended in a straightforward manner to structured curvilinear meshes with non-uniform conductivity. It has been recognized that the latter scheme has the property it proposes a discretization of the anisotropic Laplacian which can written as the compose of discrete divergence and discrete gradient operators which are mutually anti-adjoint, and follow the Support-Operator Method presented in [START_REF] Shashkov | Support-operator finite-difference algorithms for general elliptic problems[END_REF], also known as Mimetic Finite-Differences. This ensures the semi-negativity of the discrete anisotropic Laplacian, which in turn guarantees that thermal conduction yields decay of the temperature field in the L 2 norm when using an implicit scheme or an explicit scheme with suitable bound on the time-step. It has been widely used in discussed in the community as providing a robust scheme that deals particularly well with anisotropy of general structured grids. It however does not reproduce the monotonicity of thermal conduction at a discrete level, as only a non-linear discretization can achieve this property [START_REF] Sharma | Preserving monotonicity in anisotropic diffusion[END_REF].

Nevertheless, as noted earlier, several codes have made the choice of using a flux-aligned mesh, with the consequence that parallel conduction is limited to a 2D-lattice. The advent of field aligned discretizations has drawn attention to develop dedicated methods, starting with [START_REF] Bram Van Es | Finite-difference schemes for anisotropic diffusion[END_REF], followed by [START_REF] Stegmeir | The field line map approach for simulations of magnetically confined plasmas[END_REF] to introduce SOM in this context, and [START_REF] Soler | A new conservative finite-difference scheme for anisotropic elliptic problems in bounded domain[END_REF] to further ensure conservativity of the discretizations.

The final difficulty presented by thermal conduction is the magnitude of parallel conductivity in itself: combined with knowledge of the time-scales of interest, close to the acoustic transit time based on the mesh size, treatment of thermal conduction with implicit time schemes yields ill-conditioned linear systems which are problematic in terms both of accuracy and convergence with iterative methods. The accuracy issue has received little attention in the community, although theoretical investigations in the framework of Asymptotically-Preserving schemes have identified that condition number, which would scale with κ in the large conductivity limit in a standard discretizations, could be made both small and independent of conductivity in the same large conductivity limit [START_REF] Degond | An Asymptotic Preserving scheme for strongly anisotropic elliptic problems[END_REF]. Although the reduction in conditioning, and consequently the improvement in accuracy, is significant, this contribution has found little resonance in the plasma community. This is probably due to the emphasis on solution techniques for the linear problem posed by parallel conduction, and the latter theoretical investigations have left out the aspect of iterative solutions of linear systems, which is crucial for the development of edge turbulence simulations. The second point worth mentioning is the advent of efficient black-box libraries, such as Petsc [START_REF] Balay | PETSc/TAO users manual[END_REF], which have enabled modellers to test various iterative solution strategies without the effort of developing the algorithms and data structures necessary for these strategies.

Time-marching schemes

The variety in models, meshing and spatial discretization strategies is also found in time-marching strategies. Two extreme strategies can be found, in BOUT++ which uses a full implicit, time-adaptive strategy, and GBS which uses a full explicit time-marching schemes. The positioning of BOUT++ is however more unique. Other codes use semiimplicit schemes, with a view to use time-steps as close to the smallest time-scale in the models as possible. We have explained previously that models vary between codes, with two major differences in the account of magnetic fluctuations firstly, and of electron inertia secondly, but we recall once that all codes are interested in the first place in turbulence occurring on ion time and length-scales, so that both electron acoustic waves and Alfvén waves for electromagnetic codes are arguably considered out of the scope of these codes. This view is implicit in the frequent choice in the characteristic time made τ ref = R 0 /c s , which can be seen as a toroidal acoustic transit time. As in fluid modelling, explicit time-stepping has the advantage of requiring fewer operations per time-step, and in a parallel environment, of requiring a benign amount of inter-processor communication, which is favorable for the parallel scaling of the computation. It has the downside of sharp limits on the time-step imposed by stability. In particular, this can become particularly limiting when considering realistic values of certain model parameters. On the other hand, implicit schemes typically require solution of global linear systems, which induce a larger number of operations per time-step and are generally communication-intensive. As noted, the vorticity equation coupled with Ohm's law represents one of the major difficulties in fluid modelling of the edge: advection by drifts as well as parallel acoustic waves are local in nature whilst the relation between vorticity and electric potential for all models, and the relation between current and parallel vector potential for electromagnetic models involve Laplacian-like, and hence elliptic, operators that must be inverted globally. BOUT++ clearly lies at one end of the spectrum with its fully implicit, time-adaptive timeadvancement which suffers no CFL time-step restriction at the expense at a global solve of a non-linear problem at each time-step. GBS lies at the other end of that spectrum, with a fully explicit fourth-order Runge-Kutta scheme where the time-step is restricted by a CFL condition on the numerical electron parallel transit time, but the relations giving electric potential and vector potential are uncoupled, and both can be obtained by inversion of the perpendicular Laplacian. GBS further uses an approximation of the perpendicular Laplacian so as to make it local in toroidal angle, which favours parallelism for the determination of φ and A . GRILLIX uses third-order, implicit-explicit splitting method initially devised for the Navier-Stokes equation [START_REF] Em Karniadakis | High-order splitting methods for the incompressible Navier-Stokes equations[END_REF], in which the parallel current, advection and acoustic terms are treated explicitly in its EM version, and only thermal conduction is treated implicity. Thus, GRILLIX carries computational advantages similar to those of GBS, while further exploiting its aligned discretization, which allows coarse toroidal resolution and renders the CFL condition on parallel EM dynamics less stringent. Finally, SOLEDGE3X solves its ES model by treating parallel current and thermal conduction implicitly and all other terms explicitly. This allows the use of large-time steps, but at the expense of inverting a global, strongly anisotropic 3D Laplacian stemming from the time discretization of the vorticity equation:

∂ ∂t ∇ • m i n∇ ⊥ φ ZB 2 + ∇ 2 φ η = o.t. (31)
Typically, this problem is solved using Krylov methods preconditioned by Algebraic Multigrid Methods (see [START_REF] Bufferand | Progress in edge plasma turbulence modelling-hierarchy of models from 2D transport application to 3D fluid simulations in realistic tokamak geometry[END_REF] and references therein). These strategies reflect different choices in the necessary trade-off between computational time per time-step and size of the time-step when opting for explicit or implicit treatment. They are enlightened by the study given in [START_REF] Dudson | On ohm's law in reduced plasma fluid models[END_REF]: this study investigated the computational cost of modelling an isolated plasma filament with electrostatic/electromagnetic model, with or without electron inertia. We recall that this study was conducted using the code STORM [START_REF] Easy | Three dimensional simulations of plasma filaments in the scrape off layer: A comparison with models of reduced dimensionality[END_REF], using the CVODE implicit time solver from the SUNDIALS suite [START_REF] Hindmarsh | SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers[END_REF], which adapts the time-step so as to satisfy tolerance criteria, as well as reduce the number of inner iterations per time-step. The time-step used by such a solver then provides information on the fastest time-scale of interest in the simulation. Interestingly, the time-step used is the largest for EM model with electron inertia, whilst the smallest time-step is used for the electrostatic model without electron inertia. The number of inner iterations is also the smallest for the EM model without electron inertia, so that the wall clock time is around 13 times shorter than for the ES model without electron inertia.

Verification and validation

The numerical modelling of tokamak edge plasma shows specific complexities compared with other communities, such as fluid dynamics as mentioned in Ref. [START_REF] Terry | Validation in fusion research: Towards guidelines and best practices[END_REF]. Indeed, the fluid model is not unique and the multiphysical nature of the problem requires the integration of different models with different scales or representations, for example for neutrals or impurities. Moreover, in experiments, the tokamak environment makes that diagnostics are scarce with limited capability (many crucial quantities cannot be measured), and require of applying significant a-priori modeling to interpret their measurements.

Verification of codes in the literature therefore remains limited compared with computational fluid dynamics [101]. For a long time, codes were verified using test cases built on simplified physical problems for which analytical results were available, such as an instability linear growth rate, a wave phase velocity or a blob propagation speed (see for example in [START_REF] Dudson | BOUT++: A framework for parallel plasma fluid simulations[END_REF][START_REF] Easy | Three dimensional simulations of plasma filaments in the scrape off layer: A comparison with models of reduced dimensionality[END_REF]). Such a procedure does not stress systematically all terms in the equations, and thus only allows a partial verification of the code. A more rigorous method of verification is provided by the popular method of manufactured solutions. To the best of our knowledge, it was only introduced about 10 years ago to the fusion community [START_REF] Riva | Verification methodology for plasma simulations and application to a scrape-off layer turbulence code[END_REF][START_REF] Ricci | Approaching the investigation of plasma turbulence through a rigorous verification and validation procedure: A practical example)[END_REF][START_REF] Dudson | Verification of BOUT++ by the method of manufactured solutions[END_REF][START_REF] Tamain | The TOKAM3X code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF]. Using a stable target analytical solution to approach (generally a combination of trigonometric functions to make the computation of the source term trivial), this method guarantees quantitatively the correct implementation of the model equations and of the numerical scheme. In recent years, it has been used more systematically to verify codes [START_REF] Tamain | The TOKAM3X code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF][START_REF] Stegmeir | GRILLIX: a 3D turbulence code based on the flux-coordinate independent approach[END_REF][START_REF] Riva | Shaping effects on scrape-off layer plasma turbulence: A rigorous validation of three-dimensional simulations against tcv measurements[END_REF][START_REF] Giacomin | The GBS code for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary[END_REF]. More recently, a method based on the output data analysis (PoPe, Projection on Proper elements) has been proposed to go further in the verification process of simulations in any regime including turbulence [START_REF] Cartier-Michaud | A posteriori error estimate in fluid simulations of turbulent edge plasmas for magnetic fusion in tokamak using the data mining iPoPe method[END_REF]. The key idea is to quantify the numerical error by performing a data-driven identification of the mathematical model from the simulation outputs. Based on a statistical post-processing of the output database, the method provides a measure of the error by estimating the distance between the (numerical) effective and (analytical) theoretical weights of each operator implemented in the mathematical model. The method therefore not only assesses that the equations actually solved numerically are the equations that are claimed to be solved, as the MMS, but it also allows to address the question of accuracy of the numerical resolution [START_REF] Cartier-Michaud | Verification and accuracy check of simulations with PoPe and iPoPe[END_REF]. Furthermore, it presents the advantage over MMS that the verification can be performed during a simulation.

Validation of codes consists of comparing simulations and experiments on a number of common physically relevant quantities which define a set of independent observable like the average plasma density (n) or temperature (T i,e ), the ion saturation current density or the electric potential Φ. The basics of the validation methodology were initially defined in Refs. [START_REF] Terry | Validation in fusion research: Towards guidelines and best practices[END_REF][START_REF] Greenwald | Verification and validation for magnetic fusion)[END_REF][START_REF] Ricci | Approaching the investigation of plasma turbulence through a rigorous verification and validation procedure: A practical example)[END_REF].

It remains a particularly challenging task in fusion because in addition to the classical uncertainties present in numerical simulations, many specific uncertainties with experimental measurements must be taken into account as in the a priori modelling used to interpret the measures or the difficulty of setting accurately the control parameters in the tokamak.

Most of the numerical works published in the literature provide comparisons with experiments in more or less complex configurations. A synthetic vision on the state-of-the-art in the codes validation is certainly provided in the series of benchmarks published these recent years [START_REF] Militello | Multi-code analysis of scrape-off layer filament dynamics in MAST[END_REF][START_REF] Riva | Blob dynamics in the torpex experiment: a multi-code validation[END_REF][START_REF] B D Dudson | Edge turbulence in ISTTOK: a multi-code fluid validation[END_REF][START_REF] Galassi | Validation of edge turbulence codes in a magnetic X-point scenario in TORPEX[END_REF][START_REF] Oliveira | Validation of edge turbulence codes against the TCV-X21 diverted L-mode reference case[END_REF][START_REF] Moscheni | Cross-code comparison of the edge codes SOLPS-ITER, SOLEDGE2D and UEDGE in modelling a low-power scenario in the DTT[END_REF]. In these works, multi-code validation exercises are performed with respect to experimental measurements in various tokamaks in operation. The results show that all turbulence codes are able to reproduce qualitatively some key features of the plasma discharge in realistic configurations. A very nice example of qualitative agreement is given on Fig. 6 showing the structure of a filament characteristic of turbulence post-treated using a synthetic diagnostic from a STORM simulation and the corresponding experimental visualisation on MAST.

However, quantitative comparisons in medium-sized tokamaks and relevant parameter ranges remain challenging, mostly due to the lack of physics in the resolved models and the large uncertainties on experimental measurements. As an example, Fig. 7 shows recent comparisons between numerical 1D radial profiles predicted by state-of-the-art turbulence codes and TCV experimental data at the high-field-side divertor target [START_REF] Oliveira | Validation of edge turbulence codes against the TCV-X21 diverted L-mode reference case[END_REF]. If some general trends can be recovered on the numerical profiles, a full quantitative agreement between simulations and experiments is far from being achieved. The numerical results do not match within experimental uncertainty, and the profiles vary significantly between the simulations, which indicates the high sensitivity of the results to the implemented models.

More satisfactory quantitative results are, however, obtained in more academic configurations with a smaller and simpler geometry and for more tractable parameter ranges (see for example the comparisons in TORPEX in Ref. [START_REF] Riva | Blob dynamics in the torpex experiment: a multi-code validation[END_REF]).

Some examples of recent achievements

The increasing efficiency of both numerical algorithms and computational resources, in particular through massive parallelization of the codes, has allowed in recent years to continuously increase the number of degrees of freedom in the simulations. As a result, recent achievements in the literature show that medium-sized tokamak-scale turbulence simulations (plasma volumes around 15m 3 ) can now be achieved. However, the simulations remain computationally very demanding, and most of them do not integrate yet all the theoretically required physics in the model, particularly that linked to wall-plasma interactions. As a consequence, edge turbulence simulations on the ITER scale with a complete physics are still inaccessible, even if a pioneering ITER simulation has been very recently performed with GRILLIX until the formation of an initial instability at 27 microseconds [START_REF] Stegmeir | Analysis of locally-aligned and non-aligned discretisation schemes for reactor-scale tokamak edge turbulence simulations[END_REF]. Some examples of recent achievements have been picked up in the literature to illustrate the current capacity of the state-of-the-art codes. Corresponding time steps and meshes are provided in Table 3. As expected, time-step values and spatial resolutions show the impact of the numerical options (explicit/implicit, aligned/disaligned) made in the different codes to discretize the equations (See Table 2), as discussed in Sec.4. Note, however, that not all simulations considering the same range of control parameter values can damp the differences between the time steps and grids used by the different numerical schemes.

Codes

Tokamak Grids Time-step Simu time (scale, 2π/p) (ns) (ms) BOUT [START_REF] Zhu | Electromagnetic turbulence simulation of tokamak edge plasma dynamics and divertor heat load during thermal quench[END_REF] DIIID 260 × 64 × 64 -0.345 (1, 2π/5) GBS [START_REF] Giacomin | The GBS code for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary[END_REF][START_REF] Mancini | Self-consistent multi-component simulation of plasma turbulence and neutrals in detached conditions[END_REF] TCV 150 × 300 × 64 0.2 2 (1/2, 2π/1) GRILLIX [START_REF] Zholobenko | Filamentary transport in global edge-SOL simulations of ASDEX Upgrade[END_REF] ASDEX-U 440.10 3 × 16 1.2 5 (1, 2π/1) SOLEDGE3X [START_REF] Bufferand | team. Global 3D full-scale turbulence simulations of TCV-X21 experiments with SOLEDGE3X[END_REF] TCV 220 × 1300 × 64 2 0.8 (1, 2π/4) Table 3: Examples of numerical parameters used in the literature by the different codes to obtain the recent results presented below. Values depend on the configuration to simulate (geometry and parameter ranges), as well as on the numerical schemes.

As an example, in Giacomin et al. [START_REF] Giacomin | The GBS code for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary[END_REF] the grid sensitivity analysis shows that a mesh N R × N Z × N φ = 150×300×64 provides satisfactory results in a GBS simulation of a lower single-null discharge in a computational domain of halfsize of TCV and for the range of control parameters used. Regarding time scales now, the turbulence saturation time can be long with respect to very small time steps required by the numerical stability constraint imposed by the CFL condition. Typically, the turbulence saturation times are of the order of the milliseconds while the time-steps are of the order of the nanoseconds (see examples in [START_REF] Giacomin | The GBS code for the self-consistent simulation of plasma turbulence and kinetic neutral dynamics in the tokamak boundary[END_REF][START_REF] Zholobenko | Electric field and turbulence in global Braginskii simulations across the ASDEX Upgrade edge and scrape-off layer[END_REF]), with up to one order of magnitude between the codes. A full scale simulation of an attached L-mode discharge in ASDEX-U (Table 1) is shown on Fig. 8. In this simulation, the plasma is pure Deuterium and neutral gas is modelled using a simple diffusion model, which interacts with the plasma through charge exchange and ionisation reactions. The 3D plot of the instantaneous density fluctuations provided by GRILLIX [START_REF] Zholobenko | Filamentary transport in global edge-SOL simulations of ASDEX Upgrade[END_REF] shows very fine filament structures aligned along the magnetic field line in the parallel direction. They expand radially around the separatrix with a strong ballooning on the high field side as expected from interchange driven turbulence. The plot also shows that density structures develop along the divertor leg. With a resolution in the poloidal plane of 2.5ρ and 16 toroidal planes (the field-aligned approach allows a rather coarse resolution in the toroidal direction, Sec. 4.1.1 and 4.1.2) the corresponding grid is of about 7 million points. As mentioned above, due to the very small time step imposed by stability constraint (dt = 1.2ns), the simulation requires roughly 0.4 MCPUh and 2 months on 384 processor cores (8 nodes) of the Marconi-A3 SKL to reach the turbulence steady state [START_REF] Zholobenko | Electric field and turbulence in global Braginskii simulations across the ASDEX Upgrade edge and scrape-off layer[END_REF]. Similarly, a TCV simulation at full scale has been recently performed using SOLEDGE3X for a pure Deuterium plasma and neutral fluid [START_REF] Bufferand | team. Global 3D full-scale turbulence simulations of TCV-X21 experiments with SOLEDGE3X[END_REF]. The resolution of 2ρ in the poloidal plane is equivalent to the former simulation in ASDEX-U, and corresponds here to a mesh of 17 millions of points N Ψ × N θ × N φ = 220 × 1300 × 64. Fig. 9 shows density and electron temperature fluctuations in the poloidal plane with a zoom around the divertor leg. The simulation recovers in particular experimental turbulence features in the divertor such as Divertor Localized Filaments with very fine structures, quite sheared, and expanding far along the divertor leg. Such simulation ran 12 days on 3072 nodes that corresponds to a total plasma time simulated of about 0.8ms. Recent advances also allow to handle additional physics in the simulation such as external 3D magnetic perturbations [START_REF] Boinnard | Plasma turbulence simulations in a diverted tokamak with applied resonant magnetic perturbations[END_REF], electromagnetic effects [START_REF] Zhu | Electromagnetic turbulence simulation of tokamak edge plasma dynamics and divertor heat load during thermal quench[END_REF] or more sophisticated interactions of the plasma with neutral gas, impurities and machine walls [START_REF] Bufferand | Implementation of multi-component Zhdanov closure in SOLEDGE3X[END_REF][START_REF] Mancini | Self-consistent multi-component simulation of plasma turbulence and neutrals in detached conditions[END_REF]. Fig. 10 illustrates here BOUT++ electromagnetic simulations in a lower single-null plasma equilibrium with generic DIII-D tokamak H-mode parameters. Magnetic field fluctuations are shown in the poloidal plane together with radial outer-divertor heat-load profiles at different instants. The results show in particular that even for electrostatic instabilitydominated turbulence, turbulence enhances magnetic fluctuations at such a level that they completely break magnetic flux surfaces such that stochastic field lines are directly connecting pedestal top plasma to the divertor target plates or first wall, further contributing to the divertor heat-flux width broadening [START_REF] Zhu | Electromagnetic turbulence simulation of tokamak edge plasma dynamics and divertor heat load during thermal quench[END_REF]. Such very pioneering results obtained on a very short physical time should be however confirmed on a longer time in further works.

Finally, we highlight here recent simulations no longer limited to a pure Deuterium plasma but based on a multispecies fluid model using a Zhdanov closure (see Sec.3.5). Such modelling allows higher fidelity in terms of plasma Extracted from [START_REF] Zhu | Electromagnetic turbulence simulation of tokamak edge plasma dynamics and divertor heat load during thermal quench[END_REF].

interactions between neutrals and impurities through ionization, charge-exchange, recombination and molecular dissociation processes. 3D simulations of a D+T+He plasma with impurities like C or W in WEST [START_REF] Bufferand | Implementation of multi-component Zhdanov closure in SOLEDGE3X[END_REF], or a D + and D + 2 plasma interacting with two neutrals species (D and D 2 ) in TCV [START_REF] Mancini | Self-consistent multi-component simulation of plasma turbulence and neutrals in detached conditions[END_REF] have been very recently performed. In the latter in TCV, detachment conditions have been achieved for the first time in turbulence simulations at the inner target. With divertor conditions controlled by D 2 puffing, the increase in fuelling leads to the decrease of the particle and heat fluxes to the divertor target, and to the detachment of the inner target as shown on Fig. 11. It is a result of primary importance since detachment is expected to be the optimal regime for divertor operation in future reactor as ITER to harness heat exhaust.

Conclusion and Discussion

This article aims at reviewing recent advances in the literature related to fluid numerical models to perform turbulence simulation of tokamak edge plasmas. Despite the constant growth of computational resources along with significant improvements in numerical methods, fluid modeling remains standard at the edge in the international fusion community. This is in particular due to the geometrical complexity of the wall and of the magnetic equilibrium with singularities (Sec. 2), as well as to the complex plasma wall interactions which require to consider multi-species plasmas and eventually multi-physics phenomena. Moreover, with the arrival in a close future of thermonuclear reactors as ITER, there is a real need for relative fast computations (what first principles (gyro)kinetic models cannot yet offer) allowing parametric studies and the design of optimal scenarios to be used as guidelines for operation performance and safety of the machine.

Computing the moments of the Boltzmann equation leads under various assumptions to a set of complex nonlinear conservation equations for fluid quantities (Sec. 3.2) together with dedicated boundary conditions modelling at fluid level plasma wall interaction (Sec. 3.4). If most of the models were mainly focused on electrostatic turbulence in a simple quasi-neutral hydrogenic plasma, in recent years models have also considered multi-species plasma (Sec. 3.5), and incorporated electromagnetic effects allowing fluctuations of the magnetic field (Sec. 3.2). With the strong magnetization inside the tokamak, the parallel flow has a fast dynamics (characteristic time of the order of the parallel length over the sound speed c s = √ (T i + T e )/m i ), i.e. of the order of ∼ 10 -4 s in current machines) and behaves as that of a compressible gas, while the perpendicular flow mostly driven by turbulence has a much slower dynamics (characteristic time of the order of the energy confinement time ∼ 1s in current machines) and behaves as that of an incompressible fluid. Correspondingly, the turbulence is quasi-2D with very short-wavelength fluctuations in the poloidal plane (of the order of the Larmor radius) coupled to long wavelength in the parallel direction (of the order to the parallel connection length).

All these features make these equations specially challenging to resolve and stressful for the numerical methods and algorithms. Thus, works in the literature show innovative numerical methods in particular for the spatial discretizations as described in Sec.4. Probably one of the most original aspect compared to CFD lies in the meshing and discretization strategies developed to take advantage of the strong anisotropy of the flow as well as the complex geometric characteristics to be dealt with. As shown in the paper, no consensus on the optimal combination between meshing and discretization, if it exists, has been found. Indeed, if non-aligned approaches make easier handling boundary conditions, particularly when the boundary conditions mix parallel and perpendicular directions or if the wall geometry becomes very complex, aligned approaches take advantage of the turbulence anisotropy and allow a strong coarsening of the mesh in the parallel direction which is very computationally efficient. Also, when treated implicitly to avoid too sharp a limitation of the time-step, another issue set in the inversion of very ill-conditioned Laplacian-like operators, whether for the vorticity equation or for the thermal conductivity, which can be problematic in terms both of accuracy and convergence with iterative methods.

A set of reference fluid codes have been developed by the international community embedding all recent advances in numerical methods described in this review. Nowadays, all these codes have been carefully verified, and their validation with respect to experimental measurements is the focus of the teams' ongoing efforts and is making significant progress (Sec.5). In the current state-of-the-art, mean-field transport codes, 2D in their vast majority, remain the only tools capable of addressing the plasma exhaust challenge in large scale devices. Their geometrical flexibility and capability to perform simulations in a wide range of density regimes for current machines as well as future devices are obtained at the cost of compromises on the description of plasma anomalous transport, thus hindering their predictive capabilities. 3D fluid turbulence codes, on which this review was focused, pave the way to go beyond this limitation by offering a direct modelling of turbulent transport. However, their integration in a complete edge plasma physical environment (selfconsistent inclusion of neutrals and impurities, realistic magnetic and wall geometry. . . ) is at its premises, and raises physical and technical difficulties even if very encouraging medium-sized tokamak-scale 3D turbulence simulations are now obtained (Sec. 6). While achieving spatial resolution at turbulence scales requires significant resources, the range of time scales to be handled by these codes as well as the required physical time per iteration seems even more challenging. Indeed, as shown in the paper, the characteristic times evolve from the order of few nanoseconds, as imposed on the time-step by numerical stability constraints (fixed by parallel acoustic waves phase velocity) to the order of few milliseconds to reach turbulence saturation, and even longer to satisfy particles and energy balance. This time can be even larger if considering neutral physics. Additionally, and probably differently from most CFD problems, the ill-conditioned Laplacian-like operators to invert, the multiple plasma variables to take into account as well as the possible code coupling required to model plasma wall interactions lead to a significant computational time per iteration, despite the relatively high level of parallelization of the solvers running on thousand of cores in Petascales supercomputers.

Thus, despite this HPC environment, computing time seems to be for the moment a barrier for the application of these 3D turbulence codes for large-sized tokamaks as JET, JT60SA or ITER.

Perspectives in developments

This review shows all the recent progress made by the fusion community in the development of innovative numerical methods and algorithms in HPC framework to pave the route towards predictive 3D turbulence simulations of reactor relevant devices such as ITER. However, this work also highlights the gap that remains to be filled to achieve this. Thus, the perspectives for development are numerous, both in physics modelling and numerics, but always guided by the search of the best compromise between fidelity to the targeted physics objectives and computational cost. This is the purpose of fluid turbulence models, at the midpoint between more expensive, higher-fidelity, high-dimensional models ((gyro)-kinetic models) and cheaper, lower-fidelity, low-dimensional ones (1D or 0D fluid models). Finding this compromise remains an open issue since unlike most CFD configurations, it is not straightforward to estimate a priori what is the minimum set of physical ingredients and numerical precision required to achieve the desired objectives, knowing that up to now full-scale simulations with all physics and resolution of all scales are not achievable.

Below, some of the main development avenues that could be declined along the following lines:

Physics modelling:

• As mentioned in Sec.3.5, current drift-fluid models miss reactor relevant kinetic effects, such as small amplitude fluctuations at the gyro-radius scale and low collisionality effects that can be very relevant at the edge. Fullf gyro-fluid models as implemented in FELTOR appear as the most balanced compromise between embedding kinetic effects and computational cost. They could become the future standard for fluid turbulence codes. However, fullf gyro-fluid models currently have to face several theoretical issues which ought to be solved before practical implementation is attempted in an integrated production code (i.e., including neutrals, impurities. . . ).

• As mentioned in Sec. 3.5 a fusion plasma is inherently composed of multiple ion species and there are many experimental evidence that impurities have a significant effect on plasma dynamics (see an example in [START_REF] Maggi | Isotope effects on l-h threshold and confinement in tokamak plasmas[END_REF]. Thus, the development and tests of stable and efficient numerical implementations for multi-ion species driftfluid models (Zhdanov closure) as already done in SOLEDGE3X and GBS is required and should progressively replace the classical Braginskii closure.

• The key role of neutral particles and their interaction with the plasma (Sec. 3.3) calls for a more systematic self-consistent treatment of turbulence and neutral particles physics. It is clearly a key challenge to target predictive simulations both for physical and numerical modelling. Indeed, a hierarchy of neutrals models have to be considered ranging from advanced fluid models to full kinetic ones in order to find the best compromise between fidelity and computing cost. These models could be directly embedded in the codes (for examples en [START_REF] Wersal | A first-principles self-consistent model of plasma turbulence and kinetic neutral dynamics in the tokamak scrape-off layer[END_REF]) or via codes coupling to EIRENE [START_REF] Fan | Self-consistent coupling of the three-dimensional fluid turbulence code TOKAM3X and the kinetic neutrals code EIRENE[END_REF] as it is currently done in transport codes (see below).

• While full turbulence models provide the most complete transport description, in the near future their computational cost is expected to remain prohibitive for divertor design studies, which require extensive parameter scans.

Transport codes will thus remain very probably the work-horses for these studies. However, improving their predictive capability requires the development of self-consistent perpendicular transport models in the spirit of RANS in CFD to replace the constant ad hoc coefficients currently used. Following few pioneering works in the literature [START_REF] Baschetti | Self-consistent cross-field transport model for core and edge plasma transport[END_REF][START_REF] Coosemans | A self-consistent mean-field model for turbulent particle and heat transport in 2D interchange-dominated electrostatic E×B turbulence in a sheath-limited scrape-off layer[END_REF] several important theoretical questions remain to be addressed to improve the accuracy and predictive capability of these reduced models using both simulations and experimental data. For the simulation, dedicated turbulence simulations will have to be performed to provide relevant data bases. Such improvements could be based on artificial intelligence or on data assimilation as recently proposed in Ref. [START_REF] Auroux | Asymptotic behavior, non-local dynamics, and data assimilation tailoring of the reduced kmodel to address turbulent transport of fusion plasmas[END_REF] in a variational data assimilation theoretical framework.

• The current boundary conditions (BC) at the plasma-wall interface presented in Sec.3.4 and based on a single ion 1D collisionless plasma sheath model with the Maxwellian velocity distributions of plasma particles are known to be not rigorously valid for realistic plasma edge conditions (see [START_REF] Tskhakaya | One-dimensional plasma sheath model in front of the divertor plates[END_REF] and references there). Thus, further developments of a collisional sheath model including drift effects and multi-ion species are required.

Numerics:

• In (semi)-implicit time discretizations, a bottleneck in code performance is linked to the inversion of strongly anisotropic differential operators which, associated to specific boundary conditions at the wall, can be potentially very ill-conditioned. This is particularly problematic for iterative algorithms and requires the development of linear algebra solvers that are ever more efficient in terms of precision and numerical cost, but also in terms of memory requirements on HPC architectures based on GPU.

• A limiting factor for the simulations remains their computational cost. MPI communication and OpenMP performance are currently being optimized, which should allow for faster medium-sized tokamak simulations, also including magnetic flutter. For larger devices like ITER and DEMO, a larger speed-up or scalability is required, by means of GPU acceleration or 3D MPI domain decomposition. With the expected delivery of exascale machines in the next few years, progressively opens the path to 3D turbulence simulations for reactors provided strong investment is made in the numerical optimization of the codes.

• Due in large part to the increasing complexity of physical models, and the inherent introduction of random model data (particularly with each atomic physics phenomenon added to describe plasma wall interactions), it then becomes crucial and urgent to study uncertainty quantification and their propagation in these codes.

• Experimental validation of models and numerical codes is a long term process, requiring the execution and analysis of dedicated experiments in various plasma regimes. On this way a more systematic quantification of uncertainties in particular on the modelling side has to be carried out. Modern tools as uncertainties propagation and parameters sensitivity analysis as used in CFD will have to be developed. In addition, the development of robust and reliable synthetic diagnostics become a priority for codes validation as proposed recently in Ref. [START_REF] Kudashev | Development of a set of synthetic diagnostics for the west tokamak to confront 2d transport simulations and experimental data[END_REF]. Indeed, whereas models describe the physics in terms of quantities in SI units, e.g. electron density, electron and ion temperatures and plasma potential, raw data from experimental observations are in currents and voltage. These synthetic diagnostics use analytical models to convert data from simulations into signals emulating a real experimental diagnostic, thus serving as important tools bridging analytical models and experiments data. Synthetic Langmuir probes or obtaining fluctuating data in the MHz range, as well as different spectroscopy diagnostics to measure turbulence in the +100 kHz range in the SOL have to be integrated into the codes.

• Code coupling may become mandatory to improve the predictive capabilities of turbulence codes. This can occur for example when accurately taking into account plasma wall interactions (neutral and impurities with codes as EIRENE [START_REF] Reiter | The EIRENE and B2-EIRENE codes[END_REF], DEGAS2 [START_REF] Russell | A reduced model of neutral-plasma interactions in the edge and scrape-off-layer: Verification comparisons with kinetic Monte Carlo simulations[END_REF], ERO2 [START_REF] Romazanov | Beryllium global erosion and deposition at jet-ilw simulated with ero2.0[END_REF], ...), boundary conditions of the sheath (Particle in Cell Monte Carlo collision codes as BIT3 [START_REF] Tskhakaya | One-dimensional plasma sheath model in front of the divertor plates[END_REF]), when extending the simulations to the core in order to avoid artificial boundary conditions at the core/edge interface (1D code as Qualikiz [START_REF] Van De Plassche | Fast modeling of turbulent transport in fusion plasmas using neural networks[END_REF]), or when taking into account consistency between the magnetic field equilibrium and the plasma (with Grad-Shafranov codes as NICE [START_REF] Faugeras | An overview of the numerical methods for tokamak plasma equilibrium computation implemented in the nice code[END_REF]). If such coupling is rather standard in transport regimes it is still in infancy for turbulence codes. Coupling strategies will have to be developed to allow efficient data exchanges between codes including interpolation functionalities when necessary. The specificity of dealing with a time-dependent system, including fluctuations, raises questions which have not been addressed in the frame of mean-field modelling as the order of the time-discretization scheme, performance optimization, or the treatment of the time dynamics of neutrals. The global impact of the coupling on codes performances, in particular the speed up of massively parallel codes, must be also addressed to make this approach valuable.

• Simulations of plasma edge turbulence with a 3D magnetic field are still in their early stages (see for example in Refs. [START_REF] Shanahan | Fluid simulations of plasma filaments in stellarator geometries with BST-ING[END_REF][START_REF] Coelho | Global fluid simulation of plasma turbulence in a stellarator with an island divertor[END_REF]), and will need to be developed to investigate configurations relevant to fusion operations. Whether it is only small amplitude disturbances which are 3D as RMPs or ripple, or the magnetic equilibrium itself as in the stellerators, numerical developments are necessary depending on the numerical options discussed in this review and embedded in the codes.

Figure 1 :

 1 Figure1: Sketch of the diverted plasma geometry in a tokamak. (left) Toroidal topology of a tokamak plasma, which can be categorised into the core, the edge, the scrape-off layer (SOL) and the private flux region (PFR). The separatrix separates the closed flux-surface edge region from the SOL and PFR regions, where plasma streaming along open magnetic field lines is bounded by the divertor plates. The outer SOL is also limited by the first wall. In the core and edge the magnetic field lines are either closed, fill a poloidal flux-surface or volume of space (extracted from[START_REF] Held | Full-f gyro-fluid modelling of the tokamak edge and scrape-off layer[END_REF]). (right) Coordinates system and geometrical parameters R 0 be the major radius, a the minor radius on the magnetic axis and (R, Z, ϕ) a fixed cylindrical coordinates system. (ψ, θ , ϕ) curvilinear system of coordinates (extracted and modified from[START_REF] Lütjens | The CHEASE code for toroidal MHD equilibria[END_REF]).

Figure 2 :

 2 Figure 2: Illustration of a possible multi-block structured domain decomposition used in FCDM for the modelling of diverted plasmas. Extracted from [27].

Figure 3 :

 3 Figure 3: Examples of numerical grids within poloidal plane: (FCDM) (left) (extracted from Ref.[START_REF] Wiesenberger | Streamline integration as a method for structured grid generation in X-point geometry[END_REF]). The X-point has a different connectivity with respect to the other cells and a coarser resolution around it. (FCIM) (right) (extracted from Ref.[START_REF] Stegmeir | GRILLIX: A 3D turbulence code for magnetic fusion devices based on a field line map[END_REF]). The Cartesian mesh allows a homogeneous resolution and the X-point is no longer a singularity for the mesh.

Figure 4 :

 4 Figure 4: General overview of FCIM concept in tokamak geometry. (Left) Grid, i.e. set of Cartesian poloidal planes, is limited by outer (Ψ max ) and optionally inner (Ψ min ) flux surface and divertor plates. Parallel operators are discretised via field line map, i.e. field line tracing and interpolation.(Right) Scheme for toroidally staggered FCIM approach. Grid is indicated with crosses and dual grid with circles. Extracted from[START_REF] Stegmeir | GRILLIX: a 3D turbulence code based on the flux-coordinate independent approach[END_REF].

Figure 5 :

 5 Figure 5: Illustration of the working principle of the penalization method. The plasma volume is in blue, whilst the penalized volume is in blue. (a) Entire poloidal cross-secion. (b) Zoom on the left side of the bottom limiter. The wall boundary in black denotes the separation between the two (Extracted from[84]).

Figure 6 :

 6 Figure 6: Comparison between filament from visual imaging (left) and a corresponding STORM simulation plotted using the synthetic D α emission diagnostic (right). Extracted from [109].

Figure 7 :

 7 Figure 7: Comparison of averaged profiles at the high-field-side divertor target as a function of the radial distance from the separatrix (in cm). (A) Profiles of the mean plasma density (n e ), (C) electron temperature (T e ), (E) plasma potential (Φ) and (G) the parallel current density ( j ) in TCV.The experimental data from the high-field-side Langmuir probe array (HFS-LP) is indicated by the blue line, with its uncertainty indicated by the blue shaded region. The other lines give the corresponding simulated profiles from the three codes. Extracted from[START_REF] Oliveira | Validation of edge turbulence codes against the TCV-X21 diverted L-mode reference case[END_REF].

Figure 8 :

 8 Figure 8: Grillix simulation for a 3D attached L-mode discharge in ASDEX-Upgrade at scale 1. Density fluctuations in [10 19 m 3 ]. Grid of about 7 millions points, dt = 5 × 10 -5 R 0 /C s = 1.2ns chosen slightly below the stability limit. The model includes a third order hyperviscosity. Extracted from [116].

Figure 9 :

 9 Figure 9: SOLEDGE3X simulation for a 3D attached L-mode discharge in TCV-X21 at scale 1. Density (left) and electron temperature (right) fluctuations. Grid of about 17 millions points, dt = 2ns chosen slightly below the stability limit. Extracted from [117].

Figure 10 :

 10 Figure 10: BOUT++ electromagnetic simulations for a 3D H-mode discharge in DIIID at scale 1 (only one-fifth of torus is simulated). Magnetic fluctuations (left) and outer-divertor heat-load profiles at different instants (mapped back to the outboard midplane). N x × N y × N z = 260 × 64 × 64.Extracted from[START_REF] Zhu | Electromagnetic turbulence simulation of tokamak edge plasma dynamics and divertor heat load during thermal quench[END_REF].

Figure 11 :

 11 Figure 11: GBS simulations of TCV with impurities showing a detached regime at the divertor inner target. Top Left: Time-and toroidally-averaged heat flux. Top right: Time-and toroidally-averagd ion particle flux. Bottom: Degree of Detachment DOD, for flux tube at different locations, in the high-density simulation. DOD > 1 corresponds to a detached plasma[START_REF] Loarte | Plasma detachment in jet mark i divertor experiments[END_REF]. Extracted and modified from[START_REF] Mancini | Self-consistent multi-component simulation of plasma turbulence and neutrals in detached conditions[END_REF].

Preprint submitted to Computers and FluidsJuly 2023

Acknowledgements

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 -EUROfusion). Views and 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.