domain: math phys

text:
Holger Heumann

Francesca Rapetti

A finite element method with overlapping meshes for free-boundary axisymmetric plasma equilibria in realistic geometries

Existing finite element implementations for the computation of free-boundary axisymmetric plasma equilibria approximate the unknown poloidal flux function by standard lowest order continuous finite elements with discontinuous gradients. As a consequence, the location of critical points of the poloidal flux, that are of paramount importance in tokamak engineering, is constrained to nodes of the mesh leading to undesired jumps in transient problems. Moreover, recent numerical results for the self-consistent coupling of equilibrium with resistive diffusion and transport suggest the necessity of higher regularity when approximating the flux map. In this work we propose a mortar element method that employs two overlapping meshes. One mesh with Cartesian quadrilaterals covers the vacuum chamber domain accessible by the plasma and one mesh with triangles discretizes the region outside. The two meshes overlap in a narrow region. This approach gives the flexibility to achieve easily and at low cost higher order regularity for the approximation of the flux function in the domain covered by the plasma, while preserving accurate meshing of the geometric details outside this region. The continuity of the numerical solution in the region of overlap is weakly enforced by a mortar-like mapping Key-words: axisymmetric plasma equilibria in tokamaks; domain decomposition mortar method; overlapping meshes; linear and cubic

Weak Formulation

We choose a semi-circle Γ of radius ρ Γ surrounding the iron domain Ω Fe and the coil domains Ω ci . The truncated domain, we use for the computations, is denoted Ω ⊂ Ω ∞ , with boundary

]) be the space of vector fields in L 2 (Ω × [0, 2π]) 3 with curl in L 2 (Ω × [0, 2π]) 3 . We remark that: r -1 ψe φ ∈ H(curl, Ω × [0, 2π]) if and only if ψ ∈ H 1 (Ω).

Résumé : En éléments finis, les applications existantes pour le calcul d'équilibres de plasma à frontière libre en axisymétrique approchent le flux poloidale par une fonction continue à gradient discontinu, qui est localement dans chaque élément du maillage un polynôme de degré un. La position des points critiques du flux poloidale, qui est de grande importance pour les ingénieurs des tokamaks, est par consequence limitée aux seuls noeuds du maillage, et donc à l'origine de perturbations quand on modelise le passage du plasma du transitoire vers l'équilibre. De plus, des résultats numériques recents sur le couplage à l'équilibre entre le plasma et un modèle de diffusion/transport ont montré le besoin de plus de regularité dans l'approximation du flux. Dans ce travail on propose une approche par méthode d'éléments finis avec joints sur deux maillages qui se recouvrent. Un maillage structuré composé de rectangles qui couvre la chambre à vide du tokamak et un maillage non structuré de triangles pour la partie du domaine qui se trouve à l'exterieur de la chambre à vide. Les deux maillages sont superposés sur une bande etroite qui entoure la chambre à vide. Cet approche est assez flexible pour permettre d'avoir facilement et à bas coût numérique une approximation du flux poloidale plus regulière dans la chambre à vide tout en gardant une description precise des détails géometriques dans la partie externe. La continuité de la solution numérique dans la partie de recouvrement des maillages est imposée faiblement à l'aide de projections de type mortar.

Mots-clés : plasma à l'equilibre dans des tokamaks en axisymétrique; méthode de décomposition de domaine avec éléments joints (mortar); recouvrement de maillages; éléments finis linéaires et cubiques

Introduction

Computing plasma equilibria is maybe the most fundamental step in modeling for magnetic fusion applications. Main algorithmic approaches to the axisymmetric equilibrium problems were developed long time ago (we can refer to text books like [START_REF] Blum | Numerical simulation and optimal control in plasma physics[END_REF] or [START_REF] Jardin | Computational methods in plasma physics[END_REF] for details and references).

In general, one differentiates between the so-called free-boundary problem and the fixedboundary problem. The fixed-boundary problem is a semi-linear elliptic boundary value problem in the plasma domain with imposed Dirichlet data at the plasma/vacuum interface (which is assumed to be known, in this case). The most recent achievements [START_REF] Howell | Solving the Grad-Shafranov equation with spectral elements[END_REF][START_REF] Palha | A mimetic spectral element solver for the gradshafranov equation[END_REF] for this kind of equilibrium problem propose to use higher order/spectral methods to describe with high accuracy the physical field. Higher order/spectral methods [START_REF] Demkowicz | Computing with hp-adaptive finite elements[END_REF][START_REF] Schwab | p-and hp-finite element methods : theory and applications in solid and fluid mechanics[END_REF][START_REF] Karniadakis | Spectral/hp element methods for computational fluid dynamics[END_REF] are well-established approaches to approximate the solution of linear and non-linear elliptic boundary valued problems. They require isoparametric or transfinite mappings to guarantee the accuracy of the approximation as soon as the computational domain contains non-polygonal boundaries, as it occurs with the plasma modeling. The less canonical approach to higher order methods for the fixed-boundary equilibrium problem in [START_REF] Pataki | A fast, high-order solver for the grad-shafranov equation[END_REF] uses approximations of conformal mappings (see also [START_REF] Goedbloed | Conformal mapping methods in two-dimensional magnetohydrodynamics[END_REF]), nevertheless it does not cover the very important case of plasma boundaries with corners. To tackle the fixed-boundary equilibrium with prescribed curved boundary in the frame of high order methods is mainly an issue of proficiency in numerical methods for PDEs on curved domains. It is indeed important to stress that the boundary of the plasma is not known a priori. Assuming an arbitrarily detailed knowledge of the plasma boundary is not realistic for the actual physical application. The boundary is either deduced from measurements through experiments combined with reconstruction procedures or it is the output of the free-boundary equilibrium problem. In the latter case, which is the focus of this article, we have to solve a semi-linear elliptic problem in an unbounded domain where the region covered by the plasma is not known. Hence, it is of practical relevance to have numerical methods which are, to a very high degree, independent of the actual plasma boundary.

In this work we will focus on the free-boundary problem and propose an extension of the finite element (FE) method introduced in [START_REF] Blum | The self-consistent equilibrium and diffusion code SCED[END_REF] (see also [START_REF] Heumann | Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF]). A very important application for the free-boundary problem is the so-called self-consistent coupling of equilibrium with resistive diffusion and transport [START_REF] Grad | Classical diffusion in a tokamak[END_REF] that allows to simulate the evolution of the plasma equilibrium over very long time scales. Modeling in a numerically cheap and practical way the plasma/vacuum interface movement during the evolution of the plasma equilibrium is an issue of physical interest. Such simulations are essential for in silico studies of experiments in tokamaks but it turns out that the FE approach with piecewise polynomial, globally continuous approximations has two main drawbacks: 1.) The definition of the plasma boundary hinges on the critical points of the unknown flux. If the derivatives are not continuous, these points will not move in a continuous way during the evolution. 2.) The resistive diffusion and transport are described by one-dimensional equations containing metric coefficients that depend on the gradient of the solution of the equilibrium problem. These coefficients are not well-defined if the gradients are not continuous.

Hence, in this work we propose a very practical approach that allows to achieve easily continuous differentiability where it is beneficial but stays with the standard methods in the rest of the domain.

In general, it is very cumbersome to define FE spaces that allow for globally continuous differentiability. If one keeps along with triangular meshes, the reduced quintic element [START_REF] Bell | A refined triangular plate bending finite element[END_REF] is the lowest order FE with this property. Reduced quintic FEs were introduced for fusion applications firstly in [START_REF] Jardin | A triangular finite element with first-derivative continuity applied to fusion MHD applications[END_REF]. In the case of Cartesian meshes however, continuous differentiability can be easily achieved with lower order polynomials. By taking tensor products of one dimensional spline basis functions of polynomial order three we end up with a FE space of piecewise bicubic polynomials with continuous derivatives, also known as the Bogner-Fox-Schmit FE [START_REF] Bogner | The generation of interelement compatible stiffness and mass matrices by the use of interpolation formulas[END_REF][START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]. Inspired by [START_REF] Christophe | A mortar element approach on overlapping nonnested grids: application to eddy current non-destructive testing[END_REF] we present a FE method that employs two meshes, one of rectangles in the parts of the vacuum chamber domain that is accessible by the plasma, and one of triangles outside this region. This approach gives the flexibility to achieve easily and at low cost higher order regularity for the approximation of the flux function in the domain covered by the plasma, while preserving accurate meshing of the geometric details of passive structures and coils. As it is impossible to align the boundary of the mesh of rectangles with the interface of the exterior domain, we will allow for an overlap in a narrow region around the interface. The continuity of the numerical solution in the region of overlap is weakly enforced by relying on a mortar-like mapping.

The outline for the rest of the article is the following: The next two sections introduce the axisymmetric plasma equilibrium problem and present a weak formulation in a domain decomposition spirit with two distinct subdomains, one is part of the vacuum chamber accessible by the plasma and the other is the rest of the vacuum chamber together with the exterior domain. Section 4 is devoted to the presentation of the numerical method combining linear FEs for the exterior domain and the Bogner-Fox-Schmit FE in some parts of the vacuum chamber domain. Continuity is weakly enforced via a mortar-like mapping. Section 5 presents validation tests and applications from nuclear fusion science. We end with a short summary and outlook on perspectives in Section 6.

Free-Boundary Equilibrium of Toroidal Plasma

The essential equations for describing plasma equilibrium in a tokamak are force balance, the solenoidal condition and Ampère's law that read respectively

grad p = J × B, div B = 0, curl 1 µ B = J, (1)
where p is the plasma kinetic pressure, B is the magnetic induction, J is the current density and µ the magnetic permeability. In the quasi-static approximation these equations are augmented by Faraday's law in all other conducting structures, and by Ohm's laws in plasma, coils and passive structures.

For the considered setting, axial symmetry is a perfectly valid approximation. It is convenient to formulate (1) in a cylindrical coordinate system (r, ϕ, z) in order to consider only a section at ϕ = constant of the tokamak, generally referred to as poloidal section. We recall that we pass from (x, y, z) to (r, ϕ, z) by the transformation 1 x = r cos ϕ and y = r sin ϕ. Working in a poloidal section, the scalar field p does not depend on the angle ϕ, thus ∇p belongs to the poloidal (r, z)-plane. We introduce Ω ∞ = [0, ∞]×[-∞, ∞], the positive half plane, to denote the meridian plane that contains the tokamak centered at the origin. The geometry of the tokamak determines the various subdomains (see Fig. 1):

-Ω Fe ⊂ Ω ∞ denotes those parts of Ω ∞ made of iron; for an air-transformer tokamak Ω Fe = ∅; 1 We thus have a transformation in reference system, from (ex, ey, ez) to (er, eϕ, ez) ex = er cos ϕ -eϕ sin ϕ, ey = er sin ϕ + eϕ cos ϕ, and in the partial derivatives for any scalar field p, from (∂xp, ∂yp) to (∂rp, ∂ϕp), as follows ∂xp = ∂rp cos ϕ -∂ϕp sin ϕ r , ∂yp = ∂rp sin ϕ + 1 r ∂ϕp sin ϕ so that ∇p = er∂rp + eϕ 1 r ∂ϕp + ez∂zp. -Ω ci ⊂ Ω ∞ , 1 ≤ i ≤ N , denotes the intersection of the ith coil with the poloidal plane. We suppose that Ω ci has n i wire turns, total resistance R i and cross section area |Ω ci |;

Inria

-Ω L ⊂ Ω ∞ , denotes the domain bounded by the limiter, thus the domain accessible by the plasma;

-Ω p ⊂ Ω L , denotes the domain covered by the plasma.

The classical primal unknowns for toroidal plasma equilibria described by (1) are the poloidal magnetic flux ψ = ψ(r, z), the pressure p and the diamagnetic function f . The poloidal magnetic flux ψ := rA•e ϕ is the scaled toroidal component (ϕ-component) of the magnetic vector potential A, such that B = curl A, and e ϕ the unit vector for the ϕ coordinate. The diamagnetic function f = rB • e ϕ is the scaled toroidal component of the magnetic field B. It can be shown that both the pressure p and the diamagnetic function f are constant on ψ-isolines, i.e. p = p(ψ) and f = f (ψ). We refer to standard text books, e.g. [START_REF] Freidberg | Ideal Magnetohydrodynamics[END_REF], [START_REF] Blum | Numerical simulation and optimal control in plasma physics[END_REF], [START_REF] Wesson | Tokamaks. The International Series of Monographs in Physics[END_REF], [START_REF] Goedbloed | Principles of magnetohydrodynamics: with applications to laboratory and astrophysical plasmas[END_REF], [START_REF] Goedbloed | Advanced magnetohydrodynamics: with applications to laboratory and astrophysical plasmas[END_REF] and [START_REF] Jardin | Computational methods in plasma physics[END_REF] for the details and state in the following paragraphs only the final equations.

Force balance, the solenoidal condition and Ampère's law in (1) yield, in axisymmetric configuration, the following set of equations for the flux ψ(r, z):

-∇ • 1 µr ∇ψ =    rp (ψ) + 1 µ0r f f (ψ) in Ω p (ψ) ; I i /|Ω ci | in Ω ci ; 0 elsewhere in Ω ∞ , ψ(0, z) = 0 ; lim (r,z) →+∞ ψ(r, z) = 0 ; ( 2
)
where ∇ is the gradient in the half plane Ω ∞ , I i is the total current (in At, Ampère turns) in the ith coil and µ is a functional of ψ that reads

µ = µ Fe ( |∇ψ| 2 r 2 ) in Ω Fe µ 0 elsewhere , (3)
with µ 0 the constant magnetic permeability of vacuum and µ Fe the non-linear magnetic permeability of iron. Here again, we would like to stress that the plasma domain Ω p (ψ) is an unknown, which depends non-linearly on the magnetic flux ψ: the plasma domain Ω p (ψ) is a functional of the poloidal flux ψ. The different characteristic shapes of Ω p (ψ) are illustrated in Figure 1: the boundary of Ω p (ψ) either touches the boundary of Ω L (limiter configuration) or the boundary contains one or more saddle points of ψ (divertor configuration). The saddle points of ψ, denoted by (r X , z X )=(r X (ψ), z X (ψ)), are called X-points of ψ. The plasma domain Ω p (ψ) is the largest subdomain of Ω L bounded by a closed ψ-isoline in Ω L and containing the magnetic axis (r max , z max ). The magnetic axis is the point (r max , z max ) = (r max (ψ), z max (ψ)), where ψ has its global maximum in Ω L . For convenience, we introduce also the coordinates (r bdp , z bdp ) = (r bdp (ψ), z bdp (ψ)) of the point that determines the plasma boundary. Note that (r bdp , z bdp ) is either an X-point of ψ or the contact point with the limiter ∂Ω L . The equation ( 2) in the plasma domain, i.e.

-∂ r 1 µ 0 r ∂ r ψ -∂ z 1 µ 0 r ∂ z ψ = rp (ψ) + 1 µ 0 r f f (ψ), (4)
is the celebrated Grad-Shafranov-Schlüter equation [START_REF] Grad | Hydromagnetic equilibria and force-free fields[END_REF][START_REF] Shafranov | On magnetohydrodynamical equilibrium configurations[END_REF][START_REF] Lüst | Axialsymmetrische magnetohydrodynamische Gleichgewichtskonfigurationen[END_REF]. The domain of p and f f is the interval [ψ bdp , ψ max ] with the scalar values ψ max and ψ bdp being the flux values at the magnetic axis and at the boundary of the plasma:

ψ max (ψ) := ψ(r max (ψ), z max (ψ)) , ψ bdp (ψ) := ψ(r bdp (ψ), z bdp (ψ)) . (5)
The two functions p and f f and the currents I i in the coils are not determined by the model (2) and have to be supplied as data. Since the domain of p and f f depends on the poloidal flux itself, it is more practical to supply these profiles as functions of the normalized poloidal flux ψ N (r, z):

ψ N (r, z) = ψ(r, z) -ψ max (ψ) ψ bdp (ψ) -ψ max (ψ) . (6)
These two functions, subsequently termed S p and S f f , have, independently of ψ, a fixed domain [0, 1]. They are usually given as piecewise polynomial functions. Another frequent a priori model is

S p (ψ N ) = λ β r 0 (1 -ψ α N ) γ , S f f (ψ N ) = λ(1 -β)µ 0 r 0 (1 -ψ α N ) γ , (7)
with r 0 the major radius (in meters) of the vacuum chamber and α, β, γ ∈ R given parameters. We refer to [START_REF] Luxon | Magnetic analysis of non-circular cross-section tokamaks[END_REF] for a physical interpretation of these parameters. The parameter β is related to the poloidal beta [8, p. 15], whereas α and γ describe the peakage of the current profile, λ is a scaling parameter related to the total plasma current.

As we are going to present later a discretization scheme for the problem (2) that employs different approximation spaces on Ω L and its complement, we formulate the variational problem directly in a domain decomposition framework.

and

H 1 (Ω) = {u ∈ L 2 * (Ω), ∇u ∈ L 2 * (Ω) 2 } the Hilbert space endowed with the norm u 2 H 1 (Ω) = u 2 * ,Ω + |u| 2 H 1 (Ω) where |u| 2 H 1 (Ω) = ∂ r u 2 * ,Ω + ∂ z u 2 * ,Ω .
For ψ ∈ H 1 (Ω) the trace on Γ 0 vanishes in the following sense [START_REF] Haddar | Artificial boundary conditions for axisymmetric eddy current probe problems[END_REF] lim

r→0 + {r}×[-ρΓ,ρΓ]∩Ω r -2 ψ(r, z) 2 dz = 0 .
To formulate [START_REF] Albanese | On the solution of the magnetic flux equation in an infinite domain[END_REF] as variational problem in a domain decomposition framework, let us introduce the functional space

V = {(v, w) ∈ H 1 (Ω ex ) × (H 1 (Ω L ) ∩ C 0 (Ω L )), v |Γ0 = 0, v |γ = w |γ } ,
where Ω ex = Ω \ Ω L is the complement of Ω L in Ω. We require continuity in Ω L in order to have meaningful ψ max and ψ bdp that appear in the definition of Ω p and ψ N [8, Remark I.5, page 18]. It is not necessary to require differentiability to have a notion of maximum or minimum. Then, the weak formulation of (2) is:

Find (ψ ex , ψ L ) ∈ V s.t. a ex (ψ ex , v) + a L (ψ L , w) = ( I, v) ∀(v, w) ∈ V. (8)
In [START_REF] Blum | Numerical simulation and optimal control in plasma physics[END_REF], we have set

a ex (ψ, v) := Ω ex 1 µ(ψ)r ∇ψ • ∇v drdz + c(ψ, v), a L (ψ, w) := ΩL 1 µ 0 r ∇ψ • ∇w drdz -j p (ψ, w) ,
where

j p (ψ, ξ) := Ωp(ψ) rS p (ψ N ) + 1 µ 0 r S f f (ψ N ) ξ drdz , ( I, ξ) := N i=1 I i |Ω ci | Ωc i ξ drdz , (9)
and the bilinear form c(•, •) defined as

c(ψ, ξ) := 1 µ 0 Γ ψ(x)N (x)ξ(x)dS(x) + 1 2µ 0 Γ Γ (ψ(x) -ψ(y))M (x, y)(ξ(x) -ξ(y)) dS(x) dS(y) , (10)
accounts for the boundary conditions at infinity [START_REF] Albanese | On the solution of the magnetic flux equation in an infinite domain[END_REF], with x = (x r , x z ), y = (y r , y z ) and

M (x, y) = k x,y 2π(x r y r ) 3 2 2 -k 2 x,y 2 -2k 2 x,y E(k x,y ) -K(k x,y ) N (x) = 1 x r 1 δ + + 1 δ - - 1 ρ Γ and δ ± = x 2 r + (ρ Γ ± x z ) 2 .
Here, K and E are the complete elliptic integrals of first and second kind, respectively, and

k x,y = 4x r y r (x r + y r ) 2 + (x z -y z ) 2 .
We refer to [START_REF] Grandgirard | Modélisation de l'équilibre d'un plasma de tokamak[END_REF]Chapter 2.4] for the details of the derivation. The bilinear form c(•, •) follows basically from the so called uncoupling procedure in [START_REF] Gatica | The uncoupling of boundary integral and finite element methods for nonlinear boundary value problems[END_REF] for the usual coupling of boundary integral and FE methods.

In the case of vanishing plasma, S p = S f f = 0, the weak formulation ( 8) is the classical problem of non-linear magneto-statics; existence and uniqueness can be established under a monontonicity assumption for µ in the iron parts. The result follows directly from combining those for non-linear magneto-statics in [START_REF] Pechstein | Monotonicity-preserving interproximation of B-H-curves[END_REF] with the results for non-linear problems in unbounded domains [START_REF] Gatica | The uncoupling of boundary integral and finite element methods for nonlinear boundary value problems[END_REF]. If in addition we had Ω Fe = ∅ we would end up with an even simpler linear elliptic problem, for which existence and uniqueness are immediately available [START_REF] Haddar | Artificial boundary conditions for axisymmetric eddy current probe problems[END_REF][START_REF] Haddar | Axisymmetric eddy current inspection of highly conducting thin layers via asymptotic models[END_REF]. Rigorous existence and uniqueness assertion for the general case are still an open problem. See [START_REF] Temam | Remarks on a free boundary value problem arising in plasma physics[END_REF][START_REF] Berestycki | On a free boundary problem arising in plasma physics[END_REF][START_REF] Blum | Existence and control of plasma equilibrium in a tokamak[END_REF][START_REF] Nakamura | On an equilibrium of the plasma in a tokamak with a limiter[END_REF] for some theoretical work related to such results.

Finite Element Method with Overlapping Meshes

We wish to use, in the domain Ω L , FE approximations ψ h for the poloidal flux ψ that are not only continuous but have also continuous gradients ∇ψ h . This is particularly simple if we adopt a Cartesian mesh of rectangles in Ω L and use the bicubic Bogner-Fox-Schmit FE space [START_REF] Bogner | The generation of interelement compatible stiffness and mass matrices by the use of interpolation formulas[END_REF][START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]. The Bogner-Fox-Schmit rectangular FE is one of the simplest that provides continuous differentiability of the approximated solution. However, it can be applied only on a mesh of rectangles or paralellogramms, hence only to domains with boundaries that are parallel to two different straight fixed directions. As meshes with only rectangular elements are not very suitable for a realistic geometry description of the tokamak, we prefer to keep a mesh of triangles and linear Lagrangian FEs for the exterior domain Ω ex . It is clearly not possible to cover Ω L perfectly with a mesh of rectangles, we thus accept a certain overlap of the two meshes and use a mortarlike mapping to enforce continuity of traces at the interface. This modified version of the mortar element method for overlapping subdomains is related to the numerical zoom method in [START_REF] Hecht | Numerical zoom and the Schwarz algorithm[END_REF] and was introduced for non-distructive testing in [START_REF] Christophe | A mortar element approach on overlapping nonnested grids: application to eddy current non-destructive testing[END_REF], using a slightly different but equivalent formulation. The original mortar element method for overlapping subdomains and a convergence analysis of the method for a model problem can be found in [START_REF] Kuznetsov | Overlapping domain decomposition with non-matching grids[END_REF] and [START_REF] Cai | Overlapping nonmatching grid mortar element methods for elliptic problems[END_REF].

Preliminary Notation

We assume that the interface γ := ∂Ω L , between the domain Ω L , bounded by the limiter, and the exterior domain Ω ex , is polygonal and introduce a standard mesh τ ex of triangular elements that covers the domain Ω ex exterior to the domain Ω L . The boundary of Ω ex is ∂Ω ex = Γ ∪ Γ 0 ∪ γ. We assume that triangles T i of τ ex are shape regular and quasi-uniform. Next we introduce a (second, independent) mesh τ in of rectangular elements K j that covers a domain Ω in such that Ω L ⊂ Ω in . We assume that Ω in has a non-vanishing overlap with the domain Ω ex , that is,

Ω ov := Ω in ∩ Ω ex = ∅. Note that ∂Ω ov = γ ∪ γ with γ = ∂Ω in .
Both γ and γ are polygonal lines with nodes and edges from the meshes τ ex and τ in . We use N γ (resp. N γ ) to denote the set of all nodes of γ (resp. γ). Additionally we denote by N r γ (resp. N z γ ) all the nodes of γ that belong to an edge of γ that is parallel to the r-axis (resp. the z-axis). Moreover, we suppose that

Ω ci ∩ Ω in = ∅, 1 ≤ i ≤ N (together with Ω Fe ∩ Ω in = ∅ when Ω Fe = ∅).
We introduce two FE spaces over τ ex and τ in

V ex = {φ ∈ C 0 (Ω ex ) : φ |Γ0 = 0, φ |T is a linear polynomial ∀ T ∈ τ ex } V in = {ϕ ∈ C 1 (Ω in ), ∂ 2 rz ϕ ∈ C 0 (Ω in ) : ϕ |K is a bicubic polynomial ∀ K ∈ τ in }
and denote by V ex ∂ (resp. V in ∂ ) the trace space of V ex (resp. V in ) on the closed polygonal line Inria The green circles are the nodes in N z γ , the nodes of γ that belong to an edge of γ that is parallel to the r-axis, and the red crosses are the nodes in N r γ , the nodes of γ that belong to an edge of γ that is parallel to the z-axis. Both red crosses and green circles belong to N γ .

γ (resp. γ), namely, V ex ∂ = {φ |γ , φ ∈ V ex } (resp. V in ∂ = {ϕ |γ , ϕ ∈ V in }).
Note that V ex is the standard linear Lagrangian FE space and V in is known as the Bogner-Fox-Schmit FE space. Then the trace space V ex ∂ ⊂ C 0 (γ) is the span of affine linear functions defined on the mesh over γ that is induced by τ ex . The degrees of freedom of V ex can be chosen to be the nodal values at the nodes of τ ex that are not in Γ 0 , and the degrees of freedom of the trace space V ex ∂ can be chosen to be the nodal values on nodes in N γ , i.e. the nodes of τ ex that are on γ (see Figure 2, right). For V in the standard choice of degrees of freedom involves the nodal values of the function, its first order partial derivatives and its second order mixed derivatives at the nodes of τ in . Then, as degrees of freedom for the Dirichlet trace space V in ∂ ⊂ C 0 (γ) one can simply choose the nodal values at nodes of γ together with the nodal values of the partial derivative in the r-direction (resp. z-direction) if the node belongs to an edge of γ parallel to the r-direction (resp. z-direction).

The two mortar-like mappings

Since v ex |γ ∈ V in ∂ (resp. v in |γ ∈ V ex ∂ ) for v ex ∈ V ex (resp. v in ∈ V in )
to impose the transmission condition at γ and γ at the discrete level we rely on two mortar-like mappings, called nodal interpolations. More precisely, it is possible to define two operators π in :

V ex → V in ∂ and π ex : V in → V ex ∂ based on the particular choice of degrees of freedom of V ex ∂ and V in ∂ : a) For ψ in ∈ V in we define π ex ψ in ∈ V ex ∂ such that (π ex ψ in )(r i , z i ) := ψ in (r i , z i ) , ∀(r i , z i ) ∈ N γ . (11)
b) For

ψ ex ∈ V ex we define π in ψ ex ∈ V in ∂ such that (π in ψ ex )(r i , z i ) := ψ ex (r i , z i ) , ∀(r i , z i ) ∈ N γ , ∂ r (π in ψ ex )(r i , z i ) := ∂ r ψ ex (r i , z i ) , ∀(r i , z i ) ∈ N r γ , ∂ z (π in ψ ex )(r i , z i ) := ∂ z ψ ex (r i , z i ) , ∀(r i , z i ) ∈ N z γ . ( 12
)
We remark that π in is not well defined when the nodes (r i , z i ) ∈ N r γ or (r i , z i ) ∈ N z γ happen to lie on an edge of τ ex , since the gradients of functions in V ex are not single valued on edges of τ ex . In the rare cases where we need to manage this multivalued situation, we simply choose one of the possible values. The difference tends to zeros for decreasing mesh size, and hence the impact of this choice is not very crucial.

The mortar-like Galerkin formulation

By the definition of V ex and V in there are FE spaces V ex

• and V in • such that

V ex = V ex • ⊕ EV ex ∂ and V in = V in • ⊕ EV in ∂ ,
where E denotes the trivial extension operators. The elements of V ex • and V in • have vanishing Dirichlet trace on γ and γ, respectively.

We are now able to formulate the discrete variational problem:

Find (ψ ex , ψ in ) ∈ V ex × V in such that a ex (ψ ex , v) + a in (ψ in , w) = ( I, v) ∀(v, w) ∈ V ex • × V in • , ψ ex -π ex ψ in = 0 on γ , ψ in -π in ψ ex = 0 on γ , (13)
Inria where a ex (ψ, v) :=

Ω ex 1 µ(ψ)r ∇ψ • ∇v drdz + c(ψ, v) , a in (ψ, w) := Ω in 1 µ 0 r ∇ψ • ∇w drdz -j p (ψ, w) . (14)
The mappings ( I, v), j p (ψ, w) and c(ψ, v) are the same as in ( 9) and ( 10). An important difference to the method presented in [START_REF] Kuznetsov | Overlapping domain decomposition with non-matching grids[END_REF] is that we do not introduce weighting coefficients to compensate for the twofold integration over the overlapping domain and hence the consistency error can not be estimated as easily as in [13, Section 4.2]. On the other hand it is the absence of weighting coefficients that avoids the need of computing the polygonal intersections of triangles and rectangles for the assembling of a in (ψ, w) and a ex (ψ, v). With a non-zero overlap, we expect that the consistency error can nevertheless be controlled. At least the numerical experiments in [15, Section 5.] show convergence of optimal order for a eddy current problem in non-destructive testing using this kind of mortar approach.

Remark 1 Most of the theoretical work on mortar element methods [START_REF] Wohlmuth | Discretization methods and iterative solvers based on domain decomposition[END_REF][START_REF] Bernardi | A new nonconforming approach to domain decomposition: the mortar element method[END_REF] assumes the operators π ex and π in to be L 2 -projections. It is standard to define the operator

π ex : H 1 (Ω in ) → V ex ∂ and π in : H 1 (Ω ex ) → V in ∂ as follows: γ π ex ψ ξ dS = γ ψ ξ dS ∀ξ ∈ V ex ∂ ,
and

γ π in ψ ξ dS = γ ψ ξ dS ∀ξ ∈ V in ∂ .
The choice of π ex and π in as L 2 -orthogonal projections guarantees stability in the H s -norm, 0 ≤ s ≤ 1. Unfortunately, in practical implementations, this choice leads to undesired difficulties due to integration of products of FE functions defined on different meshes. This approach would require to find intersections of edges of one mesh with the elements from the other mesh. The nodal interpolations in Section 4.2 require only to find the element of the first mesh where a node of the second mesh is located.

4.4 Newton's method for [START_REF] Cai | Overlapping nonmatching grid mortar element methods for elliptic problems[END_REF] Note that a ex (•, •) (resp. a in (•, •)) in ( 13) is linear in the second argument but not in the first, due to the non-linear dependence of the physical coefficient µ (resp. of j p ) on the solution ψ. Hence, we will use Newton-type methods to find solutions of the discrete problem associated with [START_REF] Cai | Overlapping nonmatching grid mortar element methods for elliptic problems[END_REF]. This amounts to iterate the following update rule for (ψ k+1 ex , ψ k+1 in ):

d ψ a ex (ψ k ex , v)(ψ k+1 ex -ψ k ex ) + d ψ a in (ψ k in , w)(ψ k+1 in -ψ k in )+ = ( I, v) -a m ((ψ k ex , ψ k in ), (v, w)) ∀(v, w) ∈ V ex • × V in • , ψ k+1 ex -π ex ψ k+1 in = 0 on γ ψ k+1 in -π in ψ k+1 ex = 0 on γ
For the non-linear mapping a ex (•, •), taking into account (3), we have

d ψ a ex (ψ, ξ)( ψ) =c( ψ, ξ) + Ω ex 1 µ(ψ)r ∇ ψ • ∇ξ drdz -2 ΩFe µ Fe (|∇ψ| 2 r -2 ) µ 2 Fe (|∇ψ| 2 r -2 )r 3 (∇ ψ • ∇ψ) (∇ψ • ∇ξ) drdz , RR n°8916
and for the non-linear mapping a in (•, •), we have

d ψ a in (ψ, ξ)( ψ) = Ω L 1 µ 0 r ∇ ψ • ∇ξ drdz + j p (ψ; ξ, ψ) , (15)
where j p (ψ; ξ, ψ) is the approximation of the derivative

d ψ j p (ψ, ξ)( ψ) = Ωp(ψ) rS p (ψ N (ψ)) + 1 µ 0 r S f f (ψ N (ψ)) d ψ ψ N (ψ)( ψ) ξ drdz, - ∂Ωp(ψ) rS p (1) + 1 µ 0 r S f f (1) |∇ψ| -1 ( ψ -ψ(r bd (ψ), z bd (ψ)))ξ dΓ (16)
given in [START_REF] Blum | Numerical simulation and optimal control in plasma physics[END_REF]Lemma I.4], where

(d ψ ψ N (ψ)( ψ))(r, z) = ψ(r, z) -ψ N (ψ) ψ(r bdp (ψ), z bdp (ψ)) -(1 -ψ N (ψ)) ψ(r max (ψ), z max (ψ)) ψ bdp (ψ) -ψ max (ψ) .
The derivation involves shape calculus [START_REF] Murat | Sur le contrôle par un domaine géométrique[END_REF][START_REF] Delfour | Shapes and geometries[END_REF] and the non-trival derivatives:

d ψ ψ max (ψ)( ψ) = ψ(r max (ψ), z max (ψ)) and d ψ ψ bdp (ψ)( ψ) = ψ(r bdp (ψ), z bdp (ψ)) .
There are two different approaches to introduce approximations j p (ψ; ξ, ψ) of the derivative of the non-linear mapping j p (ψ, ξ)( ψ). The first replaces the integration in the analytic expression ( 16) of the derivative with standard quadrature rules. The second introduces numerical quadrature to approximate the integrals in the analytical expression (9) of the non-linear mapping and uses the analytical derivative of this approximation. While in many cases, including the derivatives of a ex (ψ, ξ), the two approaches yield the same approximation, this is not the case for j p (ψ, ξ) and we refer to [33, Section 3] for a detailed discussion on this topic. As one generally establishes convergence of numerical solutions of a discretization of a non-linear problem towards the exact solution, it is more natural to follow the second approach and to calculate analytically derivatives of discretized non-linear mappings. We use Gauss-Legendre quadrature of order 5 for integrals on rectangles. Moreover, in the approximation of

j p (ψ, ξ) = Ωp(ψ) rS p (ψ N ) + 1 µ 0 r S f f (ψ N ) ξ drdz
we do not compute exactly the intersection K ∩ Ω p (ψ) of an element K with the plasma domain Ω p (ψ), but extend S p (ψ N ) and S f f (ψ N ) by zero when ψ N > 1. Assuming that S p (•) and S f f (•) are smooth mappings from R + to R + we can use indeed standard quadrature rules on the rectangular elements K.

Let us underline that the second term on the righthand side of ( 16) vanishes whenever S p (1) = S f f (1) = 0. Moreover, the second and third term of d ψ ψ max (ψ)( ψ) make d ψ j p (ψ, ξ)( ψ) non-local in the sense that it does not vanish for ξ and ψ with disjoint supports, whenever (r max (ψ), z max (ψ)) or (r bdp (ψ), z bdp (ψ)), respectively, is in the support of ψ.

The computation of (r max (ψ in ), z max (ψ in )) and (r bdp (ψ in ), z bdp (ψ in )) for a piece wise bicubic polynomial ψ in ∈ V in is much more involved than in the case of linear Lagrangian FEs. Indeed, with piece wise linear Lagrangian FEs, the critical points are located at nodes of the mesh, whereas with high-order FEs, they occupy a position that does not coincide necessarily with a Inria node but can be at the interior of an element. For the determination of (r max (ψ in ), z max (ψ in )) we look first for the maximum of ψ in at a finite number of evenly distributed points. Such an initial guess is then refined by looking for a critical point in the neighborhood which can be accomplished with standard algorithms for constrained optimization problems, where the objective is the minimization of |∇ψ in | 2 /2 and the boundaries of the rectangular element set the constraints. In the current implementation we rely on the interior point method [START_REF] Nocedal | Numerical optimization[END_REF]Chapter 19], being the default algorithm in the function fmincon of MATLAB 2015a. For the determination of saddle points of ψ in ∈ V in we interpolate ψ in first onto the lower dimensional bilinear FE space and compute the saddle points for this representation. Then we refine the location by solving again a constrained minimization problem.

Algebraic Newton iterations

We recall the direct decomposition of V ex and V in

V ex = V ex • ⊕ EV ex ∂ and V in = V in • ⊕ EV in ∂ .
Then, if u ex and u in represent the vector of the values of degrees of freedom of ψ ex ∈ V ex and ψ in ∈ V in we have the decomposition

u ex = (u ex • , u ex ∂ ) and u in = (u in • , u in ∂ ), where u ex • (resp. u in • ) and u ex ∂ (resp. u in ∂ ) are the degrees of freedom in V ex • (resp. V in • ) and V ex ∂ (resp. V in ∂ ).
The matrix form of conditions ( 11) and ( 12) thus read, respectively:

u in ∂ = P in ∂,• u ex • + P in ∂,∂ u ex ∂ , u ex ∂ = P ex ∂,• u in • + P ex ∂,∂ u in ∂ . (17)
More precisely, if {b in •,i (r, z)} i and {b in ∂,j (r, z)} j are the sets of basis functions corresponding to the degrees of freedom in u in

• and u in ∂ we have

(P ex ∂,• ) kj = b in •,j (r k , z k ) and (P ex ∂,∂ ) kj = b in ∂,j (r k , z k ) ∀(r k , z k ) ∈ N γ .
Likewise, if {b ex •,i (r, z)} i and {b ex ∂,j (r, z)} j are the sets of basis functions corresponding to the degrees of freedom in u ex

• and u ex ∂ we have

(P in ∂,• ) kj = b ex •,j (r k , z k ) and (P in ∂,∂ ) kj = b ex ∂,j (r k , z k ) ∀(r k , z k ) ∈ N γ , (P in ∂,• ) kj = ∂ r b ex •,j (r k , z k ) and (P in ∂,∂ ) kj = ∂ r b ex ∂,j (r k , z k ) ∀(r k , z k ) ∈ N r γ , (P in ∂,• ) kj = ∂ z b ex •,j (r k , z k ) and (P in ∂,∂ ) kj = ∂ z b ex ∂,j (r k , z k ) ∀(r k , z k ) ∈ N z γ .
Hence, the assembling of ( 17) amounts mainly to find the element in which each node (r k , z k ) is located. See Figure 3 for details. The weak formulation (13) yields the following non-linear algebraic system:

A ex • (u ex • , u ex ∂ ) = F ex • A in • (u in • , u in ∂ ) = 0 u ex ∂ -P ex ∂,• u in • -P ex ∂,∂ u in ∂ = 0 u in ∂ -P in ∂,• u ex • -P in ∂,∂ u ex ∂ = 0 (18)
where

A ex • (u ex • , u ex ∂ ) is the discretization of the non-linear mapping a ex (•, •) and A in • (u in • , u in ∂ )
is the discretization of the non-linear mapping a in (•, •). Newton's method is used to solve iteratively the non-linear algebraic system [START_REF] Cuvelier | An efficient way to assemble finite element matrices in vector languages[END_REF]. At each iteration we solve a linear system of the following form:

    A ex •,• A ex •,∂ 0 0 0 0 A in •,• A in •,∂ 0 I ex ∂,∂ -P ex ∂,• -P ex ∂,∂ -P in ∂,• -P in ∂,∂ 0 I in ∂,∂         ∆u ex • ∆u ex ∂ ∆u in • ∆u in ∂     =     F ex • -A ex • (u ex • , u ex ∂ ) -A in • (u in • , u in ∂ ) P ex ∂,• u in • + P ex ∂,∂ u in ∂ -u ex ∂ P in ∂,• u ex • + P in ∂,∂ u ex ∂ -u in ∂     . The different blocks A ex •,• , A ex •,∂ (resp. A in •,• , A in •,∂
) correspond to the derivatives of the discretized non-linear mappings a ex (•, •) (resp. a in (•, •)) that we introduced in (14).

Remark 2 In the case of mortar-like mappings defined by L 2 -projection we have in (17)

P in ∂,• = (M in ) -1 Q in ∂,• , P in ∂,∂ = (M in ) -1 Q in ∂,∂ , P ex ∂,• = (M ex ) -1 Q ex ∂,• , P ex ∂,∂ = (M ex ) -1 Q ex ∂,∂ , with (M in ) ij = γ b in ∂,i b in ∂,j drdz , (M ex ) ij = γ b ex ∂,i b ex ∂,j drdz , (Q in ∂,• ) ij = γ b in ∂,i b ex •,j drdz , (Q in ∂,∂ ) ij = γ b in ∂,i b ex ∂,j dr dz , (Q ex ∂,• ) ij = γ b ex ∂,i b in •,j drdz , (Q ex ∂,∂ ) ij = γ b ex ∂,i b in ∂,j dr dz .
In contrast to the mortar-like mapping via nodal interpolation, we need to find not only the location of nodes of one mesh in elements of the other mesh but also the location of the edges.

Inria

Let us explain this for the setting in Figure 3. Basis functions and degrees of freedoms are associated to a particular node of the mesh and the basis functions vanish on all elements that do not share this particular node. If b in •,4 denotes a basis function of V in associated to the node 4 of τ in in Figure 3 and b ex ∂,3 denotes a basis function of V ex associated to the node 3 of τ ex we have that

(Q ex ∂,• ) 34 = γ83 b ex ∂,3 b in •,4 drdz + γ34 b ex ∂,3 b in •,4 drdz ,
where γ 83 and γ 34 are the edges of γ between the nodes 8, 3 and 4. It is not possible to approximate these integrals directly with quadrature as the function b in •,4 is not smooth on the edges γ 83 and γ 34 . Hence, we need to find the intersection of these edges with the elements of the mesh τ in and decompose the edges into lines where both basis functions b in •,4 and b ex ∂,3 are smooth. The edge γ 34 for example, is split into two segments e 1 and e 2 that are contained in different elements of τ in .

Similarly for the computation of

Q in ∂,∂ and Q in ∂,
• the edges of γ need to be split into segments contained in elements of the triangular mesh τ ex . The edge γ23 of γ between the nodes 2 and 3 for example, is split into three parts b 1 , b 2 and b 3 .

Geometric Coefficients

A very important output of equilibrium calculations are the so called geometric coefficients. The Grad/Hogan approach [START_REF] Grad | Classical diffusion in a tokamak[END_REF] to the simulation of evolution of plasma in a tokamak on very long timescales, asserts that the fluid model quantities such as densities or temperatures are constant on the level lines of the poloidal flux. Hence, transforming the corresponding conservation laws into a curvilinear coordinate system with one coordinate line aligned with the level lines of the poloidal flux ψ, we end up with a system of one dimensional equations, with metric coefficients due to the non-linear coordinate transformation. As the coordinate transformation depends on ψ, also the metric coefficients depend on ψ. In the tokamak literature these coefficients are referred to as geometric coefficients (see [START_REF] Blum | Plasma equilibrium evolution at the resistive diffusion timescale[END_REF] for a concise introduction to this topic focusing on numerical methods). The Grad/Hogan approach is implemented in many productive simulation tools [START_REF] Khayrutdinov | Studies of plasma equilibrium and transport in a tokamak fusion device with the inverse-variable technique[END_REF][START_REF] Crotinger | Corsica; a comprehensive simulation of toroidal magnetic-fusion devices[END_REF][START_REF] Zaitsev | The numerical solution of the self-consistent evolution of plasma equilibria[END_REF][START_REF] Honda | Simulation technique of free-boundary equilibrium evolution in plasma ramp-up phase[END_REF] that are used to study the evolution of plasma in tokamaks.

More precisely, the geometric coefficients are non-linear functionals of the following form

g c,ψN (y) = {(r,z)∈Ωp(ψ), ψN(r,z)=y} c(r, z)r |∇ψ N (r, z)| ds , (19)
that, for given smooth scalar functions c : Ω → R, are integrals along the level line {(r, z) ∈ Ω, ψ N (r, z) = y}. For the numerical approximation of g c,ψN (y) we need to find the elements that intersect with the level line {(r, z) ∈ Ω, ψ in N (r, z) = y}, where ψ in N is the normalization based on ψ in ∈ V in , the solution of the Galerkin formulation [START_REF] Cai | Overlapping nonmatching grid mortar element methods for elliptic problems[END_REF]. Since ψ in is piecewise bicubic, we have only an implicit representation of the level line. It is not possible, not even in each element K, to have a closed form expression for a parametrization s : [0, t] → Ω such that

{(r,z)∈Ωp(ψ), ψ in N (r,z)=y}∩K c(r, z)r |∇ψ in N (r, z)| ds = 1 0 |ṡ(t)| c(s(t))s r (t) |∇ψ in N (s(t))| dt .
To exploit nevertheless the high order polynomial representation of ψ in we use the Simpson

RR n°8916

quadrature rule .

{(r,z)∈Ωp(ψ), ψ in N (r,z)=y}∩K c(r, z)r |∇ψ in N (r, z)| ds = 1 0 |ṡ(t)| c(s(t))s r (t) |∇ψ in N (s(t))| dt ≈ 1 6 |ṡ(0)| c(s(0))s r (0) |∇ψ in N (s(0))| + 4|ṡ(0.
with s(0) and s(1) the two intersection points of the level line with the boundary ∂K of the element. To determine the intermediate point s(0.5) and the tangent vectors ṡ(0), ṡ(0.5) and ṡ(1) we follow the procedure outlined in [20, p. 199]: We write any point s(t) of the level line {(r, z)∈ Ω p (ψ), ψ in N (r, z) = y} ∩ K as the intersection of two level lines:

ψ in N (s(t)) = 0 , (1 -t)A(s(t)) + tB(s(t)) = 0 , ( 20
)
where the second level line is implicitly defined as convex combination of two lines, different from the level line of ψ, that intersect the level line of ψ in the two end points s(0) and s(1). In our calculations we decided for affine functions

A(s) = (s -s(0)) • (s(1) -s(0)) , B(s) = (s -s(1)) • (s(1) -s(0)) .
Hence, the intermediate point s(0.5) verifies the non-linear problem ψ in N (s(0.5)) = 0 , 0.5(s(0.5) -s(0)) • (s(1) -s(0)) + 0.5(s(0.5) -s(1)) • (s(1) -s(0)) = 0 .

Differentiating [START_REF] Demkowicz | Computing with hp-adaptive finite elements[END_REF] with respect to t we see that the tangent vectors ṡ(t) verify

∇ψ in N (s(t)) • ṡ(t) = 0 , (1 -t)∇A(s(t)) • ṡ(t) + t∇B(s(t)) • ṡ(t) = A(s(t)) -B(s(t)) ,
which is a linear problem once we know s(t).

Numerical results

In this section we present some numerical results that highlight the features of the proposed method. All the implementations and experiments were done with FEEQS.M3 . FEEQS.M is a MATLAB implementation of the methods for axisymmetric free boundary plasma equilibria that are described in [START_REF] Heumann | Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF]. The code utilizes in large parts vectorization, and therefore, the running time is comparable to C/C++ implementations (see [START_REF] Koko | Vectorized matlab codes for linear two-dimensional elasticity[END_REF][START_REF] Chen | Programming of finite element methods in Matlab[END_REF] and [START_REF] Cuvelier | An efficient way to assemble finite element matrices in vector languages[END_REF] for a review and earlier references). FEEQS.M is publicly available and a forthcoming release will contain the here introduced overlapping mesh methods for plasma equilibrium calculations.

We start with examples that show qualitatively the feasibility of the mortar-like FE method (MEM) for overlapping subdomains introduced in Section 4. Next we study numerically the dependence of the location of critical points, such as saddle points and maxima, on the values of the currents in the poloidal field coils. This is a very important application for scenario design in tokamaks. We finish this section by presenting results for the geometric coefficients, the level line integrals introduced in Section 4.6, which are necessary for simulations of transient plasmas in tokamaks.

All subsequent application examples, if not stated differently, are based on the WEST tokamak (see Figure 1 for a sketch of the different subdomains Ω L , Ω ci and Ω Fe4 . The imposed currents and the numbering of the coils can be inferred from Figure 4. For parameters in the current profile [START_REF] Blommaert | Towards automated magnetic divertor design for optimal heat exhaust[END_REF] we choose α = 0.9, β = 1.5, γ = 0.9, λ = 1806600 and R 0 = 2.4m. All the computations where performed on a MacBook Pro with a 2,8 GHz Intel Core i7 processor and 16 GB RAM, using MATLAB 2015a.

Feasibility

The initial guess of the plasma domain Ω p (ψ) for given currents in the poloidal field coils plays a crucial role in free-boundary equilibrium problems. It is common to find such initial guesses by trial and error. More sophisticated approaches are based on the formulation of inverse problems or optimal control problems, where a desired shape and position of the plasma domain is the objective and the precise values of the currents is unknown [START_REF] Blum | Numerical simulation and optimal control in plasma physics[END_REF][START_REF] Heumann | Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF]. The virtual casing principle by Shafranov and Zakharov [START_REF] Shafranov | Use of the virtual-casing principle in calculating the containing magnetic field in toroidal plasma systems[END_REF] is another fairly intuitive way of determining coil currents for a desired plasma boundary. In the present case we do not focus on this technical issue, but assume we have a good initial guess for the poloidal flux ψ, e.g., from a non-mortar formulation of the free-boundary equilibrium problem that is based on linear Lagrangian FEs. In Figure 5 we show the contour plots of the solution of the MEM for an increasing number of elements (see the table in Figure 6) of the interior rectangular mesh τ in . Not very surprisingly the solutions do not differ much and are close to the one without MEM. The visualization in Figure 6 emphasizes this observation in focusing on the plasma boundary and the data in the table of the same figure give more quantitative evidence. We need less than 10 Newton iterations to reduce the relative residual of the non-linear discrete system to values below 10 -12 .

The MEM has higher complexity, than the standard P1 method. Nevertheless, most of the additional effort for assembling the algebraic systems is negligible and can be done before the Newton loop starts. The largest additional computational effort, is most likely the inversion of the matrix in the Newton iterations, but as we work here with two dimensional problems we 13) and Bogner-Fox-Schmit (BFS) elements for an increasing number of elements in the interior rectangular mesh τ in (1st, 2nd and 3rd from the left) and a solution without MEM but with linear Lagrangian finite elements. can rely on the very efficient direct solvers, that work very well also for non-structured matrices (see the table in Figure 6 for some timings). The implementations here are not yet run time optimized, so for the moment we are not able to make general assertions about the performances. The experiments show that with minimal additional computational effort it is possible to do freeboundary equilibrium calculations with continuous derivatives.

As we combine linear Lagrangian FEs with Bogner-Fox-Schmit FEs, we can in general only expect that the error reduction in each mesh refinement step is comparable to the error reduction of a solution without MEM relying on linear Lagrangian FEs everywhere. Clearly, a genuinely high order MEM combines high order FEs on the exterior with high order FEs on the interior. We refer to [START_REF] Christophe | A mortar element approach on overlapping nonnested grids: application to eddy current non-destructive testing[END_REF]Section 5] for convergence studies on a linear problem using a similar MEM and lowest order spaces both in the exterior and the interior domains. Numerical analysis and extensive numerical experiments for a related method can be also found in [START_REF] Flemisch | A domain decomposition method on nested domains and nonmatching grids[END_REF].

To validate the convergence for the setting of this work, we consider the linear problem

-∇ • (∇u) = f in Ω , u = u 0 in ∂Ω , (21)
on a rectangular domain Ω = [-1, 1] 2 and define Ω L as the polygon with vertices (-0.125, 0.5), (0.375, 0.25), (0.375, -0.375), (0, -0.5), (-0.375, -0.375), and (-0.5, 0.25). The meshes τ in and τ ex for the interior domain Ω in and exterior domain Ω ex will be a Cartesian mesh and a triangular mesh. If h ex (resp. h in ) is the maximal diameter of elements in τ ex (resp. τ in ), and p ex (resp. p in ) the local polynomial degree of the FE spaces V ex (resp. V in ), one has optimal convergence if, for a smooth solution, the approximation error in the H 1 (Ω ex ) and H 1 (Ω in )-norms behaves as O(h p-1 ), with h = max(h ex , h in ) and p = min(p ex , p in ). This reasoning is confirmed by the numerical experiments (see Figure 7), where we took in [START_REF] Flemisch | A domain decomposition method on nested domains and nonmatching grids[END_REF] the data f and u 0 such that u(r, z) = cos(πr) sin(πz) is the solution. We do not observe any quantitative difference between the MEM using either the L 2 -projection or the nodal interpolation in the coupling condition. The use of a bicubic polynomial in the interior domain Ω in allows to achieve with the MEM a given error level with elements in Ω in larger than those used with linear Lagrangian FEs everywhere. Due to the behavior of this particular solution, we cannot expect further profit from the Bogner-Fox-Schmit FE.

Next, taking a clue from the theory of hp-FE methods we consider the data f (r, z) and ψ 0 such that ψ(r, z) = cos(πr) 4 cos(πz) 4 , for (r, z) ∈ [-0.5, 0.5] 2 , and ψ(r, z) = 0, for (r, z) on Ω \ [-0.5, 0.5] 2 , is the solution of [START_REF] Flemisch | A domain decomposition method on nested domains and nonmatching grids[END_REF]. As the solution goes very fast to zero on Ω \ Ω L we can expect here to see high order convergence due to the Bogner-Fox-Schmit FE space over Ω in . This is confirmed by the experiments (see Figure 8). We not only see higher order convergence [START_REF] Cai | Overlapping nonmatching grid mortar element methods for elliptic problems[END_REF] for an increasing number of elements in the interior rectangular mesh τ in (1st, 2nd and 3rd from the left) in comparison with ψ bdp -level lines for the solution obtained wit the standard method. Center: Table of the location of the magnetic axis and the plasma boundary defining point and the corresponding values for the solution without the MEM at the bottom line. Bottom: Some characteristic running times, the total run time (run), the time for the inversion of algebraic system in the n Newton Newton iterations (solve), the time for assembling j(•, •) in the n Newton Newton iterations (j assembling) and the numbers of degrees of freedom (V ex + V in ). w.r.t. refinement, but observe also that the computing time is much shorter than that for a low order method. This shows that the additional overhead due to coupling mappings is easily compensated. More numerical test for the MEM and the model problem ( 21) can be found in [START_REF] Truong | The mortar element method for the free-boundary toroidal plasma equilibrium problem[END_REF].

5.2 Movement of (r max , z max ) and (r bdp , z bdp )

In this example we are running a sequence of 30 simulations where all currents except one, namely the current in coil 4, are set to the values in the table of Figure 4. The current in coil 4 goes through uniform steps from I 4 = 28800 At to I 4 = 104965 At. This test case mimics the evolution of a plasma induced by the variation of currents in time. In the conforming case with linear Lagrangian FEs the magnetic axis (r max , z max ) and the plasma boundary defining point (r bdp , z bdp ) undergo a discontinuous evolution as their location is inherently restricted to vertices of the mesh. With the MEM we are able to introduce FE functions in Ω in that are not only continuous but have also continuous derivatives, hence the location of critical points is no more restricted to a finite number of points. This reasoning agrees perfectly with the observations. In Figure 9 we see that the evolution of the perturbations (∆r max , ∆z max ) and (∆r bdp , ∆z bdp ) of magnetic axis and boundary defining point evolve smoothly with the current perturbation. To highlight the influence of the continuous derivatives we compare the results with the MEM that uses bilinear FEs (Q1) instead of the bicubic Bogner-Fox-Schmit FEs. Maxima and minima of bilinear FE functions are again necessarily on vertices of the mesh, while saddle points can lie either on vertices or inside an element (see Figure 9 right). We would like to stress that the evolution of the perturbations ∆ψ max and ∆ψ bdp of the values of ψ at the magnetic axis and the boundary defining point is smooth in both cases. The discontinuous behavior of the location of critical points is not inherited to the values of ψ h at its critical points.

The visualization in Figure 10 stresses the undesired behavior that appears due to noncontinuous gradients and shows how this defect can be cured by using the Bogner-Fox-Schmit FE.

Geometric Coefficients

We are validating the computation of the geometric coefficients described in section 4.6 for the following elliptic data Inria and ψ bdp = 0 with parameters a = 1, b = 3, r max = 2, z max = 0 and ψ max = 2 for which we find closed form expression:

ψ E (r, z) = ψ max - (r -r max ) 2 a 2 + (z -z max ) 2 b 2 , (22) 2
g 1 r ,ψ E N (y) =|ψ max |abπ , g 1,ψ E N (y) =|ψ max |abπr max , g 1 r 2 ,ψ E N (y) = |ψ max |abπ r 2 max + a 2 ψ max y , g |∇ψ E | 2 r 2 ,ψ E N (y) =4π|ψ max | (b 2 -a 2 )r max (r max -r 2 max + a 2 ψ max y) -a 4 ψ max y a 3 b r 2 max + a 2 ψ max y . (23)
We project the elliptic data onto the Bogner-Fox-Schmit FE space defined over the domain [0.1, 4] × [-4.5, 4.5], then we compare the numerical values obtained for the geometric coefficients with the analytical ones. Figure 11 shows the convergence on a sequence of mesh refinements. We are computing numerically the values of the geometric coefficients at 40 equidistant values between 0 and 1 and monitor the maximal relative error. We obtain the expected high order convergence. Not very surprisingly the rate of convergence seems to correspond to the rate of convergence of the projection error for ψ E in the H 1 -norm. Compared to bilinear FEs we can achieve higher accuracy with a fixed number of unknowns. This example puts in evidence that with an appropriately chosen numerical algorithm for the approximation of the geometric coefficients, the accuracy is directly linked to the quality to the approximation of ψ. Next, the theoretical understanding of higher order FE methods on the one end and the results of the experiments in section 5.1 on the other hand, show that the MEM for the non-linear equilibrium problem with Bogner-Fox-Schmit FEs in some parts of the vacuum chamber allows to achieve a fixed accuracy with less degrees of freedoms compared to lower order FEs. So, we can conclude that the MEM with Bogner-Fox-Schmit FEs yields more accurate geometric coefficients than lower order methods. Finally, we perform a quantitative comparison of the geometric coefficients for the WEST application example. Figure 12 shows the geometric coefficients for a solution obtained by the MEM coupled either with the bilinear FE space or with the Bogner-Fox-Schmit FE space. Computations use the coarsest mesh from Figure 5 and we see that the coefficients based on MEM with the bilinear FE space suffer from small oscillations, that are due to lack of accuracy. The small oscillations near y = 0 disappear when using the Bogner-Fox-Schmit FEs.

Conclusions and Outlook

In the spirit of the hp-finite element method (hp-FEM) the most economical discretization, the one that gives the lowest error for the minimal number of unknowns, uses high polynomial degree and large elements in regions where the solution is very smooth and very small elements and low polynomial degree where the solution is singular. Regions where the smoothness deteriorates are for example the neighbourhoods of material corners (iron core, passive structures) or of discontinuous source terms (coils). It is not economical to use bicubic or even biquintic FEs there. It would be very detrimental for the efficiency of the method to increase the polynomial degree everywhere, while only on the plasma domain continuous derivatives are needed. This reasoning is even more important for iron core transformer tokamaks like WEST, where the vacuum chamber domain represents a small portion of the overall computational domain (see Figure 1). Note that the use of linear FEs everywhere outside the vacuum chamber is not very economical, but it is very flexible and easy to implement. It can be easily extended to a hp-FEM with adapted mesh sizes and polynomial degrees, once accuracy and computing time become Inria critical. As we deal here only with 2D problems and do not focus on realtime applications the computational time is currently not an issue. The implementation, not yet runtime optimized, suggests that the computational complexity of MEM with bicubic or biquintic FEs is only slightly larger than the previous approach with linear FEs everywhere. Bi-cubic FEs with isoparametric mappings everywhere, very common in fixed-boundary equilibrium calculations [START_REF] Huysmans | Isoparametric bicubic Hermite elements for solution of the Grad-Shafranov equation[END_REF][START_REF] Lütjens | Axisymmetric mhd equilibrium solver with bicubic hermite elements[END_REF], is not an option since we want that the meshes resolve the corners of the geometry (such as coils, passive structures, iron core). It should be possible to combine biquintic FEs in the vacuum chamber domain accessible by the plasma with linear FEs (or hp-FEM) in the exterior domain using techniques from discontinuous Galerkin methods that avoid overlapping meshes. It is a priori not clear which approach has better complexity but we would not expect much of a difference.

The MEM approach provides an approximation of the poloidal flux that is the solution to the discretized equilibrium problem and has continuous derivatives. It would be of course also possible, as suggested in [38, page 131], to improve a posteriori the smoothness of a numerical solution that has not continuous derivatives. But these improved solutions, also called reconstructions, in general do not solve a discretization of the equilibrium problem. Moreover, it can be fairly tricky to provide good reconstruction algorithms that avoid nonphysical oscillations. Rather then solving first a discretization of the equilibrium problem to obtain a bilinear FE approximation that is afterwards mapped onto piecewise bicubics FEs using, e.g., the method in [START_REF] Akima | A method of bivariate interpolation and smooth surface fitting based on local procedures[END_REF], we do solve directly for the unknown coefficients, the flux values and derivatives at the nodes.

We have shown that the MEM combining Cartesian and triangular meshes is a very flexible approach to introduce locally higher order regular FEs for plasma equilibrium calculations in tokamaks. While it is easy to define higher order regular FE spaces on Cartesian meshes, the triangular meshes in turn allow for an accurate resolution of design details in realistic geometries. It is possible to use similar ideas for the coupling of more complex plasma models with eddy current modeling in coils and passive structures. Moreover, we could enforce not only continuous derivatives but also continuity for higher order derivatives, using tensor products of higher order splines to define appropriate FE spaces in analogy to the Bogner-Fox-Schmit FE one. Figure 13 shows a numerical result that uses a biquintic FE space in Ω in that ensures continuity of second order derivatives. Continuity of second order derivatives goes beyond the scope of conforming FEs on unstructured meshes, but thanks to MEM we can provide now plasma equilibrium solvers with this feature. The main motivation for this work is more accurate computation of geometric coefficients and location of axis and boundary defining point for forthcoming simulations of the Grad/Hogan model. But there are many more applications in fusion science that can benefit from an MEM approach. The control of the location of the plasma boundary defining point for heat load minimization of divertor design [START_REF] Blommaert | A novel approach to magnetic divertor configuration design[END_REF][START_REF] Blommaert | Towards automated magnetic divertor design for optimal heat exhaust[END_REF] or the computation of plasma equilibria with so-called snowflake configuration [START_REF] Ryutov | Geometrical properties of a snowflake divertor[END_REF] are two of such very relevant applications.

Figure 13: Numerical results with the data from Section 5 by the MEM (13) coupled with biquintic FEs in the domain Ω L instead of the bicubic Bogner-Fox-Schmit FE. Without any post-processing, continuous second order derivatives are immediately available. The calculation is based on the second choice of the combination of meshes given in Figure 5. We show pseudoheat plots of ψ, ∂ r ψ, ∂ z ψ, ∂ r,r ψ, ∂ r,z ψ and ∂ z,z ψ (from left to right, top to bottom). Sufficiently far away from the interface γ we have consistent and smooth second order information on the poloidal flux. The non-smooth artificial interface γ introduces localized oscillations on the second order derivatives, only. Indeed, we are combining linear (vanishing second order derivatives) with bi-quintic (non-vanishing cubic second order derivatives) FEs.

Inria

Figure 1 :

 1 Figure 1: Left: Geometric description of the tokamak in the poloidal plane. Middle and right: Sketch for characteristic plasma shapes. The plasma boundary touches the limiter (middle) or the plasma is enclosed by a flux line that goes through an X-point (right).

Figure 2 :

 2 Figure 2: Left: A sketch of the computational domain Ω (ρ Γ = 5.8m) and of the different domains for the tokamak WEST. Orange, light blue and green domains are the Ω L bounded by the limiter, the coils Ω ci and the iron domain Ω Fe . The red domains, the passive structures will be important for transient simulations. Middle and Right: Two close-ups to the vacuum chamber domain with Ω L and the two meshes τ ex and τ in . The black and the red lines are γ and γ. The blue stars correspond to the degrees of freedom of V ex∂ . The red crosses and green circles signify the degrees of freedom of V in ∂ . The green circles are the nodes in N z γ , the nodes of γ that belong to an edge of γ that is parallel to the r-axis, and the red crosses are the nodes in N r γ , the nodes of γ that belong to an edge of γ that is parallel to the z-axis. Both red crosses and green circles belong to N γ .

Figure 3 :

 3 Figure 3: Nodal Interpolation (left): The node (r 3 , z 3 ) of the triangular mesh τ ex is in N γ and located in the rectangle {1 2 3 4}: it leads to 16 entries in P ex ∂,• and P ex ∂,∂ , corresponding to the values of the 16 basis functions of the Bogner-Fox-Schmit FE space not vanishing on the rectangle {1 2 3 4}. The node (r 5 , z 5 ) of the Cartesian mesh τ in is in N γ and N r γ and located in the triangle {4 9 10}: it leads to 6 entries in P in ∂,• and P in ∂,∂ , corresponding to the values of the 3 basis functions of the linear Lagrangian FE space not vanishing on the triangle {4 9 10}. L 2 -projection (right): See Remark 2.

Figure 4 :

 4 Figure 4: The values of the current (in At) for then numerical tests with the WEST tokamak. The small coils close to the Ω L are numbered 10 to 16 going from left to right and top to bottom.

Figure 5 :

 5 Figure 5: Contour plots of the solution by the MEM (13) and Bogner-Fox-Schmit (BFS) elements for an increasing number of elements in the interior rectangular mesh τ in (1st, 2nd and 3rd from the left) and a solution without MEM but with linear Lagrangian finite elements.

Figure 7 :Figure 8 :

 78 Figure 7: Error decay in H 1 (Ω)-norm of the MEM for the linear model problem (21) using either L 2 -projection (left) or nodal interpolation (right) and different mesh size ratios h in /h ex .

InriaFigure 9 :

 9 Figure 9: Evolution of the perturbations (∆r max , ∆z max ) and (∆r bdp , ∆z bdp ) of the magnetic axis and the boundary defining point w.r.t. the current perturbation ∆I 4 : the MEM (13) with the Bogner-Fox-Schmit element (left) and the MEM with bilinear finite elements instead of the Bogner-Fox-Schmit elements (right).

Figure 11 :Figure 12 :

 1112 Figure 11: Validation example: The convergence of the geometric coefficients for sequence of mesh refinements with bilinear FEs (left) and Bogner-Fox-Schmit FEs (right). We monitor the maximal relative error for g c,ψ E N (y) with c = c 1 = 1/r (blue), c = c 2 = 1 (red), c = c 3 = 1/r 2 (yellow) and c = c 4 = |∇ψ E | 2 /r 2 (purple), where the maximum is taken over 40 equidistant values of y between 0 and 1. The abscissa is the number N of degrees of freedom of the FE description of ψ,

  max (m) z max (m) ψ max (W b) r bdp (m) z bdp (m) ψ bdp (W b) Top: the ψ bdp -level lines of the solution by the MEM

		all P1			all P1		all P1
		P1/BFS			P1/BFS		P1/BFS
	P1 [33]	2.5879	0.0026	1.3968	2.2795	-0.6125	1.1652
	BFS 1	2.5843	0.0073	1.3949	2.2646	-0.6197	1.1614
	BFS 2	2.5799	0.0160	1.3890	2.2726	-0.6174	1.1582
	BFS 3	2.5842	0.0163	1.3959	2.2794	-0.6158	1.1641
		number of degrees of run(s)	solve	j assembling
		freedom V ex + V in		(s/n Newton ) (s/n Newton )
	P1 [33]	14008 + 4336	6.82	2.14/11	0.69/11
	BFS 1	14060 + 3664	6.52	2.81/8	1.91/8
	BFS 2	14060 + 14064	11.47	6.29/7	2.44/7
	BFS 3	14060 + 53344	45.43	30.85/7	5.57/7
	Figure 6:
							Inria

r

  Locations of (r max , z max ) (left) and (r bdp , z bdp ) (right) for variations ∆I 4 of the current in coil 4. Simulations are done on a mesh (top) and on a refinement of it (bottom). The solution with bilinear (Q1) FEs is compared with the numerical solution using the Bogner-Fox-Schmit (BFS) FEs.

	0.04										-0.61
									Q1							Q1
									BFS							BFS
	0.02
											-0.62
	0
	-0.02										-0.63
	-0.04
											-0.64
	-0.06
	-0.08										-0.65
	-0.1
											-0.66
	-0.12
	-0.14										-0.67
	.52	2.54	2.56	2.58	2.6	2.62	2.64	2.66	2.68	2.7	2.24	2.25	2.26	2.27	2.28	2.29	2.3
	0.04										-0.61
									Q1							Q1
									BFS							BFS
	0.02
											-0.62
	0
	-0.02										-0.63
	-0.04
											-0.64
	-0.06
	-0.08										-0.65
	-0.1
											-0.66
	-0.12
	-0.14										-0.67
	2.52	2.54	2.56	2.58	2.6	2.62	2.64	2.66	2.68	2.7	2.24	2.25	2.26	2.27	2.28	2.29	2.3
	Figure 10:

RR n°8916

http://www-sop.inria.fr/members/Holger.Heumann/Software.html Inria

an ASCII-file with the precise definition can be provided upon request RR n°8916

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments

The authors would like to thank the anonymous reviewers for their helpful comments and are grateful to Jacques Blum and Cedric Boulbe for stimulating discussions on plasma equilibrium problems. The second author warmly thanks the CASTOR team at INRIA Sophia-Antipolis for the delegation in 2015/2016, during which this work was completed.