Sebastian Minjeaud

Richard Pasquetti
email: richard.pasquetti@unice.fr

Fourier-spectral element approximation of the ion-electron Braginskii system with application to tokamak edge plasma in divertor configuration

Keywords: Magnetic confined fusion, Edge plasma, Braginskii closures, High order methods

Due to the extreme conditions required to produce energy by nuclear fusion in tokamaks, simulating the plasma behavior is an important but challenging task. We focus on the edge part of the plasma, where fluid approaches are probably the best suited, and our approach relies on the Braginskii ion-electron model. Assuming that the electric field is electrostatic, this yields a set of 10 strongly coupled and non-linear conservation equations that exhibit multiscale and anisotropy features. The computational domain is a torus of complex geometrical section, that corresponds to the divertor configuration, i.e. with an "X-point" in the magnetic surfaces. To capture the complex physics that is involved, high order methods are used: The time-discretization is based on a Strang splitting, that combines implicit and explicit high order Runge-Kutta schemes, and the space discretization makes use of the spectral element method in the poloidal plane together with Fourier expansions in the toroidal direction. The paper thoroughly describes the algorithms that have been developed, provides some numerical validations of the key algorithms and exhibits the results of preliminary numerical experiments. In particular, we point out that the highest frequency of the system is intermediate between the ion and electron cyclotron frequencies.

Introduction

Although the first investigations on magnetic confined fusion trace back to the early 50's, in Russia, researches on this topic remain of strong actuality with the present construction in Cadarache of the ITER tokamak. Thus, important efforts about the numerical simulation of the plasma behavior are worldwide carried out, and many models are used depending on the part of the plasma which is considered (the core or the edge), on the time scale which is focused on or on the particular objective which is aimed at. Among the modeling strategies, one may essentially discern kinetic (or gyrokinetic), multi-fluid, magnetohydrodynamic (MHD) or Grad-Shafranov approaches, see e.g. [START_REF] Jardin | Computational methods in plasma physics[END_REF] and references herein. For instance, a direct numerical simulation, see e.g. [START_REF] Grandgirard | a full-f global gyrokinetic Semi-Lagrangian code for ITG turbulence simulations[END_REF], needs a more accurate description of the evolution of the plasma than the one required for control considerations, see e.g. [START_REF] Blum | Reconstruction of the equilibrium of the plasma in a Tokamak and identification of the current density profile in real time[END_REF].

Because we are interested in short time scale turbulence phenomena in the scrape off layer (SOL), i.e. in the edge part of the plasma, the approach that we consider relies on the multi-fluid model that results from the moments of the Boltzmann equation expressed for the ions and the electrons, together with Braginskii like closures, as generally obtained using the Chapman-Enskog methodology [START_REF] Balescu | Transport processes in plasmas[END_REF][START_REF] Braginskii | Transport processes in a plasma[END_REF][START_REF] Goedbloed | Principles of Magnetohydrodynamics[END_REF][START_REF] Huba | NRL Plasma formulary, The office of naval research[END_REF][START_REF] Meier | A general nonlinear fluid model for reacting plasma-neutral mixtures[END_REF]. Fluid approaches are here relevant, because the temperature in the SOL is less high, i.e. generally of order of millions of degrees rather than hundreds of millions in the core, and also since the geometry may be much more complex. At this point, several variants are possible, leading to the MHD equations [START_REF] Anderson | TERPSICHORE: A threedimensional ideal MHD stability program[END_REF][START_REF] Czarny | MHD stability in x-point geometry : simulation of ELMS[END_REF][START_REF] Lütjens | The XTOR code for nonlinear 3D simulations of MHD instabilities in tokamak plasmas[END_REF] or to multi-fluid equations, and several simplifications are conceivable. For instance, in the codes BOUT, GBS or TOKAM3X [START_REF] Dudson | BOUT++: A framework for parallel plasma fluid simulations[END_REF][START_REF] Ricci | Simulation of plasma turbulence in the scrape-off layer conditions: the GBS code, simulation results and code validation, plasma[END_REF][START_REF] Tamain | Tokam-3d: A fluid code for tranport and turbulence in the edge plasma of tokamaks[END_REF], the so-called drift velocity approximation is used. It mainly consists in solving only for the components, of the ion and electron velocities, that are parallel to the magnetic field, thus allowing a drastic reduction in the number of unknowns. Other approaches, based on further simplifications of the physical model are also possible, see e.g. [START_REF] Bonnin | Code development for ITER edge modeling -SOLPS5.1[END_REF][START_REF] Buffereand | Near wall plasma simulation using penalization technique with the transport code SolEdge2D-Eirene[END_REF][START_REF] Michoski | Discontinuous Galerkin methods for plasma physics in the scrape-off layer of tokamaks[END_REF].

Like in [START_REF] Knight | CENTORI: A global toroidal electromagnetic two-fluid plasma turbulence code[END_REF][START_REF] Shumlak | Advanced physics calculations using a multi-fluid plasma model[END_REF], we do not use the velocity drift approximation, in order to rely on a more relevant model to describe plasma turbulence phenomena and also to provide a different and so complementary approach to more usual investigations. In a first step, we only consider electrically charged particles, i.e. the ions and the electrons, but a fluid modeling of neutral particles could be rather easily implemented in this framework. We also assume that at the short time-scales that we consider, the magnetic fluctuations can be neglected, i.e. the model makes use of the electrostatic assumption.

The resulting set of partial differential equations (PDEs) is solved using a high order numerical method in view to capture the instabilities that may develop in the plasma. Our numerical approach combines Fourier expansions in the toroidal direction with a spectral element method (SEM) in the poloidal plane. Using Fourier expansions more or less allows to address a set of two-dimensional problems rather the 3D one, which of course is efficient to reduce the computational time. In the frame of MHD and extended MHD, comparable approaches may be found elsewhere [START_REF] Guermond | Nonlinear magnetohydrodynamics in axisymmetric heterogeneous domains using a Fourier / finite element technique and an interior penalty method[END_REF][START_REF] Sovinec | the NIMROD team, Nonlinear magnetohydrodynamics simulation using high order finite elements[END_REF], and the codes NIMROD [START_REF] Sovinec | the NIMROD team, Nonlinear magnetohydrodynamics simulation using high order finite elements[END_REF] and M3D-C1 [START_REF] Ferraro | Calculations of two-fluid magnetohydrodynamic axisymmetric steady states[END_REF][START_REF] Jardin | A high-order implicit finite element method for integrating the two-fluid magnetohydrodynamic equations in two dimensions[END_REF] make also use of a high order approximation in the poloidal plane. The SEM is essentially a high order finite element method, that makes use in its 2D native form of unstructured quadrangular meshes, and of the tensorial product of Lagrange polynomials based on the Gauss-Lobatto-Legendre nodes, so that interpolation and quadrature points coincide [START_REF] Canuto | Spectral methods: Application to complex geometries and applications to fluid dynamics[END_REF][START_REF] Deville | High-order methods for incompressible flows[END_REF][START_REF] Karniadakis | Spectral hp element methods for CFD[END_REF]. The SEM approximation is well adapted to the complex geometries that are met when considering a tokamak in divertor configuration. Such a configuration is indeed characterized by the presence of an "X point" on the magnetic line (called the separatrix) that separates the closed and the open magnetic lines.

In Section 2, we go into the details of the physical model and outline the expected numerical difficulties. In Section 3, we describe the numerical method, which should be able to handle strongly anisotropic non self adjoint problems and also to support severe constraints on the time-step. Numerical validation exercises are provided in Section 4. They illustrate the high order convergence of our algorithm on toy problems that reflect the expected difficulties. In Section 5, we first define the initial and boundary conditions of the problem that we want to address, and especially focus on the non-trivial so-called Bohm condition that should be enforced at the targets. Preliminary numerical results are then provided and discussed. They point out the stabilizing effect of the Braginskii closures and that the highest frequency of the system is intermediate between the ion and electron cyclotron frequencies, consistently with the theoretical study presented in Appendix C. We conclude in Section 6.

The physical model

The physical model is essentially a two-fluid ion-electron model as may obtained from the moments of the Boltzmann equation together with the Braginskii closures. The physical assumptions we make to obtain this system of PDEs are given in Section 2.1. The complete set of equations is presented in Section 2.2, whereas the Braginskii closures are recalled in Section 2.3. Finally, in Section 2.4, we describe some of the difficulties associated to the present model.

Framework and physical assumptions

The modeling relies on the following assumptions:

(i) A two-fluid ion-electron modeling, based on the equations of conservation of density, momentum and energy for each species, together with the Braginskii closures, is relevant. (ii) The plasma is electrically neutral:

s n s e s = 0 (n s : species density, e s : species charge, s ∈ {i, e}). (iii) The magnetic field B is axisymmetric and its fluctuations in time are negligible with respect to the magnetic field imposed by the coils. (iv) The electric field is electrostatic: E = -∇U (E: electric field, U : electric potential).

(v) The species temperature T s and pressure p s are linked by state laws of the form p s = n s T s and one has ε s = p s /(γ -1) (ε s : internal energy, γ: ratio of specific heats).

Note that the assumptions (iii) and (iv) are coherent with the Maxwell-Faraday equation : ∇×E = -∂ t B. More precisely, the assumption (iv) may be deduced from the assumption (iii) and the Maxwell-Faraday equation since they give: ∇ × E = 0.

In addition to the coupling due to the Braginskii closures, the assumptions (ii) and (iv) induce a strong coupling between the conservation equations for ions and for electrons. The electric potential U may be seen as an additional unknown associated to the constraint of electric neutrality, which translates into a divergence free constraint on the current. Indeed, using the particles conservation equations:

∂ t n s + ∇ • (n s u s ) = 0, (1)
where u s is the species velocity, by first multiplying the ion and the electron conservation equation by the ion and electron electric charge, respectively, and then summing, one finds that ensuring the electroneutrality is equivalent to ensure that the current, j = s n s e s u s , is divergence free: ∇ • j = 0.

The complete set of equations

Let us now define the species mass densities ρ s = n s m s (m s : species mass), the species momentums q s = ρ s u s and let us set w s = e s /m s . Due to the electroneutrality assumption and for the binary configuration we consider, the ion mass density ρ i and the electron mass density ρ e are simply proportional to the total mass density ρ = s ρ s , i.e. ρ s = α s ρ, and one can check that α s = (1 -(w i /w e ) sgn(es) ) -1 (sgn stands for the sign function). The PDEs system that we consider may then write:

∂ t ρ + ∇ • (q i + q e ) = 0, ∂ t q s + ∇ • (q s ⊗ u s + p s I + Π s ) = -w s α s ρ∇U + w s q s × B + R s , ∇ • (w i q i + w e q e ) = 0, ∂ t ε s + ∇ • (ε s u s + ϕ s ) = -p s ∇ • u s -Π s : ∇u s + Q s .
The terms not yet defined are those for which the Braginskii closure is needed, namely: the collision terms R s and the viscous stress tensors Π s in the momentum conservation equations, the collision terms Q s and the heat flux densities ϕ s in the internal energy equations. These terms are all defined in Section 2.3. Such a PDEs system shows ten non-linear and coupled scalar PDEs and ten unknown scalar fields. Of course, the equations for the internal energies, ε s , yield evolution equations for the pressures, p s , by multiplication by (γ -1).

Braginskii closures

The Braginskii closures were proposed in the celebrated paper [START_REF] Braginskii | Transport processes in a plasma[END_REF]. Here we simply provide expressions that are not associated to any specific coordinate system, especially the local coordinate system such that one of the axis is aligned on the magnetic field. These expressions are those used in our numerical implementation. It should be noticed that variants or alternatives to the Braginskii model do exist, see e.g. [START_REF] Balescu | Transport processes in plasmas[END_REF][START_REF] Goedbloed | Principles of Magnetohydrodynamics[END_REF]. In practice, it would not be an heavy task to implement them. Collision terms. The collision terms that appear in the momentum equations compensate. They write R e = R and R i = -R, with:

R = R u + R T , R u = - m e n e τ e (0.51u + u ⊥ ), R T = -0.71n e ∇ T e -1.5 n e ω ce τ e b × ∇T e ,
where u = u e -u i . The collision terms in the energy equations write:

Q i = 3 m e m i n e τ e (T e -T i ) Q e = -Q i -R • u
Stress tensors. The stress tensors decompose in three terms:

Π s = -η Π -η ⊥ Π ⊥ + η ∧ Π ∧ ,
with, omitting the subscript s for the sake of simplicity:

Π = 3 2 b ⊗ b - 1 3 I b ⊗ b : W , Π ⊥ = 1 4 W + 3(W b ⊗ b + b ⊗ bW ) + 1 2 (I -15b ⊗ b)(b ⊗ b : W ) , Π ∧ = 1 4 M b W (I + 3b ⊗ b) -(I + 3b ⊗ b)W M b ,
where W = ∇u s + (∇u s ) t -2/3∇ • u s I is the strain rate tensor typical of Newtonian fluids and M b is the antisymmetric matrix such that for any vector u, M b u = b × u. The coefficients write, for the ions:

η = 0.96n i T i τ i , η ⊥ = 1.25 η (ω ci τ i ) 2 , η ∧ = 1.04 η ω ci τ i ,
and for the electrons :

η = 0.73n e T e τ e , η ⊥ = 2.79 η (ω ce τ e ) 2 , η ∧ = -1.37 η ω ce τ e .
Thermal fluxes. Like the stress tensors, the ion thermal fluxes show three components:

ϕ i = -κ i ∇ T i -κ i ⊥ ∇ ⊥ T i + κ i ∧ b × ∇T i ,
where :

κ i = 3.9 n i T i τ i m i , κ i ⊥ = 0.51 κ i (ω ci τ i ) 2 , κ i ∧ = 0.64 κ i ω ci τ i .
The electron thermal fluxes write, again with u = u e -u i : (ω ce τ e ) 2 , κ e ∧ = 0.78

ϕ e =
κ i ω ci τ i .

Expected difficulties

Let us underline some difficulties associated to the present modeling:

• The PDEs system clearly associates equations that usually correspond to description of compressible flows on the one hand, and to the description of incompressible fluids on the other. This is especially obvious when deriving equations for the total momentum and for the current, see [START_REF] Bonnement | Towards a Fourier-SEM solver of fluid models in tokamaks[END_REF]. The governing equations for the momentum are then typical of compressible flows, whereas those for the current are typical of incompressible ones, with the electric potential associated to the divergence free constraint.

• The problem exhibits some multi-scale features. As well known, the particles trajectories show a spiral motion around the magnetic lines, which "Larmor radius" is very small with respect to a characteristic size of the plasma, see e.g. [START_REF] Chen | Introduction to plasma physics and controlled fusion. 1. Plasma physics[END_REF][START_REF] Delcroix | Physique des plasmas[END_REF]. This multi-scale feature in space is associated to a multiscale feature in time, because the "cyclotron frequency" is much smaller than e.g. the discharge time.

• The electrons are about 2000 times lighter than the ions. In the frame of MHD approaches, this justifies the generalized Ohm law, which is obtained when neglecting the inertial terms in the electron momentum equation.

• The magnetic field is very strong in tokamaks. This leads, in our numerical algorithm, to a strong anisotropic problem on the electrical potential U , see Section 3.1.2.

• The Braginskii closures introduce a complex nonlinear coupling between the equations together with strong anisotropy phenomena.

-The diffusion phenomena are essentially aligned on the magnetic lines, which is mainly associated to the fact that the magnetic field is strong. From the expressions given in Section 2.3, one can observe that

η ⊥ η ∝ κ ⊥ κ ∝ 1 (ω c τ ) 2 , η ∧ η ∝ κ ∧ κ ∝ 1 ω c τ .
Some values, typical of the core and of the edge parts of the plasma, of the anisotropy factor ω c τ are provided in Table 1.

Plasma core Plasma edge T (K)

1. [START_REF] Günter | Modelling of heat transport in magnetised plasmas using non-aligned coordinates[END_REF] 2, for both the ions and the electrons and for two different temperature levels typical of the plasma edge. Such large values may however be nonphysical because if increasing too much the temperature and thus the mean free path, the plasma becomes nearly collisionless and so the Braginskii closures are less relevant. This issue is e.g. addressed in [START_REF] Epperlein | Plasma transport coefficients in a magnetic field by direct numerical simulation of the Fokker-Planck equation[END_REF]. From numerical point of view, large values of κ and η imply strong stability constraints within an explicit time marching, since the time-step should scale as their inverses.

T

Numerical method

The time discretization is described in Section 3.1 and the formulation of the space discretization is explained in Section 3.2. The anisotropy feature is discussed in Section 3.3. More specific details about the numerical algorithm are given in Section 3.4.

Time scheme

Due to the strong couplings and non-linearities of the PDEs system, one may prefer an explicit timescheme rather than an implicit one. However, the Lorentz forces should be treated implicitly to avoid numerical instabilities. Indeed, the approximation with an Euler forward scheme of the Ordinary Differential Equation (ODE): ∂ t q = q × B is unstable, since the exact solution of the ODE is such that:

q • ∂ t q = 0 , |q| = Constant b • ∂ t q = 0 , q = Constant whereas the Euler forward approximation yields: (q n+1 ⊥ ) 2 = (q n ⊥ ) 2 + (δt|B||q n ⊥ |) 2 ,
with δt for the time step. Another source of difficulty is the divergence free constraint on the current. In [START_REF] Bonnement | Towards a Fourier-SEM solver of fluid models in tokamaks[END_REF], we suggested to use an implicit-explicit Runge Kutta method (IMEX RK [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF]) for the PDEs system of conservation equations, associated to a projection method, at the end of each RK step, for the divergence free constraint. On the basis of the Helmholtz decomposition of any (sufficiently smooth) vector field into gradient and rotational parts, the projection step consisted in, first, computing the electric potential that allows to enforce the divergence free constraint and, then, updating in consequence the ion and electron momentums. However, it has turned out that such a rather natural approach was not satisfactory, yielding time-step dependent results. Some explanations on a simplified model are provided in Appendix C.2. To overcome this difficulty, hereafter we suggest to use a Strang splitting scheme that allows to address the Lorentz force and the divergence free constraint consistently in the same implicit sub-step.

The splitting scheme

Let us recall that, for the PDE: ∂ t u + L(u) = 0, with L ≡ L 1 + L 2 , the Strang splitting consists in solving, with n for the time index:

∂ t u + L 1 (u) = 0 , t ∈ (t n , t n+1/2 ), u(t n ) u n -→ ũn+1/2 u(t n+1/2 ) , ∂ t u + L 2 (u) = 0 , t ∈ (t n , t n+1 ), u(t n ) ũn+1/2 -→ u n+1/2 u(t n+1 ) , ∂ t u + L 1 (u) = 0 , t ∈ (t n+1/2 , t n+1 ), u(t n+1/2 ) u n+1/2 -→ u n+1 u(t n+1 ) .
The splitting that we suggest makes use of the structure above. We consider for the System I:

∂ t ρ + ∇ • (q i + q e ) = 0, ∂ t q s + ∇ • (q s ⊗ u s + p s I + Π s ) = R s , ∂ t p s + ∇ • (p s u s + (γ -1)ϕ s ) = (γ -1)(-p s ∇ • u s -Π s : ∇u s + Q s ), (System I)
and we include the Lorentz force and the divergence free constraint in the System II:

∂ t q s = -w s α s ρ∇U + w s q s × B, ∇ • (w i q i + w e q e ) = 0 . (System II)
To integrate each system we suggest to use RK schemes with large absolute stability regions: An explicit RK (ERK) scheme is used for System I, whereas a Diagonally Implicit RK (DIRK) scheme is used for System II. The explicit treatment is chosen to avoid the full coupling of the PDEs, since this would require the resolution of a huge and non-linear algebraic system at each time-step. However, as mentioned in Section 2.4, if treated explicitly the diffusion terms imply strong stability constraints on the time-step δt. To overcome this difficulty without impacting dramatically the computational time, a sub-time cycling technique is implemented to integrate the System I with a sub-time-step δt/s, where s is the number of subcycles.

The Strang splitting is second order accurate as soon as the schemes used for each step are at least second order accurate, which is the case in our implementation by using RK schemes of order 3. However, it is worth noting that System II will be supplemented with boundary condition to ensure it is well posed. This may lead to enforce some artificial boundary conditions on the electric potential, possibly resulting in some loss of accuracy of the potential approximation. This phenomenon is well known for projection method in incompressible fluid mechanics, which approximates the pressure at an order lower than for the velocity, see e.g. [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF].

Resolution of System II. Computation of the potential

The System II can be recasted to an elliptic problem on the electric potential. At each DIRK step, one has indeed to compute q i , q e and U such that:

q s = qs + aδt(-w s α s ρ∇U -w s BM b q s ) , s ∈ {i, e} ∇ • s w s q s = 0,
where qs is the previous RK estimate of q s , a the diagonal coefficient (Butcher tableau) of the DIRK scheme, δt the time-step and M b the antisymmetric matrix introduced previously such that b × q s = M b q s . Then, expressing the current and taking its divergence yields:

∇ • A∇U = s w s ∇ • A s qs , (2)
where A = ρ s A s a δt w 2 s α s and A s = (I + a δt w s BM b ) -1 . A calculation, postponed to Appendix A, shows that the matrix A s can be expressed as follows:

A s = (I -γ s BM b + γ 2 s B ⊗ B)/(1 + γ 2 s B 2 ), with γ s = a δt w s . (3)
Then, it readily follows from the equality:

x, x + γ s x × B + γ 2 s B(B • x) = (x, x) + γ 2 s (B • x) 2
, that A s is coercive and, since aδtw 2 s α s ρ > 0, A is also coercive. Thus, U solves a non self adjoint elliptic equation which is well posed, from the Lax-Milgram theorem, in the Sobolev space H 1 if Dirichlet conditions are considered or in H 1 /R with admissible Neumann conditions.

Once the potential U is known, one can compute the ion and electron momentums from:

q s = A s (q s -aδtw s α s ρ∇U ) . (4)
The following remarks can be expressed:

• The matrices A s depend on the time-step and on the RK step number, through the coefficient a, except if the DIRK scheme is such that its diagonal coefficient (in the Butcher tableau) is constant (this is the case in our implementation).

• The problem (2) exhibits an anisotropy structure, see Section 3.3.

• A natural boundary condition that yields a well posed problem, is obtained by applying the divergence theorem to both sides of (2) and stating that the equality is verified everywhere on the boundary

(A∇U - s w s A s qs ) • n = 0, ( 5
)
where n is the unit outwards vector normal to the boundary Γ of the computational domain. By summing the equation ( 4), expressed for the ions and the electrons and weighted by the coefficients w s , one easily checks that the boundary condition ( 5) is equivalent to the condition j • n = 0, i.e. an impermeability condition for the current, which of course is consistent with the divergence free constraint.

Fourier-SEM spatial approximation

The tokamak is represented in cylindrical coordinates by the domain Ω = Ω × [0, 2π], with (r, z) ∈ Ω and θ ∈ [0, 2π]. The Fourier-SEM method consists of using a Fourier approximation in the toroidal direction θ and the SEM in the poloidal plane Ω. Hence, let us introduce an integer K θ which stands for the number of grid points (or rz-planes) used for the discretization in the toroidal direction, an unstructured quadrangular mesh of the poloidal plane Ω and an integer N which stands for the polynomial approximation degree in each element. We denote by ϕ m (r, z) the SEM basis function associated to the grid point of index m in Ω. Recall that, on the reference element, SEM basis functions are obtained by tensorial product of the Lagrange polynomials (of degree N ) based on the N + 1 Gauss-Lobatto-Legendre 1D nodes.

For the sake of simplicity we first focus on the toy PDE:

∂ t u + ∇ • f = s,
where f is typically a diffusion flux term, to which an integration by parts is to be applied. The advection flux term is assumed included in the source term s. The weak formulation of this PDE writes:

Ω v∂ t u dΩ + Γ vf • dΓ - Ω f • ∇v dΩ = Ω vs dΩ.
In the Fourier space, u k ûk (r, z)e ikθ , with |k| ≤ K θ /2 and i 2 = -1. The ûk and û-k are complex conjugate since u is a real function. It is natural to use, in each spectral element, basis functions of the form ϕ m (r, z)e ilθ . Then, since (i) in the cylindrical coordinate system dΩ = rdθd Ω and (ii) the L 2 hermitian product, in θ ∈ (0, 2π), of e ikθ and e ilθ equals 2πδ k-l,0 (δ, Kronecker symbol), one obtains the expression that should hold for any Fourier mode k and any SEM basis function ϕ m (r, z):

Ω r ϕ m ∂ t ûk d Ω + Γ r ϕ m f k • d Γ - Ω r f k • ∇ϕ m d Ω = Ω r ϕ m ŝk d Ω ,
where, if using the usual dot product:

∇ϕ m = ( ∂ r ϕ m , -ikϕ m /r, ∂ z ϕ m ).
Things are slightly more complicated in the vectorial case

∂ t u + ∇ • f = s,
where f is now a second order tensor. After multiplication by the vector test function v and integration by parts, one finds:

Ω v • ∂ t u dΩ + Γ v • f dΓ - Ω f : ∇v dΩ = Ω v • s dΩ.
Taking into account the cross terms in the divergence operator, when expressed in the cylindrical coordinate system, and specializing the components of v to the r, θ and z components of the PDE one obtains:

Ω v r ∂ t u r dΩ + Γ v r (f dΓ) r = Ω (f r. • ∇v r + f θθ v r /r) dΩ + Ω v r s r dΩ, Ω v θ ∂ t u θ dΩ + Γ v θ (f dΓ) θ = Ω (f θ. • ∇v θ -f rθ v θ /r) dΩ + Ω v θ s θ dΩ, Ω v z ∂ t u z dΩ + Γ v z (f dΓ) z = Ω f z. • ∇v z dΩ + Ω v z s z dΩ,
where for example, f r. = (f rr , f rθ , f rz ). The corresponding expressions in Fourier space follow easily.

Finally, one should address the elliptic potential equation [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF], that may write, with obvious notation:

∇ • A∇U = ∇ • s.
In weak form and assuming the natural boundary condition [START_REF] Bonnin | Code development for ITER edge modeling -SOLPS5.1[END_REF], one obtains:

Ω A∇U • ∇v dΩ = Ω s • ∇v dΩ .
One notices that tensor A only depends on θ via the density ρ. Let us assume, in a first step, that ρ(r, θ, z) ≈ ρ0 (r, z), then one obtains for each Fourier mode k:

Ω rA ∇ Ûk • ∇ϕ m d Ω = Ω rŝ k • ∇ϕ m d Ω.
Thus, we have substituted a set of 2D like problems to the initial 3D one. The methodology is then standard, setting up a stiffness matrix, the right hand side and then solving. However, on the contrary of what will be obtained for a Poisson equation, it should be noticed that the real and imaginary parts of each mode are coupled by the non-diagonal terms of the tensor A.

The space variations of ρ mainly occur in the poloidal plane, so that one only expects small fluctuations in the toroidal direction. To take them into account one may think to use a fixed point procedure, but this remains to be included in our numerical implementation.

Anisotropic diffusion

As already mentioned, the Braginskii closures yield strongly anisotropic problems. When using a mesh non-aligned on the magnetic field, difficulties are then generally observed [START_REF] Van Es | Finite-difference schemes for anisotropic diffusion[END_REF][START_REF] Günter | Modelling of heat transport in magnetised plasmas using non-aligned coordinates[END_REF][START_REF]Three discontinuous Galerkin schemes fo rthe anisotropic heat conduction equation on non-aligned grids[END_REF]. Typically, one should handle expressions of the form:

ϕ = -κ ∇ T -κ ⊥ ∇ ⊥ T + κ ∧ (b × ∇T ) .
To this end, we introduce a diffusion tensor such that, see [START_REF] Bonnement | Towards a Fourier-SEM solver of fluid models in tokamaks[END_REF] for details:

ϕ = -K∇T , K = (κ -κ ⊥ ) b ⊗ b + κ ⊥ I + κ ∧ M b
that can be easily handled thanks to the weak formulation on which relies the SEM.

For the sake of completeness, let us mention that [START_REF] Bonnement | Towards a Fourier-SEM solver of fluid models in tokamaks[END_REF]:

• The ∇ • κ ∧ (b × ∇T ) term behaves like a transport term with velocity

u ∧ = ∇κ ∧ × b + κ ∧ ∇ × b
It is thus associated to curvature effects or to the inhomogeneity of the temperature dependent κ ∧ coefficient.

• Using the Fourier-SEM approach, satisfactory numerical results have been obtained. Especially, as in [START_REF] Fischer | Anisotropic diffusion in a toroidal geometry[END_REF] we have considered a test-case proposed by the Center for Extended Magnetohydrodynamics Modeling (CEMM), namely the "evolution of a Gaussian pulse". For this unsteady diffusion problem in a torus and using a mesh non aligned on the magnetic field lines, very satisfactory results have been obtained for the extreme case κ ⊥ /κ = 0. Moreover, comparisons carried out with P 1 finite elements have outlined the importance of using a high order method, if not using a mesh aligned on the magnetic lines [START_REF] Bonnement | Modélisation numérique par approximation fluide du plasma de bord des tokamaks (projet ITER)[END_REF].

At this point, it is of interest to point out that independently of the Braginskii closures, the elliptic PDE that the potential solves, see [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF], also shows a strong anisotropy. Indeed, the diffuse flux A∇U can be recasted as follows:

A∇U = A ∇ U + A ⊥ ∇ ⊥ U + A ∧ (b × ∇U ), with A = ρa δt s w 2 s α s , A ⊥ = ρa δt s w 2 s α s 1 + γ 2 s B 2 , A ∧ = ρa δt s γ s w 2 s α s 1 + γ 2 s B 2 γ s .
This points at the anisotropy ratio:

A ⊥ A = s w 2 s α s 1 + γ 2 s B 2 s w 2 s α s .
Clearly, this ratio depends on the time-step through the coefficients γ s . In the limit δt = 0, one has A ⊥ /A = 1. More interesting, one can estimate this ratio for the following extreme values of the time-step: δt = 2π/ω ci and δt = 2π/ω ce , that are associated to the ion and electron cyclotron angular frequencies. We find, in the first case A ⊥ /A ≈ 5.01 10 -5 and in the second case A ⊥ /A ≈ 9.24 10 -2 . Hence, for the values of the time step used in the simulations, the potential solves a strongly anisotropic elliptic PDE.

Computational details

• In order to evaluate in the most accurate way the advective and diffusive fluxes, it is a priori preferable to formulate them in terms of the computational variables, i.e. the density and the ion and electron momentum and pressures. Thus, from p s = ρ s T s and q s = ρ s u s , one obtains for the temperature and velocity gradients:

∂ j T s = T s ∂ j p s p s - ∂ j ρ ρ , ∂ j u s i = 1 ρ (∂ j q s i /α s -u s i ∂ j ρ).
Let us recall that in the cylindrical coordinate system the velocity gradient tensor shows cross terms:

∇u s =   ∂ r u s r 1 r (∂ θ u s r -u s θ ) ∂ z u s r ∂ r u s θ 1 r (∂ θ u s θ + u s r ) ∂ z u s θ ∂ r u s z 1 r ∂ θ u s z ∂ z u s z   ,
and that its trace equals the divergence of the velocity.

Concerning the advective fluxes one uses:

∇ • p s u s = u s • ∇p s + p s ∇ • u s , ∇ • (q s ⊗ u s + p s I) = ∇q s u s + ∇ • u s q s + ∇p s .
• To solve the algebraic system associated to the potential, we use a static condensation technique: The unknowns associated to the inner grid-points of each element are eliminated, so that one only solves for the unknowns located on the edges of the elements. The algorithm is the one described in [START_REF] Lazar | Fekete-Gauss spectral elements for incompressible Navier-Stokes flows: The two-dimensional case[END_REF].

The system is then of smaller dimension, i.e. O(N ) rather than O(N 2 ) (N , polynomial approximation degree) and better conditioned, i.e. O(N ) rather than O(N 3 ) [START_REF] Melenk | On condition numbers in hp-FEM with Gauss-Lobatto based shape functions[END_REF].

• A BICGSTAB algorithm is used to solve this non-symmetric algebraic system for the potential. In a matrix free implementation of the method, we are limited to the use of a Jacobi preconditioner. When assembling the stiffness matrix (in sparse manner), the best results have been obtained with the algebraic multigrid preconditioner HYPRE used through the PETSc software. It requires only a few iterations of the BICGSTAB algorithm.

• In the explicit part, a spectral vanishing viscosity (SVV) technique is used for stabilization, see e.g. [START_REF] Xu | Stabilized spectral element computations of high-Reynolds number incompressible flows[END_REF]. This allows to address stiff problems together with preserving the "spectral accuracy" of the Fourier-SEM approximation.

• Isoparametric mappings may be used to take into account the curvature, in the poloidal plane, of the boundary of the computational domain. The technique is based on a bending procedure similar to the one proposed in [START_REF] Lazar | Fekete-Gauss spectral elements for incompressible Navier-Stokes flows: The two-dimensional case[END_REF] for triangular elements.

Numerical validation

In this Section we focus on the key algorithms of the present Fourier-SEM method and provide numerical validations by comparison to exact solutions. First, the accuracy of the method is quantified for an anisotropic diffusion problem inspired from a test case proposed by the CEMM. Second, the efficiency of the algorithm proposed for the resolution of System II is checked for the toy problem described in Appendix C.

Anisotropic diffusion with exact solution

The computational domain is a torus with axis at (r = r 0 , z = 0). Given a function ψ(r, z) that only depends on the distance to the torus axis, as in [START_REF] Fischer | Anisotropic diffusion in a toroidal geometry[END_REF] we define the magnetic field as:

B = ∇ × ( ψ r e θ ) + 1 r e θ = 1 r (∂ z ψe r + e θ -∂ r ψe z ) , (6)
and we choose: ψ = C ψ ln(1 + ((r -r 0 ) 2 + z 2 )/a 2 ), C ψ = a 2 /(2r 0 q 0 ) with r 0 = 3 and a = q 0 = 1. The resulting field lines are then helically wrapped on toroidal surfaces, i.e. of constant distance to the torus axis.

In [START_REF] Fischer | Anisotropic diffusion in a toroidal geometry[END_REF] tests were carried out for a steady axisymmetric diffusion problem, using the diffusion coefficients κ = 1 and κ ⊥ = 10 -9 .

Here we consider the extreme anisotropic case such that κ ⊥ /κ = 0, i.e. with diffusion tensor K = κ B ⊗ B/B 2 . Such a problem is generally ill posed in the steady case, but is well posed in the unsteady one:

∂ t u + ∇ • (K∇u) = s.
Clearly, if s = 0 and if the initial condition is constant on the toroidal surfaces, then, if assuming coherent boundary conditions, the solution of the diffusion problem remains constant. One can make this test case more convincing by choosing an appropriate source term s so that the exact solution is unsteady and non axisymmetric. Details are provided in Appendix B. Computations have been carried out for κ = 1 in the torus of cross section (3, 4) × (-1, 1) and in the time interval (0, 1), with exact solution and source term:

u(r, θ, z, t) = exp -2 (r -r 0 ) 2 + z 2 sin θ sin πt 2 , s = - κ B 2 r 4 sin πt 2 ∂ 2 θ u + π 2 cos πt 2 u. (7)
This exact solution is used to define the initial condition and Dirichlet boundary conditions. The discretization of the square section is uniform and makes use of n el ∈ {12 2 , 18 2 , 24 2 } elements with polynomial approximation of degree N ∈ {2, 4, 6, 8} in each of them. In toroidal direction there is no Fourier approximation error since the variation of u is sinusoidal (in practice the θ-step equals π/2). The time-step δt = 2 10 -5 is small, so that the RK4 approximation error is negligible. For n el = 18 2 , Fig. 1 (left) shows the evolution of the error (in the max norm) for the different values of N . Their variations, at the final time t = 1, with respect to the polynomial degree N and for the different values of n el are shown in Fig. 1 (right). A semi-log plot is used to clearly demonstrate the exponentially fast decrease of the error. Thus, although the mesh is not aligned on the magnetic field lines very accurate results are obtained if using a high order approximation. One also observes that the expected increase in time of the error can be made very slow. In order to validate the algorithm used to approximate System II, we provide numerical simulations in a simplified geometry (slab configuration with a constant magnetic field) and with a constant density ρ, so that the exact solution of the problem can be analytically computed, see Appendix C.1. In Cartesian coordinates, denoting by (e x , e y , e z ) the canonical basis of R 3 , we assume that the magnetic field B = Be y is constant and aligned with the y-axis, and that the initial data (ion and electron momentums q i,0 , q e,0 ) belong to the (x, z)-plane. In this case, System II ensures that q i , q e ∈ span{e x , e z } for all time and the 3D problem is reduced to the following 2D one:

Numerical approximation of System II

       ∂ t q (x) s = -w s α s ∂ x Ψ -w s Bq (z) s , s ∈ {i, e}, ∂ t q (z) s = -w s α s ∂ z Ψ + w s Bq (x) s , s ∈ {i, e}, ∂ x w i q (x) i + w e q (x) e + ∂ z w i q (z) i + w e q (z) e = 0.
where q (x) s and q (z) s denotes the two components of the ion or electron momentum and Ψ stands for ρU . We supplement the system with the boundary condition

w i q (x) i
+ w e q (x) e n (x) + w i q (z) i + w e q (z) e n (z) = 0, where (n (x) , n (z) ) is the outward normal to the computational domain Ω. The computational domain is the unit square Ω = (0, 1) 2 . We assume that e i = -e e and m i = 1000m e . The initial data are defined from the initial total momemtum q 0 and the initial current j 0 by

q i,0 = 1 w i -w e
-w e q 0 + j 0 , q e,0 = 1 w i -w e w i q 0 -j 0 , and we use the following expressions for q 0 and j 0 q

(x) 0 = ∂ z Ψ 0 B , q (z) 0 = - ∂ x Ψ 0 B , j (x) 0 = 2 sin(πx) cos(πz), j (z) 0 = -2 cos(πx) sin(πz),
with Ψ 0 = Be -3(x-0.5) 2 -9(z-0.5) 2 . We can easily check that these initial data are compatible with the divergence free constraint and the boundary condition. The exact solution is given in Appendix C.1.

We perform numerical simulations using different uniform meshes with n el ∈ {4 2 , 8 2 , 16 2 , 24 2 }, different degrees of polynomial approximation N ∈ {2, 4, 6, 8} and different time steps δt = 2π m e /m i /B/n δt with n δt ∈ {80, 240, 720}. The variations of the max norm of the error for q (x) e , at the final time t = 2π m e /m i /B, with respect to the polynomial degree N are shown in Fig. 2, for the different values of n el and for n δt = 240, 720. We obtain similar results for all the variables. The errors decrease when N increases until they reach a saturation value due to the time discretization. A semi-log plot is used to emphasize the exponentially fast decrease of the errors with respect to N . Obviously, the saturation values are the same for all meshes and depend only on the time step. They are 1.57 10 -4 , 5.83 10 -6 and 2. 16 10 -7 for the time steps obtained with n δt = 80, 240, 720, respectively. It demonstrates that the method is of order 3 in time, consistently with the third order DIRK scheme used for the time integration.

Preliminary investigations

We consider the JET tokamak which is located at Culham, UK, and operated by the Culham Centre for Fusion Energy. In Section 5.1, we describe precisely the physical problem we solve, that is the geometry of the domain and of the magnetic field, the reference values used to get a dimensionless model, and the initial and boundary conditions. Some numerical results are given in Section 5.2. Our aim is not to carry out a realistic simulation which would require to use a parallel and optimized version of our code, but to illustrate the influence of the Braginskii closure and underline that the highest frequency of the system is intermediate between the ion and electron cyclotron frequencies, consistently with the theoretical study presented in Appendix C.

The physical problem 5.1.1. Geometry and magnetic field

The geometry of the edge part of the plasma, in the poloidal plane, together with a mesh aligned on the magnetic surfaces, was provided to us in the frame of collaborations with partners working on the MHD JOREK code [START_REF] Czarny | MHD stability in x-point geometry : simulation of ELMS[END_REF]. Fig. 3 shows the computational domain and the mesh, which is composed of quadrangular elements. Using such a mesh, it is possible to implement the SEM using any polynomial approximation degree N > 1 in each element.

The magnetic field is also shown in Fig. 3. It is obtained with the JOREK software when solving the Grad-Shafranov equilibrium ∇p = j × B. Such a computation provides the density ρ, the (total) pressure p, and the magnetic potential ψ. For the poloidal component B p , one defines

B p = 1 r -∂ z ψ, 0, ∂ r ψ ,
and one assumes, for the toroidal component: B θ = C/r (r: dimensionless distance at the axis of the tokamak, C: a constant equal to 3 in our simulations).

Boundary conditions

At the internal surface of the computational domain, i.e. against the plasma core, we assume Dirichlet conditions for all variables, except the potential. They will be supposed constant in time and so defined by the initial conditions. Elsewhere, no boundary conditions are enforced for the density since the fluid is supposed to be only outgoing.

At the outer surface, excluding the targets, we assume free-slip like conditions for the ion and electron momentums. Thus, in the explicit part of the algorithm (System I, see Section 3), we enforce these vector fields to be tangent to the boundary. At this part of the outer surface, natural boundary conditions, i.e. homogeneous Neumann conditions, are used for the pressures.

At the targets, i.e. at the remaining parts of the outer surface where the magnetic lines impact the boundary, we use the so-called Bohm boundary conditions. In their simplest form, this means that one should have M ≥ 1, where M is the ion parallel Mach number. More elaborate conditions are e.g. provided in [START_REF] Loizu | Boundary conditions for plasma fluid models at the magnetic presheath entrance[END_REF][START_REF] Stangeby | The plasma boundary of magnetic fusion devices[END_REF]. On the basis of investigations carried out on a one-dimensional "minimal transport model" [START_REF] Isoardi | Penalization modeling of a limiter in the tokamak edge plasma[END_REF], we control the ion velocity by varying the density and preserve the ion and electron momentums [START_REF] Minjeaud | High order approximation of a tokamak edge plasma transport minimal model with Bohm boundary conditions[END_REF]. Moreover, in this simplified frame we also assume that the temperatures are preserved, so that the pressures vary like the density. Thus, with (c, sound velocity):

M 2 = (u i • b) 2 c 2 = (q i • b) 2 γα i ρ 2 s α s T s /m s ,
the algorithm is the following: If M < 1, then set ρ := ρ/M and p s := p s /M .

Initial conditions

To set up the initial conditions, we again start from a Grad-Shafranov equilibrium obtained with the JOREK software. The difficulty comes from the fact that one has to go from single fluid (axisymmetric) results, say the density ρ and the total pressure p, to two-fluid ones. Taking into account the electroneutrality assumption s n s e s = 0 yields: ρ = s n s m s = n i (m i -m e e i /e e ) and thus the value of n i . Then, knowing the total pressure p = s p s and assuming equal the ion and electron "temperatures", T i and T e , one obtains p = s n s T s = n i T i (1 -e i /e e ), which yields the value of T i (= T e ). Visualizations of the density and of the temperature at the initial time are given in Fig. 4. To compute the momentums at t = 0, one first defines the current. When taking into account the axisymmetry of the magnetic field B and the assumption B θ = C/r, from the Ampère theorem one finds that the current j is azimuthal:

j ∝ ∇ × B = (∂ z B r -∂ r B z )e θ ,
so that ∇ • j = 0, i.e. j is identically divergence free, and j • ∇p = 0, as desired. Then, we compute j(t = 0), in such a way that the Grad-Shafranov equilibrium relation ∇p = j × B is at best verified. From the minimization, at each grid-point, of the Euclidean norm of the residual one obtains:

j θ = 1 B 2 p (B z ∂ r p -B r ∂ z p) .
Then, we define the ion velocity, assuming : (i) u i = 0 inside the separatrix and (ii) u i colinear to the magnetic field in the SOL:

u i = u b.
The parallel velocity u is defined by enforcing M = ±1 at the targets, and by assuming linear variations along the magnetic lines. One can then compute the ion momentum q i . The electron momentum q e is finally computed from the expression j = s w s q s .

Numerical results

To point out the influence of the Braginskii modeling, we compare the results obtained in the Euler approximation (simulation S0) to those obtained using the Braginskii closures (simulations S1 and S2). To make easier the computation, (i) as e.g. in [START_REF] Shumlak | Advanced physics calculations using a multi-fluid plasma model[END_REF] we use the heavy electron approximation m i /m e = 100, and (ii) the parallel transport coefficients of the ion and of the electron have been divided by 100 and by 500, respectively, in simulation S1, and by 10 and 50, in simulation S2.

Using the sequential version of the code, calculations have been done in an axisymmetric configuration, using 3417 spectral elements with polynomial approximation in each variable of degree N = 3. Then, the number of grid-points equals 31194 and one has 10 unknowns per grid-point. The dimensionless time-step is taken equal to δt = 5 10 -6 (reference time: 2.044 10 -5 s). The S0 and S1 computations have been done with 5 subcycles, see details in Section 3.1.1, whereas 50 subcycles have been used for the simulation S2. It should be noticed that the present space discretization is probably not sufficient. Indeed, the dimensionless ion Larmor radius is O(1/B), with B ≡ |B| = O(10 3 ) and the dimensionless surface of the computational domain equals S = 0.7816. Then, SB 2 = O(10 6 ) 31194 grid-points are required to go at the scale of the ion Larmor radius. This is why the following results should be considered as preliminary ones. Moreover, it would be required to go farther in time to provide relevant comparisons with experimental or numerical results of the literature. We aim at addressing such more realistic simulations in next future using a parallel version of the code.

In all simulations, the density and pressures evolve weakly, so that the visualizations look like those in Fig. 4. More attentively, the main variations are observed in the immediate vicinity of the targets, where the Bohm conditions apply. Visualizations of the ion velocity and of the current at t = 0.02 are provided for the Euler computation in Fig. 5. As expected, ions are outgoing at the targets and the current is essentially confined inside the separatrix. On animations of these fields, one observes high frequency oscillations of the poloidal components. Later on, these oscillations, that we do not associate to a turbulence phenomenon, die down and the system evolves with a much slower dynamics. Qualitatively, differences between the Euler and Braginskii computations are weak: when using the Braginskii closures, the vector fields are better aligned on the magnetic field, the oscillations are of weaker amplitude and more rapidly damped, especially for simulation S2.

Time variations are provided in Fig. 6, which shows the evolution at the point (r = 1.25, z = 0.056) of the z-component of the ion momentum. The stabilizing effect of the Braginskii closures is noticeable, and is clearly more important when increasing the values of the parallel transport coefficients, i.e. is more important for simulation S2 than for S1. From these curves, especially from the zoom in Fig. 6 (right), one can discern that the oscillations in the poloidal plane are characterized by two dominant frequencies: simplified theoretical analysis, it turns out to be useless to capture the very high frequency associated to the electron cyclotron motion. The maximum value of the magnetic field max B ≈ 2655 can be associated to the maximum frequency f ≈ 4226. Thus, contrarily to the space discretization, the time step δt = 5 10 -6 used in the computations was certainly small enough. For the true ratio m i /m e ≈ 1836, one can expect that it is enough to capture a frequency f = m i /m e f ic ≈ 42.8 max(f ic ). Visualizations of the potential at time t = 0.02 are provided for the simulations S0 and S1 in Fig. 7. One can again observe the stabilizing effect of the Braginskii closures, since the space variations are of weaker amplitude for the Braginskii computation. On animations one observes a rotation motion of the potential extrema.

Profiles of different quantities, along horizontal and (approximately) vertical mesh lines issued from the X-point, at time t = 0.1105 and for the simulation S1 are shown in Fig. 8. Some of them, namely the total density, total pressure and toroidal ion velocity profiles have not significantly evolved during the simulation. Clearly, the z-profiles do not show the pedestral shape typical of the High (H) mode regime and thus point out a Low (L) mode regime. Comparisons with experimental or numerical results of the literature, e.g. with the quasi-steady state regimes reported for the lower X-point configuration in [START_REF] Aydemir | Shear flows at the tokamak edge and their interaction with edge-localized modes[END_REF], could be of great interest, but it would be required to use a finer grid and to go farther in time to provide a relevant analysis.

Conclusion

A Fourier-SEM solver has been developed to solve the Braginskii system with the electrostatic assumption, and preliminary results have been presented for the JET tokamak configuration. A parallel version of the code, that makes use of a domain decomposition both in the poloidal plane and in the azimuthal direction, is close to be operational. This will allow, in next future, to carry out more relevant simulations and to provide comparisons with various other codes, including MHD codes. Some difficulties however remain: especially, we have used the heavy electron approximation and the parallel transport coefficients have been lowered; Also, one may think that more involved Bohm conditions are required (temperature dependent, more stable at the targets ends) to obtain fully reliable results.

A major extension to the present model would be to include the Maxwell equations. The electrostatic approximation, which assumes that the magnetic field is frozen, is only valid for the shortest time-scales. This may be satisfactory to exhibit the L-H transition but excludes the MHD underlying physics, which is of major importance to describe the instabilities (Edge-Localized Modes) that usually develop from the separatrix inside the SOL.

An extension of interest would be to include neutrals, since neutrals mainly provide source or sink terms to the physics of charged particles. This could be rather easy to do if a fluid approximation is also used for neutrals. The coupling with neutral transport code like EIRENE [START_REF] Reiter | The EIRENE and B2-EIRENE codes[END_REF] can also be considered, since such a coupling has already been realized for two fluid turbulence codes [START_REF] Bufferand | Numerical modeling for divertor design and in support to the WEST project[END_REF].

Another extension of the present work could include the coupling with multi-species or single species codes, e.g. to avoid the drift velocity approximation only where required. Concerning the modelling, the route to go from the kinetic equations to various fluid approximations, especially MHD and extended MHD, is well known. From the numerical implementation point of view it would be of course simpler to look for a coupling with codes using similar numerical methods, i.e. also based on finite element type approximations in the poloidal plane and Fourier expansion in toroidal direction, e.g. JOREK or NIMROD.

We aimed at (i) learning some information about the influence of the divergence-free constraint ∇ • (w i q i + w e q e ) = 0 on the behavior of q i and q e (see Section Appendix C.1), (ii) understanding how the solution of this toy model can be approximated by a projection method (see Section Appendix C.2). Let us assume that e i = -e e and introduce the ratio ε = m e /m i . The (dimensionless) values of the coefficients w s and α s can be expressed in terms of ε: w i = 1, w e = -1/ε, α i = 1/(1 + ε), α e = ε/(1 + ε), and the above system can be recasted using the total momentum q = q i + q e , the current j = w i q i + w e q e and Ψ = ρU :

       ∂ t q = j × B, ε∂ t j + ∇Ψ = q -(1 -ε)j × B, ∇ • j = 0. (C.1)
To make the calculation easy, we assume that the problem is solved in a smooth bounded 2D domain Ω embedded in the (x, z)-plane, that B = Be y is uniform and aligned with the y-axis and that the total momentum q 0 and the current j 0 at initial time satisfy q 0 , j 0 ∈ span{e x , e z }, where (e x , e y , e z ) stands for the canonical basis of R 3 . The system ensures that q, j ∈ span{e x , e z } for all time. Note also that, in this simplified (2D) geometry, we have, for any field u:

∇ • u × B = B ∇ × u, ∇ × u × B = -B ∇ • u, and u × B • τ = -B u • n, (C.2)
where the vectors colinear to e y are assimilated to scalar quantities, the vector n = (n x , n z ) stands for the (unit) outward normal vector at the boundary of the domain Ω, and τ = (n z , -n x ) is a tangent vector at the boundary of Ω.

Finally, the system is supplemented with the boundary condition j • n = 0 on the whole boundary. The first observation is that the boundary condition imposed on j implies a boundary condition for q. Indeed, from the first equation of (C.1) and the last equality of (C.2) we obtain:

∂ t q • τ = -B j • n = 0, on the boundary.
It means that the boundary condition q • τ = q 0 • τ (q 0 for the initial state) is implicitly enforced.

Due to these boundary conditions, the vectors j and q are entirely determined by their divergence and rotational (considered as a scalar). Let us denote d j = ∇ • j, r j = ∇ × j, d q = ∇ • q and r q = ∇ × q. By taking the divergence and the rotational of the two first equations of system (C.1) and bearing in mind that ∇ • j = 0, we find:

∂ t d q = B r j , (C.3) ∂ t r q = 0, (C.4) ε ∂ t r j = -B d q , (C.5) d j = 0. (C.6)
Thus, the system of PDEs (C.1) is replaced by a set of systems of ODEs (a system for each point of Ω). Each of them can be readily solved. Combining the time derivative of (C.5) and the equation (C.3), we obtain:

∂ tt r j = - B 2 ε r j .
This allows to find the exact expressions of r j and d q :

r j = r j 0 cos B √ ε t - 1 √ ε d q 0 sin B √ ε t , d q = √ εr j 0 sin B √ ε t + d q 0 cos B √ ε t .
Finally, from their divergence and rotational we can go back to j and q:

j = j 0 cos B √ ε t + 1 B √ ε q 0 × B -∇ Ψ0 sin B √ ε t , q = q 0 + j 0 √ ε B sin B √ ε t - 1 B 2 q 0 × B -∇ Ψ0 cos B √ ε t -1 × B,
where Ψ0 is defined as the solution of:

-∆ Ψ0 = -∇ • q 0 × B , ∇ Ψ0 • n = q 0 × B • n, on the boundary. (C.7)
Hence, both ion and electron quantities have a periodic behavior with the same angular frequency B/ √ ε. This unexpected analysis points out the influence of the electroneutrality, since if discarding the divergence free constraint we would obtain oscillations with angular frequency B for the ions, and B/ε for the electrons. For numerical applications, this means that it is not needed to resolve the electron cyclotron frequency to obtain accurate approximations of system (C.1).

Appendix C.2. Projection method

Following the same lines, we can analyze the projection method for the problem (C.1). First, let us note that the term (1 -ε)j × B in the right-hand-side of the second equation of (C.1) plays no role in the previous analysis (at least in the derivation of the exact expression for q and j). Indeed, in our simplified geometry this term is irrotational (since the current is divergence free) and, consequently, could be included in the definition of the potential Ψ. For this reason, in the sequel we introduce an additional coefficient µ and write the system as follows: (C.8)

       ∂ t q = j × B,
The system (C.1) considered previously corresponds to µ = ε. However, since the solutions j and q of the system (C.8) are the same for any values of µ, for the sake of simplicity we will consider the value µ = 1 at the end of this section.

The projection method starts from the exact initial states q 0 = q 0 and j 0 = j 0 . Then, assuming given the states q k and j k at the time t k = kδt, the algorithm to define q k+1 and j k+1 writes in two steps:

• in the first step, we compute predicted quantities q and ĵ discarding the divergence free constraint, that is we solve the system:    ∂ t q = j × B, ε∂ t j = q -(1 -µ)j × B, (C.9) on one time step. No boundary condition are needed at this step.

• then, starting from the predicted quantities q and ĵ, we compute the electric potential Ψ k+1 , the end-of-step current j k+1 and total momentum q k+1 by solving on one time step:

    
∂ t q = 0, ε∂ t j + ∇Ψ = 0, ∇ • j = 0.

(C.10)

The boundary condition j • n = 0 is enforced at this step.

As in the previous section, these two steps can be written using the variables d j = ∇ • j, r j = ∇ × j, d q = ∇ • q and r q = ∇ × q. The prediction step writes:

   ∂ t r q = -B d j , ∂ t d q = B r j , ε∂ t r j = -B d q -(1 -µ)d j , ε∂ t d j = B r q -(1 -µ)r j , (C.11) whereas the correction step writes ∂ t r q = 0, ∂ t d q = 0, ∂ t r j = 0, d j = 0 .

Thus, the correction step consists of setting d j to zero while preserving the quantities r q , d q and r j . It means that if we manage to eliminate d j from the system (C.11), we will obtain the system effectively solved for r q , d q , and r j when using the projection method. To this end, we take the time derivative of the equations on r q and r j , and replace ∂ t d j by its expression as a function of r q and r j . Finally, the projection method consists of solving the following system at each time step:

                   ∂ tt r q = - B 2 ε r q -(1 -µ)r j , ε∂ tt r j = - B 2 ε εr j -(1 -µ) r q -(1 -µ)r j , ∂ t d q = B r j , d j = 0 .
(C.12)

The initial conditions at time t k = kδt associated to this system are defined from the given state j k , q k . Since we have now second order equations in time, it is needed to define initial conditions for ∂ t r q and ∂ t r j . Because at the beginning of the time step d j k is zero, we obtain the following initial conditions from system (C.11): r q (t k ) = r q k , r j (t k ) = r j k , d q (t k ) = d q k , ∂ t r q (t k ) = 0,

∂ t r j (t k ) = - B ε d q k .
This system can be compared to the system (C.3)-(C.6) putting in light the splitting error. First, one observes that, contrarily to what happens for exact solution of (C.8), here the term (1 -µ)j × B plays a role. Even if the system (C.12) can be explicitly solved in the case where µ = ε, it is quite cumbersome to carry out the explicit form of the splitting errors in this case. However, we can easily go further if µ = 1, since in this case the system writes:

                   ∂ tt r q = - B 2 ε r q , ∂ tt r j = - B 2 ε r j , ∂ t d q = B r j , d j = 0 .
so that, the quantities r j , d q , d j may be accurately computed, whereas a splitting error is observed on r q : r q k = r q 0 cos B √ ε δt k .

Thus, for small time step δt, the splitting error at a given time T > 0 behaves like B 2 /(2ε)T δt. It is of order 1, as expected for the (non incremental) projection method, but more significantly the constant B 2 /(2ε) is in practice large. Because the magnetic field is strong, it turns out that the projection method for system (C.8) leads to an important splitting error even for the very small time step δt = 2πε/B, that is associated to the electron cyclotron period.

First, it iss

  needed to introduce the cyclotron angular frequencies ω cs = Be s /m s , where B = |B|, and the collision times τ s = m /(n s e 4 ). Let us also define u = (u • b)b, with b = B/B, for the projection of any vector field u onto the magnetic field lines, and u ⊥ = uu . Similar notations are used for the gradient operator, i.e. ∇ T = (∇T • b)b and ∇ ⊥ T = ∇T -∇ T for any scalar field T .

=

  -κ e ∇ T e -κ e ⊥ ∇ ⊥ T e + κ e ∧ b × ∇T e , κ e = 3.2 n e eT e τ e m e , κ e ⊥ = 1.47 κ e

-

  The parallel transport coefficients η and κ may show very high values. Assuming the following reference values (i.e. those used in Section 5): T = 100 eV , n = 10 19 m -3 , l = 2 m, m = m proton and e = e, and if we set u = T /m = 97861 ms -1 and t = l /u = 2.04 10 -5 s, then one obtains the results shown in Table

Figure 1 :

 1 Figure 1: Time-variation of the error for n el = 18 2 and different values of N (at left). Error with respect to N for different values of n el (at right).

2 :

 2 Max norm of the error for q (x) e at time t = 2π me/m i /B with respect to N for different values of n el , at left for n δt = 240, at right for n δt = 720.

5. 1 . 2 .

 12 Nondimensionalization Dimensionless values are hereafter systematically used. The values of reference quantities are chosen as follows: mass: m ≡ m i = 1.673 10 -27 kg (proton mass), electric charge: e ≡ -e e = 1.6022 10 -19 C, length: l = 2 m, temperature: T = 100 eV , particles density: n = 10 19 m -3 .

Figure 3 :

 3 Figure 3: Computational domain and mesh (at left); Visualization of the magnetic field (at right).

Figure 4 :

 4 Figure 4: Density ρ (at left) and temperature T i = Te (at right), at t = 0.

Figure 5 :

 5 Figure 5: Ion velocity (at left) and current density (at right) for the Euler computation, at t = 0.02. Colors are associated to the azimuthal component and arrows to the poloidal one.

f 1 ≈Figure 6 :

 16 Figure 6: Time-variation of the z-component of the ion momentum for the Euler and for the two Braginskii computations (at left) and zoom for t ∈ (0, 0.02) (at right).

Figure 7 :

 7 Figure 7: Potential for the Euler (at left) and Braginskii S1 (at right) computations at t = 0.02.

Figure 8 :

 8 Figure8: Profiles of the ion velocity components (top) and of the potential, the total density and pressure (bottom), along an horizontal line (at left) and an (approximately) vertical line (at right) issued from the X-point, at time t = 0.1105. The X is point is located at (x = 1.527, z = -0.639). For the x-profiles small quantities are multiplied by 10.

Appendix C. 1 .

 1 Exact solutions of model (C.1)

  ε∂ t j + ∇Ψ = q -(1 -µ)j × B, ∇ • j = 0.

Table 1 :

 1 Values of the product ωcτ for the ions and the electrons and for temperature-density pairs typical of the core and of the edge of the plasma

		10 8	5.8 10 5
	n (m -3 )	10 20	10 19
	ω ce τ e	3.39 10 7	1.2 10 5
	ω ci τ i	1.12 10 6	3.96 10 3

Table 2 :

 2 Dimensionless temperature, density, cyclotron frequency, collision time and parallel transport coefficients for the momentum and energy conservation equations.

		n	f c	τ	η	κ
	ions	0.298 0.5	291.8	2.21	0.316	1.28
	electrons 0.298 0.5 5.36 10 5 0.036 0.004	31.9
	ions	3.05 5.9	303	6.13 106.02 430.71
	electrons 3.05 5.9 5.58 10 5 0.101	1.33	10709

Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement number 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The authors wish also to thank T. Goudon and J.L. Guermond for fruitful discussions.

Appendix A. Proof of [START_REF] Aydemir | Shear flows at the tokamak edge and their interaction with edge-localized modes[END_REF] We aimed at deriving an explicit expression of A s = (I + a δt w s BM b ) -1 . Let us set γ s = a δt w s and consider u = A s v then we have:

Hence, from the last equality we can deduce that:

so that we have:

This prove (3).

Appendix B. Anisotropic diffusion test case

If the initial condition is constant on the magnetic surfaces associated to [START_REF] Buffereand | Near wall plasma simulation using penalization technique with the transport code SolEdge2D-Eirene[END_REF] and if the diffusion is purely parallel, i.e. such that κ ⊥ /κ = 0, then the exact solution, say u o , remains equal to the initial condition since ∇ • K∇u o = 0.

We assume now the following polar variation of the exact solution: u = u o f (θ), where f (θ) is a smooth 2π-periodic function. Since B θ = 1/r, one finds for the diffusion flux density ϕ = -K∇u:

where ϕ o = -K∇u o and f is the first derivative of f . Now one must look at the ∇ • ϕ term. Clearly, one has ∇ • ϕ o = 0 and since on the toroidal surfaces (where the diffusion takes place) B 2 r 2 is constant, one finds:

with f " for the second derivative of f . Let us assume that the exact solution is time dependent such that: u = u o f (θ)g(t). Then the source term is:

θ) = sin θ and g(t) = sin(πt/2) one obtains [START_REF] Blum | Reconstruction of the equilibrium of the plasma in a Tokamak and identification of the current density profile in real time[END_REF].

Appendix C. Study of a simplified model

We address in this section the following toy model:

∂ t q s = -w s α s ρ∇U + w s q s × B, s ∈ {i, e}, ∇ • (w i q i + w e q e ) = 0.