# Tokamak free-boundary plasma equilibrium computations in presence of non-linear materials

Cedric Boulbe, Blaise Faugeras, Guillaume Gros, Francesca Rapetti

► To cite this version:

Cedric Boulbe, Blaise Faugeras, Guillaume Gros, Francesca Rapetti. Tokamak free-boundary plasma equilibrium computations in presence of non-linear materials. Journal of Scientific Computing, 2023, 96, 10.1007/s10915-023-02265-8. http://dx.doi.org/10.1007/s10915-023-02269-3

HAL Id: hal-03423469

https://hal.science/hal-03423469v3

Submitted on 6 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Tokamak free-boundary plasma equilibrium computations in presence of non-linear materials

Cédric Boulbe · Blaise Faugeras · Guillaume Gros · Francesca Rapetti

Received: date / Accepted: date

**Abstract** We consider the axisymmetric formulation of the equilibrium problem for a hot plasma in a tokamak. We adopt a non-overlapping mortar element approach, that couples $C^0$ piece-wise linear Lagrange finite elements in a region that does not contain the plasma and $C^1$ piece-wise cubic reduced Hsieh-Clough-Tocher finite elements elsewhere, to approximate the magnetic flux field on a triangular mesh of the poloidal tokamak section. The inclusion of ferromagnetic parts is simplified by assuming that they fit within the axisymmetric modeling and a new formulation of the Newton algorithm for the problem solution is stated, both in the static and quasi-static evolution cases.

**Keywords** Tokamak · equilibrium · non-linear materials · reduced Hsieh-Clough-Tocher finite element · mortar element method · Newton method

**Mathematics Subject Classification (2020)** 65N30 · 65D05

## 1 Introduction

Theoretical and computational plasma physics is a wide subject with applications ranging from low temperature plasmas for lighting, thrusters and materials processing to hot plasmas for fusion; from ultra-cold plasmas to particle accelerators; from beams to pulsed power; and from intense kinetic non-equilibrium plasmas to high power microwaves. Each application is characterized by a proper space-time scaling, mathematical model and computational approach. In this work, we are interested in simulating the equilibrium of a plasma for fusion reaction in a tokamak [4]. We push forward the method proposed in [10] to compute a plasma equilibrium in tokamak devices that include ferromagnetic parts. The choice of an iron-transformer tokamak is due to Paul-Henri Rebut, a French physicist, working on nuclear fusion. From 1970 to 1973, Rebut contributed to the creation of TFR (Tokamak of Fontenay-aux-Roses), then of JET (Joint European Torus) and of

<footnote>Cédric Boulbe · Blaise Faugeras · Guillaume Gros · Francesca Rapetti
Dep. de Mathématiques, Univ. Côte d'Azur, CNRS, INRIA CASTOR Team, Nice (FR).
E-mail: cedric.boulbe, blaise.faugeras, guillaume.gros, francesca.rapetti at univ-cotedazur.fr</footnote>
Tore Supra (after the discontinuation of TFR). Tore Supra later became WEST (Tungsten (W) Environment in Steady-state Tokamak). The tokamak is a sort of huge transformer where the plasma current is the secondary circuit coupled to the primary one represented by the current in the coils that generate the poloidal field. In a tokamak with iron, the magnetic field lines are better conveyed (than by the air) leading to an increase in the poloidal flux thus generating a longer fusion reaction (at that time, the technology of supra-conducting coils to generate high intensity fields was not so well developed yet). However, the presence of iron makes numerical computations more involved. Indeed, the magnetic permeability $\mu$ depends non linearly on the magnetic induction. Moreover, the presence of iron parts (an internal kernel with an external arm) breaks the toroidal symmetry of the physical parameter distribution despite the plasma equilibrium being an axisym- metric phenomenon. For the TFR, that worked from 1973 to 1986, the presence of iron caused an instability on the horizontal displacement of the plasma, as de- scribed in [6,25]. Tokamaks of new generation, such as ITER (under construction in Cadarache, France), are iron-free: thanks to modern technologies, the magnetic induction in the plasma can easily reach 10 teslas (and this would have not been possible with iron parts saturating at lower intensities). Iron-tokamaks such as WEST and JET are still used by scientists to make experiments.

In this work, we rely on the existing literature (see, e.g., [4, 5, 10, 11, 16, 17]) to make some steps forwards. We approximate the solution of both the static and quasi-static evolution problems by adopting a non-overlapping domain decompo- sition mortar-like approach, coupling $C^0-C^1$ finite elements, on a triangular mesh of the poloidal section of a iron tokamak. Due to the presence of iron parts and to the mortar-like coupling, the Newton's algorithm adopted to solve the non-linear problems, has a new iteration, written in terms of the coupling and the Jacobian matrices.

We start in Section 2 by recalling from books such as [4,18] how to derive the famous Grad-Shafranov-Schlüter equation [14,19,26], to solve for the numer- ical simulation of the axisymmetric equilibrium of the plasma. Then, in Section 3, we treat the static case: we adapt the finite element approach proposed in [10] involving highly regular approximations of the poloidal magnetic flux field $\psi$, to the iron case. Already in [17], it has been remarked that finite elements providing piece-wise polynomial approximations of $\psi$ that are only $C^0$ have two main draw- backs: 1.) The definition of the plasma boundary hinges on the critical points of the unknown flux $\psi$. If the derivatives of $\psi$ are not continuous, these points will neither be correctly calculated nor move in a continuous way during the plasma evolution. Indeed, with classical piece-wise finite elements, critical points of $\psi$ are necessarily located at mesh nodes. 2.) The resistive diffusion and transport of the heat in plasma are described by one-dimensional equations containing metric coef- ficients that depend on the gradient of the solution $\psi$ of the equilibrium problem. Many plasma characteristics (e.g. the so-called safety factor or the average cur- rent density profile), important to quantify stability or for monitoring during the experiment, are defined as integrals involving the gradient $\nabla\psi$ of the poloidal flux $\psi$ (see [4]). These coefficients are not well-defined if the gradients are not continu- ous. Differently to other approaches in the recent literature (see for example [17,21, 24]), to solve the axisymmetric formulation of the free-boundary plasma equilib- rium in a tokamak, we rely on a non-conforming domain decomposition approach
that couples $C^0$ piece-wise linear Lagrange finite elements in a region that does not contain the plasma and $C^1$ piece-wise cubic reduced Hsieh-Clough-Tocher [9] finite elements elsewhere. This approach gives the flexibility to achieve easily and at low cost higher order regularity for the approximation of the flux function $\psi$ in the domain covered by the plasma, thus resolving the cited drawbacks, while preserving accurate meshing of the geometric details in the rest of the computational domain. The continuity of the numerical solution at the coupling interface is weakly enforced by mortar projection [2]. We write the matrix problem and the modified Newton method to solve it. Section 4 is dedicated to the quasi-static evolution case: the new discrete problem includes a circuit equation (see [16], appendix A) and passive terms. The Newton iteration in the quasi-static case also counts out new terms. Numerical results for the static and quasi-static evolution cases are then presented<fnref n="1" /> in Section 5. Few conclusions in Section 6 end the paper.

## 2 The direct equilibrium problem

In a plasma for nuclear fusion, the charged particles (essentially, tritium and deuterium) at an extremely high temperature (ten times larger than that in the Sun) endure a fusion reaction, that is they stitch together, against the Coulomb repulsion, yielding production of energy, helium and neutrons. No material on Earth can support the temperature of such a hot mixture but due to the fact that the involved particles are charged, they can be confined in a toroidal chamber with magnetic field, tokamak in Russian. An additional iron structure can be installed in a tokamak to increase the poloidal flux thus generating a longer reaction. To keep up a fusion reaction we have, among many other tasks, to control the plasma in order to maintain it in equilibrium. A comprehensive survey of the (direct and inverse) mathematical problems associated with this equilibrium and of their low-order $C^0$ finite element modeling is described in [16] and the therein references. Here, we focus on the direct problem of computing a static equilibrium or the quasi-static evolution of a plasma in a tokamak by a mortar element approach coupling $C^0$ piece-wise linear with $C^1$ piece-wise cubic finite elements.

### 2.1 Mathematical properties for modeling the plasma at the equilibrium

The description of the plasma as a fluid that carries electrical currents and magnetic fields is surely simplified (e.g., kinetic effects are ignored) but it enables the derivation and understanding of some of its most basic properties. In particular, the equations of magneto-hydro-dynamics (MHD) may be used to describe how the magnetic configuration of a tokamak holds the plasma in equilibrium. These are the continuity, momentum and energy equations in the plasma domain for the volume charge density $\rho$ (with dimensions in the SI system<fnref n="2" /> as $[L]^{-3}[T][I]$), the

<footnote n="1">¹ All numerical simulations are here performed with the software NICE (see [11]).</footnote>

<footnote n="2">² In the Standard International (SI) unit system, mass $M$ (kg), length $L$ (m), time $T$ (s) and current intensity $I$ (A) are base dimensions (resp., units).</footnote>
fluid velocity $\mathbf{v}$ (as $[L][T]^{-1}$) and the pressure $p$ (as $[M][L]^{-1}[T]^{-2}$), respectively,

$$ 
\begin{aligned}
\partial_t \rho + \operatorname{div}(\rho \mathbf{v}) &= 0 && \text{(continuity equation)}, \\
\rho \frac{d\mathbf{v}}{dt} = \mathbf{J} \times \mathbf{B} - \operatorname{grad} p && \text{(momentum equation)}, \\
\frac{d}{dt} \left( \frac{p}{\rho \gamma} \right) &= S && \text{(energy equation)},
\end{aligned}
\quad (1) $$

where $S$ is a source collecting several terms, $\gamma > 1$ and $\frac{d}{dt}$ denotes the material time derivative, together with the magneto-quasi-static Maxwell's equations in the whole domain

$$ 
\begin{aligned}
\mathbf{E} + \mathbf{v} \times \mathbf{B} &= \mathbf{0} && \text{(ideal Ohm's law)}, \\
\operatorname{curl}\mathbf{E} = -\partial_t \mathbf{B} && \text{(Faraday's law)}, \\
\operatorname{curl}(\frac{1}{\mu}\mathbf{B}) &= \mathbf{J} && \text{(Ampère's theorem)}, \\
\operatorname{div}(\mathbf{B}) &= \mathbf{0} && \text{(solenoidality condition)}.
\end{aligned}
\quad (2) $$

for the electric field $\mathbf{E}$ (as $[M][L][T]^{-3}[I]^{-1}$), the magnetic induction field $\mathbf{B}$ (as $[M][T]^{-2}[I]^{-1}$), the current density $\mathbf{J}$ (as $[I][L]^{-2}$), where $\mu$ (as $[M][L][T]^{-2}[I]^{-2}$) the magnetic permeability and $\partial_t \mathbf{B}$ the time derivative of $\mathbf{B}$. Non-ferromagnetic parts of the tokamak have $\mu = \mu_0$, where $\mu_0$ is the magnetic permeability of the vacuum. Suitable boundary conditions close the MHD system. These conditions translate into mathematical terms the following facts: (1) the plasma is confined inside a perfectly conducting wall, (2) the wall separates the plasma from a vacuum region, and (3) the plasma is surrounded by external coils. We will detail these conditions when stating the final form of the problem to solve. The MHD system of equations (1), (2) is labeled as *ideal* since all resistive, viscous, conductive and diffusive terms have been neglected. Taking into account all these effects is mathematically and physically far from trivial (to this purpose, one can see the work done in [15]) and goes beyond the purpose of the present analysis.

In this work we consider the equilibrium of the plasma, we remain at the diffusion time scale (the slowest one) in a device with characteristic length of meters. The equations describing an ideal MHD equilibrium are respectively, force balance (between the kinetic force and the magnetic force), Ampère's theorem and the solenoidality condition, that are

$$ \operatorname{grad} p = \mathbf{J} \times \mathbf{B}, \quad \operatorname{curl}(\frac{1}{\mu}\mathbf{B}) = \mathbf{J}, \quad \operatorname{div}(\mathbf{B}) = 0. \quad (3) $$

We know that the distribution of iron structures in these tokamaks is not at all axisymmetric. Therefore, to fit within the axisymmetric modeling we make the following assumption: all cross-sections of the considered iron tokamak are identical. Moreover, when writing integrals, we will omit the integration element if this is not misleading.

Under the assumption of perfect axial symmetry of the device geometry and physical parameters' distribution, let $\mathbf{e}_\varphi$ is the unit vector for the toroidal coordinate $\varphi$ in the coordinate system $(r, \varphi, z)$. In these coordinates, $r$ measures the distance from the tokamak axis, $\varphi$ is the toroidal angle and $z$ is the height along the tokamak axis. The magnetic induction field $\mathbf{B}$ can be decomposed into the sum of $\mathbf{B}_t$, a vector in the same direction as $\mathbf{e}_\varphi$, and $\mathbf{B}_p$, a vector with a direction
orthogonal to $\mathbf{e}_\varphi$, and both $\mathbf{B}_t$ and $\mathbf{B}_p$ independent of $\varphi$. We express $\mathbf{B}$ in terms of the poloidal flux function $\psi$ and of another function $f$ as given in [4], that is

$$ \mathbf{B} = \frac{1}{r} \mathrm{grad} \psi \times \mathbf{e}_\varphi + \frac{f}{r} \mathbf{e}_\varphi = \mathbf{B}_p + \mathbf{B}_t. \quad (4) $$

The first term, $\mathbf{B}_p$, in (4) is the poloidal component of $\mathbf{B}$ that lies in the cross-section plane $(r, z)$ also called poloidal section ($\varphi = \text{constant}$) of the tokamak. The second term $\mathbf{B}_t$ in (4) is the toroidal component of $\mathbf{B}$ and $f$, referred to as the poloidal current flux, is such that $f \mathbf{e}_\varphi = r \mathbf{B}_t$. To fulfill the solenoidality condition on $\mathbf{B}$, a magnetic vector potential $\mathbf{A}$ such that $\mathbf{B} = \mathrm{curl} \mathbf{A}$ is introduced. The Coulomb gauge condition $\mathrm{div} \mathbf{A} = 0$, which is typically used when the propagation velocity of the perturbations of the magnetic field lines is smaller than the speed of light, is imposed on $\mathbf{A}$ to ensure uniqueness. We will see that, in an axisymmetric formulation as the one we consider here, the Coulomb gauge on $\mathbf{A}$ is automatically satisfied. In particular, $\mathbf{B}_p = \mathrm{curl} \mathbf{A}_t$ and the poloidal magnetic flux $\psi$ thus represents the scaled toroidal component of the vector potential $\mathbf{A}$, namely $\psi \mathbf{e}_\varphi = r \mathbf{A}_t$. We recall from books (see, e.g., [4, 18]) two important properties for the mathematical modeling of the plasma at the ideal MHD equilibrium. Their proofs are here detailed for completeness. Property 1 states that the poloidal magnetic flux $\psi$ is a key quantity in modeling plasma in tokamaks.

**Property 1.** The lines of both the current density $\mathbf{J}$ and magnetic induction $\mathbf{B}$ are on surfaces of constant value for $\psi$ (and $p$). They are called magnetic surfaces.

*Proof* In cylindrical coordinates we have $\mathrm{grad} p = (\partial_r p, \frac{1}{r} \partial_\varphi p, \partial_z p)^\top$. As we assume $\mu$ constant in the plasma domain, we obtain

$$ \mathrm{curl} \left( \frac{1}{\mu} \mathbf{B}_t \right) = \mathrm{curl} \left( \frac{f}{\mu r} \mathbf{e}_\varphi \right) = \frac{1}{\mu r} (-\partial_z f) \mathbf{e}_r + \frac{1}{\mu r} (\partial_r f) \mathbf{e}_z \quad (5) $$

being $\mathbf{e}_r$ and $\mathbf{e}_z$ the unit vectors for the poloidal coordinates $r, z$, respectively. Let us consider the force balance identity in (3), then

$$ 0 = \mathrm{grad} p \cdot \mathbf{B}_p = \mathrm{grad} p \cdot \mathrm{curl} \mathbf{A}_t = \frac{1}{r} \partial_r p (-\partial_z \psi) + \frac{1}{r} \partial_z p (\partial_r \psi) $$

since<fnref n="3" /> $\mathrm{curl} \mathbf{A}_t = \frac{1}{r}(-\partial_z \psi, 0, \partial_r \psi)^\top$. Hence, $p$ is constant over surfaces where $\psi$ is constant, so $p = p(\psi)$. Moreover, for the axisymmetry assumption,

$$ 0 = \mathrm{grad} p \cdot \mathbf{J}_p = \mathrm{grad} p \cdot \mathrm{curl} \left( \frac{1}{\mu} \mathbf{B}_t \right) = \frac{1}{\mu r} \partial_r p (-\partial_z f) + \frac{1}{\mu r} \partial_z p (\partial_r f) $$

and thus $f$ is constant over surfaces where $p$ and $\psi$ are constant, so $f = f(\psi)$. □

Magnetic field line twist while developing in the toroidal direction, and charged particles remain localized closely to the magnetic surfaces, while moving in the plasma mixture (see, e.g., [18] for more details on the magnetic field properties in plasmas). The helicity of a divergence-free vector field is a standard measure for the extent to which the field lines wrap and coil around one another (see,

<footnote n="3">³ From $\psi \mathbf{e}_\varphi = r \mathbf{A}_t$, we get $\mathbf{A}_t = (0, A_\varphi, 0)^\top$ where $A_\varphi$ is equal to $\frac{1}{r} \psi$. In cylindrical coordinates, for this vector $\mathbf{A}_t$ we have $\mathrm{curl} \mathbf{A}_t = (-\partial_z A_\varphi, 0, \frac{1}{r} \partial_r (r A_\varphi))^\top$.</footnote>
e.g., [7]). The connection between twists and knots is analyzed in [23,22] and the topological interpretation of helicity in terms of linking numbers is given, e.g., in [3,20] and references therein. The magnetic flux surfaces are hence strictly connected to the value of $\psi$, they are nested surfaces, that do not intersect with any material of the tokamak, and ensure the confinement of charged particles, namely the confinement of plasma inside a tokamak. The center of the plasma, where pressure has a maximum, is called the magnetic axis. The last surface in proximity of which charged particles can still move, without striking against the limiter or the divertor of the device, defines the domain containing the plasma, say $\Omega_p(\psi)$. Property 2 relates the poloidal magnetic flux $\psi$ to the current in the plasma domain $\Omega_p(\psi)$.

**Property 2.** The flux $\psi$ verifies the Grad-Shafranov-Schlüter equation in $\Omega_p(\psi)$.

*Proof* We have

$$ 
\begin{aligned}
(-\partial_z f \mathbf{e}_r + \partial_r f \mathbf{e}_z) &= f'(\psi) (-\partial_z \psi \mathbf{e}_r + \partial_r \psi \mathbf{e}_z) && (\text{as } f = f(\psi)) \\
&= f'(\psi) \operatorname{grad} \psi \times \mathbf{e}_\varphi && (\text{eq. (4)}) \\
&= f'(\psi) r \mathbf{B}_p \\
&= f'(\psi) r \operatorname{curl} \mathbf{A}_t.
\end{aligned}
$$

We consider the force balance equation in (3). We obtain

$$ \operatorname{grad} p = (\mathbf{J}_p + \mathbf{J}_t) \times (\mathbf{B}_t + \mathbf{B}_p) \quad \Rightarrow \quad \operatorname{grad} p = \mathbf{J}_p \times \mathbf{B}_t + \mathbf{J}_t \times \mathbf{B}_p $$

as

$$ \mathbf{J}_p \times \mathbf{B}_p = \operatorname{curl}(\frac{1}{\mu} \mathbf{B}_t) \times \mathbf{B}_p = \frac{1}{\mu} f'(\psi) \operatorname{curl} \mathbf{A}_t \times \operatorname{curl} \mathbf{A}_t = \mathbf{0} $$

and $\mathbf{J}_t \times \mathbf{B}_t = \mathbf{0}$ by using vector identities. Let us now consider the identity

$$ \operatorname{grad} p \cdot (\mathbf{e}_\varphi \times \operatorname{curl} \mathbf{A}_t) = (\mathbf{J}_p \times \mathbf{B}_t + \mathbf{J}_t \times \operatorname{curl} \mathbf{A}_t) \cdot (\mathbf{e}_\varphi \times \operatorname{curl} \mathbf{A}_t). \quad (6) $$

The magnetic surfaces are defined by a constant value of $p$ and $p = p(\psi)$, as we have seen with Property 1, hence $\operatorname{grad} p = f'(\psi) \operatorname{grad} \psi$. Since $\operatorname{curl} \mathbf{A}_t = \frac{1}{r}(-\partial_z \psi, 0, \partial_r \psi)^\top$, we have $(\mathbf{e}_\varphi \times \operatorname{curl} \mathbf{A}_t) = \frac{1}{r}(\partial_r \psi, 0, \partial_z \psi)^\top$ and the left-hand side of (6) gives

$$ p'(\psi) \operatorname{grad} \psi \cdot (\mathbf{e}_\varphi \times \operatorname{curl} \mathbf{A}_t) = \frac{1}{r} p'(\psi) |\operatorname{grad} \psi|^2. $$

For the terms in the right-hand side of (6) we have

$$ 
\begin{aligned}
(\mathbf{J}_p \times \mathbf{B}_t) \cdot (\mathbf{e}_\varphi \times \operatorname{curl} \mathbf{A}_t) \\
&= (\operatorname{curl}(\frac{1}{\mu} \mathbf{B}_t) \times \mathbf{B}_t) \cdot (\mathbf{e}_\varphi \times \operatorname{curl} \mathbf{A}_t) \\
&= (\frac{1}{\mu} f'(\psi) \operatorname{curl} \mathbf{A}_t \times \frac{f}{r} \mathbf{e}_\varphi) \cdot (\mathbf{e}_\varphi \times \operatorname{curl} \mathbf{A}_t) \\
&= -\frac{1}{\mu r} f f'(\psi) |\frac{1}{r} \operatorname{grad} \psi|^2 = -\frac{1}{\mu r^3} f f'(\psi) |\operatorname{grad} \psi|^2.
\end{aligned}
$$
Let us denote by $j_{\varphi} := \mathbf{J}_t \cdot \mathbf{e}_{\varphi}$ the toroidal component of the current density $\mathbf{J}$. The vector identity $(a \times b) \cdot (\alpha \times \beta) = (a \cdot \alpha)(b \cdot \beta) - (a \cdot \beta)(b \cdot \alpha)$ for generic vectors $a, b, \alpha, \beta$, yields

$$ 
\begin{aligned}
(\mathbf{J}_t \times \mathbf{B}_p) \cdot (\mathbf{e}_{\varphi} \times \operatorname{curl} \mathbf{A}_t) \\
&= (\mathbf{J}_t \times \operatorname{curl} \mathbf{A}_t) \cdot (\mathbf{e}_{\varphi} \times \operatorname{curl} \mathbf{A}_t) \\
&= (\mathbf{J}_t \cdot \mathbf{e}_{\varphi}) |\operatorname{curl} \mathbf{A}_t|^2 - (\mathbf{J}_t \cdot \operatorname{curl} \mathbf{A}_t) (\operatorname{curl} \mathbf{A}_t \cdot \mathbf{e}_{\varphi}) \\
&= j_{\varphi} \frac{1}{r^2} |\operatorname{grad} \psi|^2
\end{aligned}
$$

since $\operatorname{curl} \mathbf{A}_t \cdot \mathbf{e}_{\varphi} = 0$. We have thus obtained

$$ \frac{1}{r} p'(\psi) |\operatorname{grad} \psi|^2 = \left[ \frac{1}{r^2} j_{\varphi} - \frac{1}{\mu r^3} f f'(\psi) \right] |\operatorname{grad} \psi|^2. $$

We can simplify by $|\frac{1}{r} \operatorname{grad} \psi|^2$ in the three terms of the equation above and

$$ j_{\varphi} = \frac{1}{\mu r} f f'(\psi) + r p'(\psi) \qquad (7) $$

is the current in the plasma domain $\Omega_p(\psi)$. In cylindrical coordinates, the left-hand side of (7) in terms of $\psi$ becomes $-\operatorname{div}(1/(\mu r) \operatorname{grad} \psi) := -\Delta^* \psi$ and thus it holds

$$ -\Delta^* \psi = \frac{1}{\mu r} f f'(\psi) + r p'(\psi) \qquad (8) $$

that is the Grad-Shafranov-Schlüter equation for $\psi$ in $\Omega_p(\psi)$. $\square$

For an air-transformer tokamak, $\mu = \mu_0$ everywhere and $-\Delta^*$ in equation (8) is a linear second-order elliptic operator. For an iron-transformer tokamak, $\mu$ is a given function of $|\mathbf{B}_p|^2$, thus of $|\frac{1}{r} \operatorname{grad} \psi|^2$, in the ferromagnetic region and $-\Delta^*$ in equation (8) becomes a non-linear second-order elliptic operator. The main challenges for solving equation (8) numerically are its formulation on an infinite domain, the non-linear right-hand side, the non-linear permeability in iron and the non-linearity due to the free plasma boundary. In the following, we rely on existing literature, as for example [16, 17, 10, 11], and point out the modifications introduced by a non-conforming mortar-like approach to the whole resolution algorithm.

## 2.2 The continuous problem in the poloidal section

We start by recalling the formulation of the plasma equilibrium evolution problem as stated in [4, 5, 16, 18]. We introduce $\mathcal{D} = [0, \infty] \times [-\infty, \infty]$, the positive half plane that contains the poloidal section. The geometry of the tokamak determines various sub-domains (see Fig. 1, left) that are used to set the expression of $j_{\varphi}$ accordingly.

- $\Omega_{Fe} \subset \mathcal{D}$ denotes those parts of $\mathcal{D}$ made of iron that do not carry any current but where the magnetic permeability $\mu$ is not constant and depends non-linearly on $\psi$, namely $\mu(\psi, r, z) = \mu_{FE}(|\operatorname{grad} \psi|^2 r^{-2}) \ge \mu_0$; elsewhere, that is in $\Omega_{Fe}^c, \mu = \mu_0$;
<figure><img src="image_2.png" /><figcaption>**Fig. 1** (Taken from [16]) Left: Geometric description of the tokamak in the poloidal plane. Center and right: sketch for two characteristic plasma shapes. The plasma is depicted in red. Its boundary is a magnetic surface. The plasma is an extremely hot mixture made of charged particles that continue to move in the tokamak chamber under the effect of the magnetic forces. Those charged particles that leave the plasma, finish their run by striking into a material of the machine. The two presented configurations are those of a tokamak plasma where its boundary either touches the limiter (center) or goes through an X-point (right).</figcaption></figure>

- $\Omega_{c_i} \subset \mathcal{D}, 1 \le i \le N_c$, denotes the intersection with the poloidal plane of the $i$th field coil carrying currents and $N_c$ denotes the number of poloidal field coils used for confinement in the reactor. Here, the $i$th coil has cross section area $|\Omega_{c_i}|$ and carries a total current $I_i$ (which is a given constant, in the static modelization, or a function, solution of an electric circuit equation, in the quasi-static case);

- $S \subset \mathcal{D}$, denotes the part of domain $S = \cup_j \Omega_{ps_j}$, where are located the $N_{ps}$ passive structures $\Omega_{ps_j}$, characterised by an electric conductivity $\sigma_j \ne 0, j = 1, ..., N_{ps}$ (note that the tokamak vacuum vessel is a passive structure too);

- $\Omega_L \subset \mathcal{D}$, denotes the domain bounded by the limiter, thus the domain accessible by the plasma;

- $\Omega_p \subset \Omega_L$, denotes the domain covered by the plasma and the boundary $\partial\Omega_p$ is the outermost closed $\psi$-isocontour contained within $\Omega_L$.

The equilibrium of plasma in a tokamak has to satisfy at each instant the following non-linear initial boundary value problem: for each $t \in [0, T], T \ge 0$, find $\psi(r, z, t)$ such that

$$-\Delta^* \psi = j_{\varphi} = \begin{cases} rp'(\psi) + \frac{1}{\mu_0 r} ff'(\psi) & \text{in } \Omega_p(\psi), \\ I_i(\psi, t)/|\Omega_{c_i}| & \text{in } \Omega_{c_i}, i = 1, \dots, N_c, \quad (*) \\ -\frac{\sigma_j}{r} \partial_t \psi & \text{in } \Omega_{ps_j} \subset S, j = 1, \dots, N_{ps}, \\ 0 & \text{elsewhere,} \end{cases} \quad (9)$$

$$\psi(0, z, t) = 0 \quad \forall (0, z) \in \Gamma_0, \quad \lim_{||(r,z)||_2 \to +\infty} \psi(r, z, t) = 0,$$

$$\psi(r, z, 0) = \psi_0(r, z) \quad \forall (r, z) \in \mathcal{D} \quad \text{(initial condition)}$$

with more details about (*) in Section 5, where $I_i(\psi, t)$ is defined by (19). The plasma domain $\Omega_p(\psi)$ is unknown and depends non-linearly on the poloidal flux $\psi$
(we have a free-boundary problem). The boundary of $\Omega_p(\psi)$ either touches that of $\Omega_L$ (limiter configuration, as in Fig. 1 middle) or contains one or more saddle points of $\psi$ (divertor configuration, as in Fig. 1 right). The saddle points of $\psi$, denoted by $(r_x, z_x)=(r_x(\psi), z_x(\psi))$, are called X-points of $\psi$. The plasma domain $\Omega_p(\psi)$ is the largest sub-domain of $\Omega_L$ bounded by a closed $\psi$-isoline in $\Omega_L$ and containing the magnetic axis $(r_a, z_a)$. The magnetic axis is the point $(r_a, z_a) = (r_a(\psi), z_a(\psi))$, where $\psi$ has in $\Omega_L$ its global maximum (or minimum, depending on axis positive direction). Let $(r_b, z_b) = (r_b(\psi), z_b(\psi))$ be the point that determines the plasma boundary. Note that $(r_b, z_b)$ is either an X-point of $\psi$ or the contact point with $\partial\Omega_L$. As explained in [4,5,15], to determine the two functions $p'$ and $f f'$ it is necessary to complete (9) with additional (transport and diffusion) equations. In this work, we will assume that, up to some scaling coefficient $\lambda$, the functions $p'$ and $f f'$ are known. We thus suppose that we are given two polynomial or piece-wise polynomial functions $A(\cdot)$ and $B(\cdot)$ defined on $[0, 1]$ such that, in the plasma domain $\Omega_p(\psi)$, we have $j_\varphi = \lambda(\frac{r}{r_0}A(\psi_N) + \frac{r_0}{r}B(\psi_N))$. Here, $r_0$ the characteristic major radius (in meters) of $\Omega_L$ and $\psi_N(\psi, r, z) = (\psi(r, z)-\psi_a(\psi))/(\psi_b(\psi)-\psi_a(\psi))$ is the normalized poloidal flux. The domain of $p'$ and $f f'$ is the interval $[\psi_a, \psi_b]$ (supposing $\psi_a < \psi_b$) with the scalar values $\psi_a$ and $\psi_b$ being the flux values at the magnetic axis and at the boundary of the plasma. To solve numerically problem (10) we work in a domain $\Omega \subset D$, known as the ABB domain, named after Albanese-Blum-Barbieri [1], who first introduced it, associated with $D$ (see Fig. 2, left), delimited by a half-circle $\gamma$ of radius $\rho_\gamma > 0$ including $\Omega_L \cup \Omega_{Fe} \cup_i \Omega_{c_i} \cup_j \Omega_{ps_j}$ and the vertical segment $\Gamma_0 = \{0\}_r \times [-\rho_\gamma, \rho_\gamma]_z$.

### 3 The static problem

Let us first consider the static equilibrium problem: find $\psi(r, z)$ such that

$$ -\Delta^* \psi = j_\varphi = \begin{cases} rp'(\psi) + \frac{1}{\mu_0 r} f f'(\psi) & \text{in } \Omega_p(\psi), \\ I_i/|\Omega_{c_i}| & \text{in } \Omega_{c_i}, \quad i = 1, \dots, N_c, \\ 0 & \text{elsewhere,} \end{cases} \quad (10) \\ \psi(0, z) = 0 \quad \forall (0, z) \in \Gamma_0, \quad \lim_{||(r,z)||_2 \to +\infty} \psi(r, z) = 0. $$

Problem (10) results from (9) when we neglect the effects due to passive structures (that is $\sigma_j = 0$ for all $j = 1, \dots, N_{ps}$) and $I_i$ denotes the total current (in Ampère turns) in the $i$th coil, set independently from $\psi$.

#### 3.1 The weak form

To introduce the non-overlapping domain decomposition framework, we set $\Omega = \Omega^{\text{in}} \cup \Omega^{\text{ex}}$ where $\Omega^{\text{in}}$ is a bounded domain containing $\Omega_L$ (see Fig. 2, right) and $\Omega^{\text{ex}} = \Omega \setminus \Omega^{\text{in}}$. The boundary of $\Omega^{\text{in}}$ is denoted $\mathcal{I}$, to recall that it is an interface between the two sub-domains $\Omega^{\text{in}}$, $\Omega^{\text{ex}}$, on which we will impose the continuity of $\psi$, in a weak sense, through a mortar-like $L^2$ projection [2]. The weak formulation of (10) is: Find $\psi = (\psi_{\text{ex}}, \psi_{\text{in}}) \in \mathcal{V}$ such that

$$ a(\psi, s) := a^{\text{ex}}(\psi^{\text{ex}}, v) + a^{\text{in}}(\psi^{\text{in}}, w) = \ell(I, s) \quad \forall s = (v, w) \in \mathcal{V}_{0,\mathcal{I}} \quad (11) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2 The ABB domain (left) associated with $\mathcal{D}$ and a zoom (right) on the sub-domain $\Omega^{\text{in}}$ involved in the domain decomposition formulation. The mortar interface $\mathcal{I}$ in this case is closed. Example of a passive structure $S$ where $\sigma \neq 0$.</figcaption></figure>

where $\mathcal{V} = \{(v, w) \in \mathcal{H}^1(\Omega^{\text{ex}}) \times \mathcal{H}^1(\Omega^{\text{in}}), v|_{\Gamma_0} = 0, v|_{\mathcal{I}} = w|_{\mathcal{I}}\}$, being $\mathcal{H}^1(\Omega)$ the functional space defined as $\mathcal{H}^1(\Omega) = \{u \in L^2_*(\Omega), \nabla u \in L^2_*(\Omega)^2\}$ with $L^2_*(\Omega) = \{g : \Omega \to \mathbb{R}, \|g\|_{*,\Omega}^2 := \int_{\Omega} (\frac{g}{r})^2 r dr dz < \infty\}$ and $\nabla$ denoting the gradient operator in the poloidal variables. We have also set $\mathcal{V}_{0,\mathcal{I}} = \{(v, w) \in \mathcal{V}, v|_{\mathcal{I}} = w|_{\mathcal{I}} = 0\}$ and

$$ 
\begin{align} 
a^{\text{ex}}(\psi, v) &:= \int_{\Omega^{\text{ex}}} \frac{1}{\mu r} \nabla \psi \cdot \nabla v \, drdz + c(\psi, v), \nonumber \\ 
a^{\text{in}}(\psi, w) &:= \int_{\Omega^{\text{in}}} \frac{1}{\mu_0 r} \nabla \psi \cdot \nabla w \, drdz - J_p(\psi, w), \nonumber \\ 
J_p(\psi, w) &:= \int_{\Omega_p(\psi)} \lambda \left( \frac{r}{r_0} A(\psi_N) + \frac{r_0}{r} B(\psi_N) \right) w \, drdz, \tag{12} \\ 
\ell(I, s) &:= \sum_{i=1}^{N_c} \frac{I_i}{|\Omega_{c_i}|} \int_{\Omega_{c_i}} (\chi_{\Omega^{\text{ex}}} v + \chi_{\Omega^{\text{in}}} w) \, drdz \nonumber 
\end{align} 
$$

with $\lambda$ a scaling coefficient such that the total plasma current is

$$ I_p = \lambda \left| \int_{\Omega_p(\psi)} \left( \frac{r}{r_0} A(\psi_N) + \frac{r_0}{r} B(\psi_N) \right) dr dz \right| $$

We recall that, in the expression of $a^{\text{ex}}(\psi, v)$, the magnetic permeability can depend on $\psi$ as follows

$$ \mu = \mu_{FE} \left( \frac{|\nabla \psi|^2}{r^2} \right) \chi_{\Omega_{Fe}} + \mu_0 \chi_{\Omega^{\text{ex}} \setminus \Omega_{Fe}}, $$

where $\chi_V$ is the characteristic function of a set $V$ and $\mu_{FE}$ is a given function. Note that $\ell(I, s)$ contains the expression $\chi_{\Omega^{\text{ex}}} v + \chi_{\Omega^{\text{in}}} w$ to deal with the presence of coils in $\Omega^{\text{in}}$ and $\Omega^{\text{ex}}$. Moreover, $c(\psi, v) = \int_{\gamma} v \partial_n \psi d\Gamma$ takes into account the condition at infinity on $\gamma$ and will be discretized as explained in [16,17]. Under suitable assumptions, such as for example $\Omega_{Fe} = \emptyset$ or $I_p$ assigned, it can be proven that problem (11) has a unique solution [13], in the general case the question is theoretically open$^4$. In the next section, we propose a Newton method to solve

<footnote>$^4$ We wish to recall the fundamental contribution of Roland Glowinski to the analysis, the finite element approximation and numerical resolution by Newton-like methods of such nonlinear problems.</footnote>
the discrete problem associated with (11) when $\Omega_{Fe} \neq \emptyset$ and $\mu$ needs to be estimated from experimental data. The reconstruction of the function $\mu$ for a given ferromagnetic material is performed on the tokamak in absence of plasma. The function $\mu_r(H_p)$ representing the relative magnetic permeability ($\mu = \mu_0 \mu_r$) is experimentally determined as a function of the modulus $H_p$ of the poloidal magnetic field $H_p$. We anticipate that, as remarked in [12], if $\mu$ was directly linked to $B_p$, the Newton algorithm generally used to solve the final discrete problem could be divergent, as $B_p^2$ and thus $\mu$ vary significantly from one iteration to another. Therefore, the function $\mu_r(H_p)$ is first reconstructed by relying on the Ampère theorem and then, at each iteration $n$, we use

$$ \mu^n = g \left( \frac{|\nabla \psi_h^n|^2}{r^2 (\mu^{n-1})^2} \right). $$

To define the corresponding Jacobian matrix, we compute the derivatives w.r.t. the unknown field $\psi$ of the non-linear operators in (11), and then we evaluate them on discrete fields with special care. By involving directional derivatives, we can define $D_\psi a^{\text{ex}}(.,.)$, the differentiation operator w.r.t. $\psi$ in the direction of $\tilde{\psi}$, as follows

$$ D_{\psi} a^{\text{ex}}(\psi, s)(\tilde{\psi}) = a^{\text{ex}}(\tilde{\psi}, s) - 2 \int_{\Omega_{Fe}} \frac{g'(.)}{g^2(.)} \frac{1}{r^3} (\nabla \tilde{\psi} \cdot \nabla \psi) (\nabla \psi \cdot \nabla s) \quad (13) $$

where $(.)$ stands for $(|\nabla \psi|^2/r^2/(\mu^{n-1})^2)$. For $D_\psi a^{\text{in}}(.,.)$, the derivative w.r.t. $\psi$ of $J_p(.,.)$ is computed<fnref n="5" /> analytically on an approximation of this functional by a quadrature formula.

## 3.2 The discrete static problem

A mortar finite element approach is applied to (11) to get the discrete problem. Let $\tau^{\text{ex}}$ (resp. $\tau^{\text{in}}$) be a mesh of triangles that covers $\Omega^{\text{ex}}$ (resp. $\Omega^{\text{in}}$). The two meshes $\tau^{\text{ex}}$, $\tau^{\text{in}}$ are shape regular and quasi-uniform, with maximal diameters $h_{\text{ex}}$, $h_{\text{in}}$, respectively. We assume that $\mathcal{I}$ is a polygon with nodes and edges in $\tau^{\text{ex}}$. We wish to use in $\Omega_L \subset \Omega^{\text{in}}$, a finite element approximation $\psi_h$ for the poloidal flux $\psi$ that is not only continuous but has also component-wise continuous gradient $\nabla \psi_h$. This is possible if we use the piece-wise cubic reduced or minimal Hsieh-Clough-Tocher (rHCT) finite element space, say $\mathcal{V}^{\text{in}}$, on $\tau^{\text{in}}$ (see [9]). This regularity is not necessary in $\Omega^{\text{ex}}$ therefore we couple rHCT finite elements in $\Omega^{\text{in}}$ with continuous piece-wise linear finite elements, say $\mathcal{V}^{\text{ex}}$, on $\tau^{\text{ex}}$. The finite element space over the mesh $\tau^{\text{ex}}$ is $\mathcal{V}^{\text{ex}} = \{v \in \mathcal{C}^0(\Omega^{\text{ex}}), v|_{\Gamma_0} = 0, v|_T \in \mathbb{P}_1(T), \forall T \in \tau^{\text{ex}}\}$, whereas over $\tau^{\text{in}}$ is $\mathcal{V}^{\text{in}} = \{w \in \mathcal{C}^1(\Omega^{\text{in}}), w|_T \in P_{\text{loc}}(T), \forall T \in \tau^{\text{in}}\}$. The space $P_{\text{loc}}(T)$ reads

$$ P_{\text{loc}}(T) = \{w \in \mathcal{C}^1(T), w|_{B_i} \in \mathbb{P}_3(B_i), (\partial_n w)|_{B_i} \in \mathbb{P}_1(B_i), \forall b_i \in \partial B_i \cap \partial T\}, $$

with the triangle $T = [V_1, V_2, V_3]$ cut into three triangles $B_i = [G, V_m, V_\ell]$, having vertices in $V_m, V_\ell$ with $m, \ell \in \{1, 2, 3\} \setminus \{i\}$ and at the barycenter $G$, for each

<footnote n="5">⁵ To approximate $J_p(.,.)$ in (12) by a quadrature formula we need to know the domain $\Omega_p(\psi)$ occupied by the plasma. This domain is an unknown of the equilibrium problem, as it depends on $\psi$. An efficient technique to determine it is stated in [10].</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Barycentric subdivision of $T = [V_1, V_2, V_3] \in \tau^{\text{in}}$ into three sub-triangles $B_i$. Locally on $T$, at the three vertices $V_i$, we reconstruct the height $\psi_h(V_i)$ (black filled thick points) of $\psi$ and the tangent plane to the surface $\psi$ as generated by $\partial_r \psi_h(V_i)$, $\partial_z \psi_h(V_i)$ (empty circles around the vertices). Here **rot** $\psi_h = (\partial_r \psi_h, -\partial_z \psi_h)^t$. If we compute **rot** $\psi_h$, we get the restriction of $r B_p$ to $T$ with continuous components (the two black filled thick points) at the $V_i$ and continuous normal component (the small arrows) at the $b_i = \partial B_i \cap \partial T$ (see details in [8]).</figcaption></figure>

$i = 1, 2, 3$ as shown in Fig. 3. In the $P_{\text{loc}}(T)$ space definition, $n$ is the outward normal vector to $\partial T$, $b_i$ the edge $\partial B_i \cap \partial T$ and $(\partial_n w)|_{b_i}$ the normal derivative of $w$ along the edge $b_i$.

Note that $\Omega_{Fe} \subset \Omega^{\text{ex}}$. Let us also write $\mathcal{V}^{\text{ex}} = \mathcal{V}_{\circ}^{\text{ex}} \oplus \mathcal{E}\mathcal{V}_{\circ}^{\text{ex}}$ and $\mathcal{V}^{\text{in}} = \mathcal{V}_{\circ}^{\text{in}} \oplus \mathcal{E}\mathcal{V}_{\circ}^{\text{in}}$, where, for example, $\mathcal{V}_{\circ}^{\text{ex}}$ (resp. $\mathcal{V}_{\circ}^{\text{ex}}$) is the subspace of $\mathcal{V}^{\text{ex}}$ described by basis functions associated with dofs at nodes in $\bar{\Omega}^{\text{ex}} \setminus \mathcal{I}$ (resp., $\bar{\Omega}^{\text{ex}} \cap \mathcal{I}$) and $\mathcal{E}$ denotes the extension by zero operator. The functions in $\mathcal{V}_{\circ}^{\text{ex}}$ and $\mathcal{V}_{\circ}^{\text{in}}$ have vanishing Dirichlet trace on $\mathcal{I}$. The discrete problem to solve reads: Find $\psi_h \in \mathcal{V}_h$ such that

$$a(\psi_h, s_h) = \ell(I, s_h) \quad \forall s_h = (v_h, w_h) \in \mathcal{V}_{\circ}^{\text{ex}} \times \mathcal{V}_{\circ}^{\text{in}} \quad (14)$$

where

$$\mathcal{V}_h = \{(u_h^{\text{in}}, u_h^{\text{ex}}) \in \mathcal{V}^{\text{in}} \times \mathcal{V}^{\text{ex}}, u_{h|\Gamma_0}^{\text{ex}} = 0, \int_{\mathcal{I}} (u_h^{\text{in}} - u_h^{\text{ex}}) z_h \, d\mathcal{I} = 0, \forall z_h \in \mathcal{M}_h\},$$

with $\mathcal{M}_h = \{\xi_h \in \mathcal{C}^0(\mathcal{I}) : \xi_{h|e} \in \mathbb{P}_1(e), \forall e \in (\tau^{\text{ex}})|_{\mathcal{I}}\}$ the mortar multiplier space. The bilinear and linear forms $a(\cdot, \cdot)$, $\ell(I, \cdot)$ are defined as for the problem (11) and evaluated in (14) for functions in the discrete space $\mathcal{V}_h$.

### 3.3 The matrix problem and the Newton algorithm

Let us denote by $\{v_i^{\text{ex}}\}_{i=1,N^{\text{ex}}}$ the dual basis of $\mathcal{V}^{\text{ex}}$ for the P1 dofs associated with the $N^{\text{ex}}$ nodes $V_i \in \tau^{\text{ex}}$ and $\{w_j^{\text{in}}\}_{j=1,3N^{\text{in}}}$ that of $\mathcal{V}^{\text{in}}$ for the rHCT dofs at the $N^{\text{in}}$ nodes $V_j \in \tau^{\text{in}}$. Let $\mathbf{A}$ (resp. $\mathbf{C}$, $\mathbf{L}^{\text{in}}$, $\mathbf{L}^{\text{ex}}$) be the matrix associated with the integral expressions in (12) contained in $a(\cdot, \cdot)$ (resp., in $\mathbf{c}(\cdot, \cdot)$, in $\ell(\cdot, \cdot)$ for the coil $\Omega_{c_i}$ if this coil is in $\Omega^{\text{in}}$ or in $\Omega^{\text{ex}}$) and $\mathbf{J}(\cdot)$ (resp., $\mathbf{U}_I^{\text{in}}$, $\mathbf{U}_I^{\text{ex}}$) the vector with components resulting from $\mathbf{J}_p(\cdot)$ (resp., holding the currents $I_i$ for the coil $\Omega_{c_i}$ if it is in $\Omega^{\text{in}}$ or in $\Omega^{\text{ex}}$). To take into account iron parts in $\Omega^{\text{ex}}$, we separate the elliptic operator into the linear part and a nonlinear part, say $\mathbf{A}_0\psi + \mathbf{A}_{\mu}(\psi)$, where the vector $\psi$ gathers all dofs of $\psi_h \in \mathcal{V}_h$, the matrix $\mathbf{A}_0$ has entries

$$(\mathbf{A}_0)_{ij} = \int_{\Omega^{\text{ex}} \setminus \Omega_{Fe}} \frac{1}{\mu_0 r} \nabla v_i^{\text{ex}} \cdot \nabla v_j^{\text{ex}} \, dr \, dz, \quad i, j = 1, N^{\text{ex}},$$
and the vector $\mathbf{A}_{\mu}(\psi)$ has components

$$ (\mathbf{A}_{\mu}(\psi))_i = \int_{\Omega_{F_e}} \frac{1}{\mu(\psi_h)r} \nabla v_i^{\text{ex}} \cdot \nabla \psi_h dr dz, \quad i = 1, N^{\text{ex}}. $$

Equation (14) in its fully discretized form reads $\mathbf{e}(\psi) = \mathbf{0}$ with

$$ \mathbf{e}(\psi) := (\mathbf{A}_0 + \mathbf{C})\psi + \mathbf{A}_{\mu}(\psi) - \mathbf{J}(\psi) - \mathbf{L}^{\text{in}}\mathbf{U}_I^{\text{in}} - \mathbf{L}^{\text{ex}}\mathbf{U}_I^{\text{ex}} \quad (15) $$

where, for $k = 1, 3N^{\text{in}}$, we have

$$ (\mathbf{J}(\psi))_k = \mathbf{J}_p(\psi_{N,h}, w_k), \quad \text{and} \quad \mathbf{L}_{i,k}^{\text{in}} = \frac{1}{|\Omega_{c_i}|} \int_{\Omega_{c_i}} w_k^{\text{in}} dr dz, $$

for those indices $i = 1, N_c$ such that $\Omega_{c_i} \subset \Omega^{\text{in}}$. For the indices $i$ such that $\Omega_{c_i} \subset \Omega^{\text{ex}}$, the definition of $\mathbf{L}_{i,k}^{\text{ex}}$, with $j = 1, N^{\text{ex}}$, is similar to that of $\mathbf{L}_{i,k}^{\text{in}}$, just replacing $w_k^{\text{in}}$ by $v_j^{\text{ex}}$. Newton's iterations for problem (15) are

$$ \psi^{n+1} = \psi^n - [\mathbf{e}_\psi(\psi^n)]^{-1} \mathbf{e}(\psi^n), \quad (16) $$

with

$$ [\mathbf{e}_\psi(\psi)] = D_\psi[(\mathbf{A}_0 + \mathbf{C})\psi + \mathbf{A}_\mu(\psi)] - D_\psi \mathbf{J}(\psi). $$

Let $\mathbf{u}^{\text{ex}}$ and $\mathbf{u}^{\text{in}}$ gather the values of dofs for $\psi_h^{\text{ex}} \in \mathcal{V}^{\text{ex}}$ and $\psi_h^{\text{in}} \in \mathcal{V}^{\text{in}}$, respectively. We have $\mathbf{u}^{\text{ex}} = (\mathbf{u}_o^{\text{ex}}, \mathbf{u}_\partial^{\text{ex}})$ and $\mathbf{u}^{\text{in}} = (\mathbf{u}_o^{\text{in}}, \mathbf{u}_\partial^{\text{in}})$ where $\mathbf{u}_o^{\text{ex}}$ (resp. $\mathbf{u}_o^{\text{in}}$) and $\mathbf{u}_\partial^{\text{ex}}$ (resp. $\mathbf{u}_\partial^{\text{in}}$) are for dofs in $V_o^{\text{ex}}$ (resp. $V_o^{\text{in}}$) and $V_\partial^{\text{ex}}$ (resp. $V_\partial^{\text{in}}$). The mor- tar coupling condition in $\mathcal{V}_h$ links the block $\mathbf{u}_\partial^{\text{ex}}$ to the block $\mathbf{u}_\partial^{\text{in}}$ by the matrix relation $\mathbf{P}\mathbf{u}_\partial^{\text{ex}} = \mathbf{D}\mathbf{u}_\partial^{\text{in}}$ with $(\mathbf{P})_{i,j} = \int_{\mathcal{I}} v_{\partial,i}^{\text{ex}} v_{\partial,j}^{\text{ex}} d\mathcal{I}$, for all $i, j = 1, N_\partial^{\text{ex}}$, and $(\mathbf{D})_{i,k} = \int_{\mathcal{I}} v_{\partial,i}^{\text{ex}} w_{\partial,k}^{\text{in}} d\mathcal{I}$, for all $i = 1, N_\partial^{\text{ex}}$ and $k = 1, N_\partial^{\text{in}}$. The inclusion of the coupling condition matrix form into the algebraic system associated with the discrete problem (14) is done by introducing the reduced variable, $\mathbf{X}$, such that

$$ \psi = \begin{pmatrix} \mathbf{u}_o^{\text{ex}} \\ \mathbf{u}_\partial^{\text{ex}} \\ \mathbf{u}_o^{\text{in}} \\ \mathbf{u}_\partial^{\text{in}} \end{pmatrix} = \begin{bmatrix} \mathbf{I} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{P}^{-1}\mathbf{D} \\ \mathbf{0} & \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{I} \end{bmatrix} \begin{pmatrix} \mathbf{u}_o^{\text{ex}} \\ \mathbf{u}_o^{\text{in}} \\ \mathbf{u}_\partial^{\text{in}} \end{pmatrix} = \mathbf{Q}\mathbf{X}. $$

Equation (15) rewritten in terms of $\mathbf{X}$ becomes

$$ \mathbf{e}(\mathbf{X}) := \mathbf{Q}^{\top} [(\mathbf{A}_0 + \mathbf{C})\mathbf{Q}\mathbf{X} + \mathbf{A}_{\mu}(\psi) - \mathbf{J}(\psi) - \mathbf{L}^{\text{in}}\mathbf{U}_I^{\text{in}} - \mathbf{L}^{\text{ex}}\mathbf{U}_I^{\text{ex}}]. \quad (17) $$

For $\mathbf{J}(\psi) = \mathbf{J}(\mathbf{Q}\mathbf{X}) = \mathbf{H}(\mathbf{X})$ we get

$$ D\mathbf{X}\mathbf{H}(\mathbf{X})d\mathbf{X} = D_\psi\mathbf{J}(\psi)\mathbf{Q}d\mathbf{X} = \mathbf{J}\mathbf{a}_c\psi(\psi)\mathbf{Q}d\mathbf{X}, $$

with $\mathbf{J}\mathbf{a}_c\psi(\psi)$ the matrix representing the derivative of $\mathbf{J}(\psi)$ w.r.t. $\psi$. For the vector $\mathbf{A}_{\mu}(\psi) = \mathbf{A}_{\mu}(\mathbf{Q}\mathbf{X}) = \mathbf{G}(\mathbf{X})$, we obtain $D\mathbf{X}\mathbf{G}(\mathbf{X})d\mathbf{X} = \mathbf{A}_{\mu,\psi}(\psi)\mathbf{Q}d\mathbf{X}$ with $[\mathbf{A}_{\mu,\psi}]_{i,j}$ given in (13) by setting $\psi_h, v_i^{\text{ex}}, v_j^{\text{ex}}$ at the place of $\psi, \tilde{\psi}$ and $s$, respectively. Newton's iterations for problem (17) read

$$ \mathbf{X}^{n+1} = \mathbf{X}^n - [\mathbf{e}_{\mathbf{X}}(\mathbf{X}^n)]^{-1} \mathbf{e}(\mathbf{X}^n) \quad (18) $$

where

$$ [\mathbf{e}_{\mathbf{X}}(\mathbf{X}^n)] = \mathbf{Q}^{\top} [(\mathbf{A}_0 + \mathbf{C}) + \mathbf{A}_{\mu,\psi}(\psi^n) - \mathbf{J}\mathbf{a}_c\psi(\psi^n)]\mathbf{Q}. $$

In the next Section, we present the new iteration of the Newton's algorithm that is used to solve the discrete version of the quasi-static evolution problem in iron tokamaks. The evolution is quasi-static in the sense that it happens slowly enough for the system to go from one physical equilibrium to the successive.
## 4 The quasi-static evolution problem

In a tokamak, the poloidal field system is made of L circuits; each circuit, labeled by i with $1 \le i \le L$, includes $N_i$ coils (out of $N_c$) and $M_i$ power supplies (out of M). We refer to Appendix B in [5] and Appendix A in [16] for the technical details to get the circuit equation

$$ \vec{I}_i(\psi, t) = S_i \vec{V}_i(t) + R_i \vec{\Psi}_i(\partial_t \psi), \quad i = 1, \dots, L, \qquad (19) $$

expressing the current in the ith circuit. In (19), $\vec{I}_i$ is the vector of size $M_i + N_i$ containing the currents at the $M_i$ supplies and in the $N_i$ coils of the ith circuit; $\vec{V}_i \in \mathbb{R}^{M_i}$ is the one of the tensions applied to the supplies; finally, $\vec{\Psi}^i(\partial_t \psi) \in \mathbb{R}^{N_i}$ is the vector such that

$$ \vec{\Psi}^i(\partial_t \psi) = \left( \frac{1}{|\Omega_{c_i, 1}|} \int_{\Omega_{c_i, 1}} \partial_t \psi dr dz, \dots, \frac{1}{|\Omega_{c_i, N_i}|} \int_{\Omega_{c_i, N_i}} \partial_t \psi dr dz \right)^\top, $$

with $\Omega_{c_i,j}$ denotes the jth coil section in the ith circuit. The two matrices $S_i \in \mathbb{R}^{(M_i+N_i) \times M_i}$ and $R_i \in \mathbb{R}^{(M_i+N_i) \times N_i}$, which multiply $\vec{V}_i$ and $\vec{\Psi}_i(\partial_t \psi)$, respectively, contain information about the ith circuit, namely the physical characteristics of its components and their relative connections. The current density in (★) for problem (9) flowing in $\Omega_{c_i,l}$ is

$$ (j_{\varphi})_{\ell} = \frac{1}{|\Omega_{c_i, \ell}|} \left( S_i \vec{V}_i(t) + R_i \vec{\Psi}_i(\partial_t \psi) \right)_{\ell}, \quad \ell = 1, \dots, N_i. $$

For problem (9) in the tokamak WEST, the coils $\Omega_{c_i}$ are in $\Omega^{ex}$ (i.e., $L^{in}U^{in}_I = 0$). The weak form of (9) reads: given $T \ge 0$, find the function $\psi: t \in [0, T] \mapsto \psi(t) = (\psi_{ex}(t), \psi_{in}(t)) \in \mathcal{V}$ such that $\psi(0) = \psi_0$ and, $\forall s = (v, w) \in \mathcal{V}_{0,\mathcal{I}}$, it holds

$$ a(\psi(t), s) - j_{ps}(\dot{\psi}(t), s) - j_c(\dot{\psi}(t), s) = \ell^{\text{ex}}(S_i \vec{V}_i(t), s), \qquad (20) $$

where $\dot{\psi}(t)$ is the time derivative of $\psi(t)$ and

$$ \begin{align} a(\psi(t), s) &:= a^{\text{ex}}(\psi_{\text{ex}}(t), v) + a^{\text{in}}(\psi_{\text{in}}(t), w), \nonumber \\ j_{\text{ps}}(\dot{\psi}(t), s) &:= -\sum_{j=1}^{N_{\text{ps}}} \int_{\Omega_{\text{ps},j}} \frac{\sigma_j}{r} \dot{\psi}(t) (\chi_{\Omega^{\text{ex}}} v + \chi_{\Omega^{\text{in}}} w) dr dz, \nonumber \\ j_c(\dot{\psi}(t), s) &:= \sum_{i=1}^{L} \sum_{j=1}^{N_i} \frac{(R_i \vec{\Psi}_i(\dot{\psi}(t)))_j}{|\Omega_{c_i,j}|} \int_{\Omega_{c_i,j}} \chi_{\Omega^{\text{ex}}} v dr dz, \tag{21} \\ \ell^{\text{ex}}(S_i \vec{V}_i(t), s) &:= \sum_{i=1}^{L} \sum_{j=1}^{N_i} \frac{(S_i \vec{V}_i(t))_j}{|\Omega_{c_i,j}|} \int_{\Omega_{c_i,j}} \chi_{\Omega^{\text{ex}}} v dr dz, \nonumber \end{align} $$

with $a^{\text{ex}}(.,.), a^{\text{in}}(.,.),$ the bilinear forms defined in (12). Problem (20) has to be discretized in space and in time.

The semi-discrete evolution problem is obtained by applying to (20), in space, the mortar finite element approach detailed in Section 3.2. The semi-discrete problem thus reads: given $T \ge 0$, find the function $\psi_h: t \in [0, T] \mapsto \psi_h(t) = (\psi_h^{\text{ex}}(t), \psi_h^{\text{in}}(t)) \in \mathcal{V}_h$ such that $\psi_h(0) = \psi_{0,h}$ and, $\forall s_h = (v_h, w_h) \in \mathcal{V}_h^{\text{ex}} \times \mathcal{V}_h^{\text{in}}$, it holds

$$ a(\psi_h, s_h) - j_{ps}(\dot{\psi}_h(t), s_h) - j_c(\dot{\psi}_h(t), s_h) = \ell^{\text{ex}}(S_i \vec{V}_i(t), s_h) \qquad (22) $$
with $\psi_{0,h}$ a representation of $\psi_0$ in the space $V_h$. The bilinear forms appearing in (22) are defined as for (20) and evaluated for functions in the discrete spaces. We finally rely on an implicit Euler scheme to fully discretize in time the semi-discrete problem. The discrete problem thus reads: Find $\psi_h^k = (\psi_h^{ex,k}, \psi_h^{in,k}) \in V$ approximating $\psi_h(t_k)$, for $t_k \in [0, T]$, such that, $\forall s_h = (v_h, w_h) \in V_o^{ex} \times V_o^{in}$, it holds

$$a_{evol}(\psi_h^k, s_h) := a_{evol}^{\text{ex}}(\psi_h^k, v_h) + a_{evol}^{\text{in}}(\psi_h^k, w_h) = \ell^{\text{ex}}(S_i \vec{V}_i^k, s_h) \quad (23)$$

with:

$$\begin{align*}
a_{evol}^{\text{ex}}(\psi_h^k, v) &:= a^{\text{ex}}(\psi_h^k, v) - \frac{1}{\Delta t}(j_{\text{ps}}(\psi_h^k, v) - j_{\text{ps}}(\psi_h^{k-1}, v)) \\
&\phantom{:=} - \frac{1}{\Delta t}(j_c(\psi_h^k, v) - j_c(\psi_h^{k-1}, v)), \\
a_{evol}^{\text{in}}(\psi_h^k, w) &:= a^{\text{in}}(\psi_h^k, w) - \frac{1}{\Delta t}(j_{\text{ps}}(\psi_h^k, w) - j_{\text{ps}}(\psi_h^{k-1}, w)), \\
j_{\text{ps}}(\psi^k, v) &:= -\sum_{j=1}^{N_{\text{ps}}} \int_{\Omega_{c_{i,j}}} \frac{\sigma_j}{r} \psi^k v \, dr \, dz, \tag{24} \\
j_c(\psi^k, v) &:= \sum_{i=1}^{L} \sum_{j=1}^{N_i} \frac{(\mathbf{R}_i, \vec{\psi}_i(\psi^k))_j}{|\Omega_{c_{i,j}}|} \int_{\Omega_{c_{i,j}}} \chi_{\Omega^{\text{ex}}} v \, dr \, dz, \\
\ell^{\text{ex}}(S_i \vec{V}_i^k, s) &:= \sum_{i=1}^{L} \sum_{j=1}^{N_i} \frac{(S_i \vec{V}_i^k)_j}{|\Omega_{c_{i,j}}|} \int_{\Omega_{c_{i,j}}} \chi_{\Omega^{\text{ex}}} v \, dr \, dz.
\end{align*}$$

Hence, for the evolution problem, if we denote $\mathbf{X}_{\text{before}}$ the solution $\mathbf{X}$ computed at the previous time step, equation (22) in its fully discretized form reads

$$\begin{align}
\mathbf{e}(\mathbf{X}) := \mathbf{Q}^{\top} \left[ \begin{aligned}
& (\mathbf{A}_0 + \mathbf{C} - \frac{1}{\Delta t}(\mathbf{M}_{\text{ps}} + \mathbf{LRL}^{\top}))\mathbf{Q}\mathbf{X} \\
& + \frac{1}{\Delta t}(\mathbf{M}_{\text{ps}} + \mathbf{LRL}^{\top})\mathbf{Q}\mathbf{X}_{\text{before}} \\
& + \mathbf{A}_{\mu}(\psi) - \mathbf{J}(\psi) - \mathbf{LSV}
\end{aligned} \right], \tag{25}
\end{align}$$

where $\mathbf{L}$ stands for $\mathbf{L}^{\text{ex}}$, and $\mathbf{M}_{\text{ps}}$ (resp. $\mathbf{LRL}^{\top}$) is the matrix associated with the bilinear form $j_{\text{ps}}(\cdot, \cdot)$ (resp. $j_c(\cdot, \cdot)$). Therefore, we have

$$[\mathbf{e}_{\mathbf{X}}(\mathbf{X}^k)] = \mathbf{Q}^{\top}[\mathbf{A}_0 + \mathbf{C} - \frac{1}{\Delta t}(\mathbf{M}_{\text{ps}} + \mathbf{LRL}^{\top}) + \mathbf{A}_{\mu, \psi}(\psi^k) - \mathbf{Jac}_{\psi}(\psi^k)]\mathbf{Q}. \quad (26)$$

It does not exist, to our knowledge, analytical solutions for the free-boundary equilibrium problems considered in this paper. We provide nevertheless some numerical evidence of convergence for the proposed method.

## 5 Numerical results

The initial guess of the plasma domain $\Omega_p(\psi)$ for given currents in the poloidal field coils plays a crucial role in free-boundary equilibrium problems. Here, we find such initial guesses by solving 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. 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., $\mathbf{X}^0$ could be a non-mortar formulation of the free-boundary equilibrium problem involving piece-wise linear FEs everywhere, as explained in [11]).
<table><caption>**Table 1** Convergence history of Newton iterations for WEST: iteration number *n* and residual relative error ||**X**<sup>*n*</sup> − **X**<sup>*n*−1</sup>||/||**X**<sup>*n*−1</sup>||, with either Ω<sub>Fe</sub> = ∅ or Ω<sub>Fe</sub> ≠ ∅.</caption><thead><tr><th>n</th><th>if Ω<sub>Fe</sub> = ∅</th><th>if Ω<sub>Fe</sub> ≠ ∅</th></tr></thead><tbody><tr><td>1</td><td>1.77919 × 10<sup>-2</sup></td><td>9.78267 × 10<sup>-3</sup></td></tr><tr><td>2</td><td>4.35470 × 10<sup>-4</sup></td><td>8.56241 × 10<sup>-4</sup></td></tr><tr><td>3</td><td>4.05152 × 10<sup>-7</sup></td><td>6.17721 × 10<sup>-4</sup></td></tr><tr><td>4</td><td>2.80505 × 10<sup>-11</sup></td><td>3.46301 × 10<sup>-4</sup></td></tr><tr><td>5</td><td></td><td>9.38432 × 10<sup>-5</sup></td></tr><tr><td>6</td><td></td><td>3.09165 × 10<sup>-5</sup></td></tr><tr><td>7</td><td></td><td>8.76154 × 10<sup>-6</sup></td></tr><tr><td>8</td><td></td><td>3.11734 × 10<sup>-6</sup></td></tr><tr><td>9</td><td></td><td>5.16426 × 10<sup>-7</sup></td></tr><tr><td>10</td><td></td><td>2.58976 × 10<sup>-10</sup></td></tr><tr><td>11</td><td></td><td>6.29602 × 10<sup>-14</sup></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Poloidal cross section showing ψ-isolines for a WEST typical equilibrium, supposing either Ω<sub>Fe</sub> = ∅ (left) or Ω<sub>Fe</sub> ≠ ∅ (right). The red curve is the plasma boundary and the blue one is the coupling interface.</figcaption></figure>

We start with the static problem. The good convergence of the Newton iterations applied to (17), with error threshold fixed to 10<sup>-10</sup>, is shown in Table 2. We remark that the presence of iron parts slows down considerably the convergence speed. Fig. 4 shows a typical WEST poloidal flux map calculated by NICE. A zoom on the distribution of ψ<sub>h</sub> in Ω<sup>in</sup> is proposed in Fig. 5, for the case Ω<sub>Fe</sub> ≠ ∅. The X-point and plasma axis are highlighted in the small pictures. With the adopted cubic C<sup>1</sup> FEs in Ω<sup>in</sup>, these points do not coincide with nodes of the computational mesh, thus assuming a more physically meaningful position than the one that could be computed with low-order C<sup>0</sup> FEs.

For the quasi-static evolution problem, simulations are done with Δt = 0.001 s over the time interval [0, 0.03]. In Table 2, we report the relative residual for the Newton’s algorithm at three time steps of the simulation with Ω<sub>Fe</sub> ≠ ∅ (the error threshold is still fixed to 10<sup>-10</sup>). In the quasi-static evolution case the
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Left: Zoom on the inner vacuum vessel region inside the coupling interface (the $\psi$-isolines in $\Omega^{in}$ cross without discontinuity the interface). Right: Zoom around the plasma axis (red point, top) and X-point (red point, bottom).</figcaption></figure>

<table><caption>**Table 2** Convergence history of Newton iterations for WEST: iteration number $n$ and residual relative error $||X^n - X^{n-1}||/||X^{n-1}||$, for the times $t = 0.01$ s, $t = 0.02$ s, $t = 0.03$ s.</caption><thead><tr><th>n</th><th>t = 0.01</th><th>if t = 0.02</th><th>if t = 0.03</th></tr></thead><tbody><tr><td>1</td><td>8.32182 × 10<sup>-3</sup></td><td>1.11224 × 10<sup>-2</sup></td><td>1.39098 × 10<sup>-2</sup></td></tr><tr><td>2</td><td>2.37455 × 10<sup>-4</sup></td><td>3.01691 × 10<sup>-4</sup></td><td>4.68653 × 10<sup>-4</sup></td></tr><tr><td>3</td><td>8.37148 × 10<sup>-7</sup></td><td>1.27151 × 10<sup>-6</sup></td><td>4.72572 × 10<sup>-6</sup></td></tr><tr><td>4</td><td>9.49314 × 10<sup>-12</sup></td><td>3.72105 × 10<sup>-11</sup></td><td>8.28915 × 10<sup>-10</sup></td></tr><tr><td>5</td><td></td><td></td><td>2.57022 × 10<sup>-14</sup></td></tr></tbody></table>

convergence of the Newton's method is faster than for the static problem since, at each time step, we start close to an equilibrium. Figures 6-9, present the separatrix computed by a full P1 approach (blue line) and by a coupled P1-HCT one (red line), at four time steps of the simulation over a given mesh. In all the cases, the blue and the red lines are very close one from another, stating that the mortar coupling in the quasi-static evolution case is correctly realised. The coarser mesh zones visible in the figures are parts of the limiter. In Figures 7-9 the volume bounded by the separatrix looks smaller in the C0-C1 (red) case because the extremum value of $\psi$ on the limiter is taken at a point of its boundary, not necessarily coincident with a point of the mesh.

Finally, the computed behavior of the plasma looks more realistic when a C0-C1 approach is adopted. Indeed, as shown in Figure 10, the C0-C1 simulation yields a displacement of the magnetic axis without the staircase effects that appear with the full C0 simulation.
<figure><img src="image_7.png" /><figcaption>Fig. 6 Comparison between P1 and P1-HCT plasma boundary computations at $t = 0.001$ s.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 Comparison between P1 and P1-HCT plasma boundary computations at $t = 0.01$ s.</figcaption></figure>

## 6 Conclusions

We have focused on the numerical computation of a MHD equilibrium for a hot plasma in iron tokamaks, such as WEST or JET, still in activity nowadays. A short
<figure><img src="image_9.png" /><figcaption>Fig. 8 Comparison between P1 and P1-HCT plasma boundary computations at $t = 0.02$ s.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9 Comparison between P1 and P1-HCT plasma boundary computations at $t = 0.03$ s.</figcaption></figure>

An overview on the mathematical complexity to treat the magnetic induction **B** in tokamak plasmas has preceded the equations. We have underlined that axisymmetric plasma equilibrium simulations need to rely on accurate reconstructions of the
<figure><img src="image_11.png" /><figcaption>**Fig. 10** Time evolution of the magnetic axis coordinates $r_a$ (left) $z_a$ (right): comparison of the results obtained by a full C0 (P1, in blue) approach and a coupled C0-C1 (P1CT, in red) approach.</figcaption></figure>

poloidal magnetic flux $\psi$ and of its gradient $\nabla\psi$, at least in the part of the tokamak cross section that is accessible to the plasma. Therefore we have considered the $C^1$, piece-wise cubic, reduced Hsieh-Clough-Tocher (rHCT) finite elements with $C^0$ piece-wise linear Lagrange ones by relying on a non-overlapping mortar element method and to discretize, in such a cross section, the Grad-Shafranov-Schlüter's equation. At the discrete level a Newton's method is proposed to solve the coupled nonlinear problem and numerical evidence of its convergence is given. The presence of material non-linearities has been taken into account in the algorithm by suitably modifying the computation of the Jacobian matrix.

## Acknowledgements

The authors are grateful to Jacques Blum for his precious enlightenments on plasma physics and to Holger Heumann for stimulating conversations on the derivation and approximation of the Grad-Shafranov-Schlüter's equation.
FR thanks INRIA for the delegation during which this work was accomplished.

## Declarations

The authors declare no competing interests.

GG is a PhD student from the Université Côte d'Azur.

FR received funding from the French National Research Agency under Grant Agreement SISTEM (ANR-19-CE46-0005-03).

CB and BF have worked within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 — EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.
References

1. Albanese, R., Blum, J., Barbieri, O.: On the solution of the magnetic flux equation in an infinite domain. In: EPS. 8th Europhysics Conference on Computing in Plasma Physics (1986), pp. 41–44 (1986)

2. Bernardi, C., Maday, Y., Patera, A.: A new nonconforming approach to domain decomposition: The mortar element method. In: Nonlinear Partial Differential Equations and Their Applications. in H. Brézis and J.-L. Lions (eds.), Collège de France Seminar XI. (1992)

3. Blackman, E.G.: Magnetic helicity and large scale magnetic fields: a primer. J. Fluid Mechanics **188**, 59–91 (2015)

4. Blum, J.: Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks. Series in Modern Applied Mathematics. Wiley Gauthier-Villars, Paris (1989)

5. Blum, J., Heumann, H., Nardon, E., Song, X.: Automating the design of tokamak experiment scenarios. J. Comp. Phys. **394**, 594–614 (2019)

6. Blum, J., Le Foll, J., Thooris, B.: The self-consistent equilibrium and diffusion code SCED **24**, 235–254 (1981)

7. Cantarella, J., DeTurck, D., Gluck, H., Teytel, M.: Influence of geometry and topology on helicity. Geophysical Monograph-American Geophysical Union **111**, 17–24 (1999)

8. Christiansen, S.H., Hu, K.: Generalized finite element systems for smooth differential forms and Stokes’ problem. Numer. Math. **140**, 327–371 (2018)

9. Clough, R., Tocher, J.: Finite element stiffness matrices for analysis of plates in bending. In: Proc. Conf. Matrix Methods in Struct. Mech. Air Force Inst of Tech., Wright Patterson A.F Base, Ohio (1965)

10. Elarif, A., Faugeras, B., Rapetti, F.: Tokamak free-boundary plasma equilibrium computation using finite elements of class $C^0$ and $C^1$ within a mortar element approach. J. Comput. Phys. **439**, 110388 (2021)

11. Faugeras, B.: An overview of the numerical methods for tokamak plasma equilibrium computation implemented in the NICE code. Fusion Eng. Design **160**, 112020 (2020)

12. Glowinski, R., Marrocco, A.: Analyse numérique du champ magnétique d'un alternateur par éléments finis et sur-relaxation ponctuelle non linéaire. Computer Meth. Appl. Mech. Engng. **3**, 55–85 (1974)

13. Glowinski, R., Marrocco, A.: Sur l'approximation par éléments finis d'ordre 1, et la résolution par pénalisation-dualité, d'une classe de problèmes de Dirichlet non linéaires. C.R.A.S. Serie A, Paris **278**, 1649–1652 (1974)

14. Grad, H., Rubin, H.: Hydromagnetic equilibria and force-free fields. In: 2nd UN Conf. on the Peaceful Uses of Atomic Energy, vol. 31, p. 190 (1958)

15. Heumann, H.: A galerkin method for the weak formulation of current diffusion and force balance in tokamak plasmas. J. Comp. Phys. **442**, 110483 (2021)

16. Heumann, H., Blum, J., Boulbe, C., Faugeras, B., Selig, G., Ané, J.M., Brémond, S., Grandgirard, V., Hertout, P., Nardon, E.: Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications. J. Plasma Physics **81**, 905810301 (2015)

17. Heumann, H., Rapetti, F.: A finite element method with overlapping meshes for free-boundary axisymmetric plasma equilibria in realistic geometries. J. Comput. Phys. **334**, 522–540 (2017)

18. Jardin, S.: Computational Methods in Plasma Physics. Chapman & Hall/CRC Computational Science. CRC Press, Boca Raton (2010)

19. Lüst, R., Schlüter, A.: Axialsymmetrische magnetohydrodynamische gleichgewichtskonfigurationen. Z. Naturforsch. A **12**, 850–854 (1957)

20. MacTaggart, D., Valli, A.: Magnetic helicity in multiply connected domains. J. Plasma Physics **85**, 775850501 (2019)

21. Minjeaud, S., Pasquetti, R.: Fourier-spectral elements approximation of the two fluid ion-electron braginskii system with application to tokamak edge plasma in divertor configuration. J. Comput. Phys. **321**, 492–511 (2016)

22. Moffatt, H., Ricca, R.: Helicity and the calugareanu invariant. pp. 411–429 (1992)

23. Moffatt, H.K.: The degree of knottedness of tangled vortex lines. J. Fluid Mechanics **35**, 117–129 (1969)

24. Ratnani, A., Crouseilles, N., Sonnendrücker, E.: An isogeometric analysis approach for the study of the gyrokinetic quasi-neutrality equation. J. Comput. Phys. **231**, 373–393 (2012)
25. Rebut, P.: Instabilités non magnétohydrodynamiques dans les plasmas à densités de courant élevé. J. Nuclear Energy Part C 4, 159 (1963)

26. Shafranov, V.: On magnetohydrodynamical equilibrium configurations. Soviet Journal of Experimental and Theoretical Physics 6, 545 (1958)