# Numerical simulation of tokamak plasma equilibrium evolution

Guillaume Gros, Blaise Faugeras, C. Boulbe, J.F. Artaud, Rémy Nouailletas, Francesca Rapetti

## ► To cite this version:

Guillaume Gros, Blaise Faugeras, C. Boulbe, J.F. Artaud, Rémy Nouailletas, et al., Numerical simulation of tokamak plasma equilibrium evolution. RR-9548, INRIA. 2024. hal-04589897

HAL Id: hal-04589897

https://inria.hal.science/hal-04589897

Submitted on 28 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
Numerical simulation of tokamak plasma equilibrium evolution

Guillaume Gros, Blaise Faugeras, C. Boulbe, J.F. Artaud, Rémy Nouailletas,

Francesca Rapetti

► To cite this version:

Guillaume Gros, Blaise Faugeras, C. Boulbe, J.F. Artaud, Rémy Nouailletas, et al. Numerical sim- ulation of tokamak plasma equilibrium evolution. Journal of Computational Physics, 2025, 529, pp.113849. 10.1016/j.jcp.2025.113849 . hal-04589897v2

HAL Id: hal-04589897

https://inria.hal.science/hal-04589897v2

Submitted on 27 Feb 2025

**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
informatics mathematics
Inria

Numerical simulation of tokamak plasma
equilibrium evolution

G. Gros, B. Faugeras, C. Boulbe, J.-F. Artaud, R. Nouailletas, F. Rapetti

RESEARCH
REPORT

Nº 9548

May 2024

Project-Team CASTOR

ISSN 0249-6399 ISRN INRIA/RR-9548-FR+ENG
Numerical simulation of tokamak plasma equilibrium evolution

G. Gros<sup>a,*</sup>, B. Faugeras<sup>a</sup>, C. Boulbe<sup>a</sup>, J.-F. Artaud<sup>b</sup>, R. Nouailletas<sup>b</sup>, F. Rapetti<sup>a</sup>

<sup>a</sup>Université Côte d'Azur, CNRS, INRIA Castor, LJAD, F-06108 Nice, France

<sup>b</sup>CEA, IRFM, F-13108 Saint Paul-lez-Durance, France

**Abstract**

This paper focuses on the numerical methods recently developed in order to simulate the time evolution of a tokamak plasma equilibrium at the resistive diffusion time scale. Starting from the method proposed by Heumann in 2021 for the coupling of magnetic equilibrium and current diffusion, we introduce a new space discretization for the poloidal flux using coupled $C^0$ and $C^1$ finite elements. This, together with the use of cubic spline functions to represent the poloidal current function in the resistive diffusion equation, enables to restrain numerical oscillations which can occur with the original method. In order to compute consistently the plasma resistivity and the non-inductive bootstrap current terms needed in the resistive diffusion equation we add to the model an evolution equation for electron temperature in the plasma. It is also used to evolve the pressure term in the simulation. These numerical methods are implemented in the plasma equilibrium code NICE. A free plasma displacement is simulated and comparison with experimental results from the WEST tokamak are used to validate the simulation. The code is also coupled to a magnetic feedback controller making it possible to simulate a prescribed plasma scenario. The results for an X-point formation scenario in the WEST tokamak are presented as an illustration of the efficiency of the developed numerical methods.

**Keywords:** tokamak, plasma, equilibrium, resistive diffusion, heat transport, FEM, control, scenario simulation

<footnote>*Corresponding author</footnote>

<footnote>Email address: guillaume.gros@univ-cotedazur.fr (G. Gros)</footnote>
<figure><img src="image_2.png" /></figure>
# 1. Introduction

2

4

6

8

10

12

14

16

18

20

22

24

26

28

30

32

34

Numerical simulation of tokamak plasma equilibrium evolution at the timescale of current and heat diffusion also known as resistive diffusion timescale, is essential to analyse existing Tokamaks discharges, to prepare new plasma scenarios and to develop and test plasma controllers [1]. The goal of this paper is to present the model and numerical methods recently developed in order to achieve such simulations. We refer to text books [2, 3] for the basics of tokamak plasma modelling. In particular the notion of average over magnetic surfaces which is used in this work can be found in these references.

Since the first studies [4, 5, 6] introducing this technique of averaging over magnetic surfaces allowing to write a set of 1D averaged transport equations simulating the diffusion of the plasma pressure and of the toroidal magnetic field inside the plasma, several papers have been published which describe the numerical simulation, and the associated codes, of the coupling between this set of 1D transport equations and the 2D Grad-Shafranov equilibrium equation [7, G2M], [8, SCED], [9, 10], [11, TSC], [12, CORSICA], [13, DINA-CHCRONOS], [14, TOPICS], [15, ASTRA-SPIDER], [16, JINTRAC-CREATENL], [17], [18, FEEQS-METIS]. From a high perspective such a coupling goes as follows. Given the pressure profile $p$ and the poloidal current function $f$ in the plasma as well as the voltages in poloidal field coils circuits, the 2D equilibrium code evolves the poloidal flux $\psi(r, z, t)$ and the equilibrium configuration. On the other hand, given an equilibrium configuration the 1D transport code evolves the pressure profile and the poloidal current function. This is commonly done by including in the transport system the so-called resistive diffusion equation, a 1D evolution equation for $\psi(\rho, t)$ where $\rho$ is a label for magnetic surfaces. Using the average Grad-Shafranov equation one can deduce the evolution of $f$ from the evolution of $\psi(\rho, t)$. Among the numerous numerical difficulties arising with such a strategy the most challenging is certainly to ensure the consistency between the two representations of the same quantity, $\psi(r, z, t)$ and $\psi(\rho, t)$ which evolve with 2 different equations. Different techniques and algorithms to tackle this issue can be found in the previous references. It is necessary to iterate between the equilibrium equation and the resistive diffusion equation at each time step. Moreover a good choice for the boundary condition to be applied at the plasma boundary in the resistive diffusion equation is not very clear (see for e.g. [17, Section
<figure><img src="image_3.png" /></figure>

# Numerical simulation of tokamak plasma
equilibrium evolution

G. Gros*, B. Faugeras*, C. Boulbe*, J.-F. Artaud†, R. Nouailletas†, F. Rapetti*

Project-Team CASTOR

Research Report n° 9548 — May 2024 — 38 pages

**Abstract:** This paper focuses on the numerical methods recently developed in order to simulate the time evolution of a tokamak plasma equilibrium at the resistive diffusion time scale. Starting from the method proposed by Heumann in 2021 for the coupling of magnetic equilibrium and current diffusion, we introduce a new space discretization for the poloidal flux using coupled $C^0$ and $C^1$ finite elements. This, together with the use of cubic spline functions to represent the poloidal current function in the resistive diffusion equation, enables to restrain numerical oscillations which can occur with the original method. In order to compute consistently the plasma resistivity and the non-inductive bootstrap current terms needed in the resistive diffusion equation we add to the model an evolution equation for electron temperature in the plasma. It is also used to evolve the pressure term in the simulation. These numerical methods are implemented in the plasma equilibrium code NICE. A vertical displacement event is simulated and comparison with experimental results from the WEST tokamak are used to validate the simulation. The code is also coupled to a magnetic feedback controller making it possible to simulate a prescribed plasma scenario. The results for an X-point formation scenario in the WEST tokamak are presented as an illustration of the efficiency of the developed numerical methods.

**Key-words:** tokamak, plasma, equilibrium, resistive diffusion, heat transport, FEM, control, scenario simulation

RESEARCH CENTRE
SOPHIA ANTIPOLIS – MÉDITERRANÉE

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex

<footnote>* Université Côte d'Azur, CNRS, INRIA CASTOR, LJAD</footnote>

<footnote>† IRFM, CEA, Cadarache</footnote>
36

4.3]) and such algorithms are known not to be very robust.

38

40

42

44

46

In a recent paper [19] the author proposes an alternative approach for the simulation of plasma equilibrium evolution by solving together with the equilibrium equation an equation for the evolution of the poloidal current function $f$. The idea of directly solving an equation for $f$ can already be found e.g. in [20, 7] but using a surface averaged 1D equation which makes it very difficult to be solved simultaneously with the 2D equilibrium equation. In [19] the equation for $f$ is still essentially 1D since the unknown is of the form $f(\psi(r, z, t), t)$ but magnetic surface averaging is avoided using a weak formulation. This enables to solve simultaneously and consistently the two equations for $\psi(r, z, t)$ and $f(\psi(r, z, t), t)$ as a system of two non-linear equations.

48

50

52

54

56

58

This contribution develops on the method proposed in [19] for the coupling of magnetic equilibrium and current diffusion. The original implementation of [19] is based on a classical $C^0$ piecewise linear P1 finite element method (FEM) for the discretization of the poloidal flux $\psi(r, z, t)$ and Legendre polynomials up to degree 10 for the discretization of the poloidal current function $f$. We introduce a new space discretization for the poloidal flux using coupled $C^0$ and $C^1$ FEM [21, 22]. This, together with the use of cubic spline functions to represent $f$, enables to restrain numerical oscillations which can occur with the original method (see Appendix A). Moreover the use of a high order $C^1$ FEM to represent the plasma enables a smooth evolution of the plasma boundary and magnetic axis which is not possible with P1 FEM.

60

62

64

We add to the model an evolution equation for electron temperature in the plasma. This enables us to evolve the pressure term in the simulations and to compute consistently the plasma resistivity and the non-inductive current terms called bootstrap current needed in the resistive diffusion equation for $f$.

66

68

70

72

The model is presented in Section 2 and the discretization methods in Section 3. These methods are implemented in the code NICE [23] and numerical results are provided in Section 4. A free plasma displacement simulation for the WEST tokamak [24, 25, 26] is conducted and comparison with experimental results from shot 59695 are used to validate the simulation. This section ends with results for a simulation with zero plasma resistivity to further validate the proposed approach. Eventually the code is coupled to a magnetic feedback controller [27] enabling the simulation of a prescribed plasma scenario. The results for an X-point formation scenario in the WEST
**Résumé** : Cet article se concentre sur les méthodes numériques récemment développées dans le but de simuler l'évolution temporelle de l'équilibre du plasma dans un tokamak à l'échelle de temps de la diffusion résistive. Avec comme point de départ la formulation proposée par Heumann en 2021 pour le couplage de l'équilibre magnétique avec la diffusion du courant, nous introduisons une nouvelle discrétisation spatiale pour le flux poloidal utilisant un couplage d'élément finis de classes $C^0$ et $C^1$. Ceci ajouté à l'utilisation de splines cubiques pour approcher la fonction de courant poloidal dans l'équation de diffusion résisitive nous permet de restreindre les oscillations numériques qui peuvent apparaître dans la méthode originelle. Afin de calculer de manière consistente la résistivité du plasma et les courants non-inductifs dits de bootstrap qui sont nécessaires à la résolution de l'équation de diffusion résistive, nous ajoutons au modèle une équation d'évolution pour la température électronique dans le plasma. Ceci permet également de faire évoluer le terme de pression dans la simulation. Ces méthodes numériques sont implémentées dans le code d'équilibre NICE. Une instabilité verticale est simulée et comparée à des résultats expérimentaux d'une décharge du tokamak WEST. Le code est également couplé à un contrôleur magnétique rendant possible la simulation d'un scenario prescrit. Les résultats pour un scenario de formation du point X dans le tokamak WEST sont présentés et illustrent l'efficacité des méthodes numériques développées dans ce travail.

**Mots-clés** : tokamak, plasma, équilibre, diffusion résistive, transport de la chaleur, éléments finis, contrôle, simulation de scenario
tokamak are presented as an illustration of the efficiency of the developed numerical methods. Results for a vertical displacement event (VDE) are also obtained turning off the vertical position magnetic control.

76

## 2. Continuous model

The model for plasma equilibrium evolution presented here is made up of three main state equations: the 2D quasi-static equilibrium evolution equation for the poloidal flux $\psi$ described and used in [10, 2, 28, 29, 19, 22], the resistive current diffusion equation given as an evolution equation for the poloidal current function $f(\psi, t)$ provided in [19] and an equation for the electron energy transport inside the plasma. The knowledge of the electron temperature $T_e$ enables us to compute consistently the pressure term needed in the equilibrium equation as well as the parallel resistivity $\eta^{\parallel}$ and the bootstrap current which appear in the evolution equation for $f$.

78

80

82

84

86

### 2.1. Equilibrium

88

90

92

Let us introduce a cylindrical $(r, \phi, z)$ coordinate system and assume axisymmetry. Maxwell's equations together with force balance in the plasma reduce to the following equation for the poloidal flux $\psi(r, z)$ in the poloidal section of the tokamak contained in the positive half plane $D = [0, \infty] \times [-\infty, \infty]$:

$$ \mathcal{L}\psi = j_{\phi}, \qquad (1) $$

94

where $j_{\phi}$ is the toroidal component of the current density $\mathbf{j}$, and the second order elliptic differential operator $\mathcal{L}$ is defined by

$$ \mathcal{L}\psi := -\nabla \cdot \left( \frac{1}{\mu(\psi)r} \nabla\psi \right), \qquad (2) $$

96

with $\nabla$ the 2D operator in the $(r, z)$-plane and $\mu(\psi)$ is the magnetic permeability. The magnetic field decomposes in a poloidal and toroidal component $\mathbf{B} = \mathbf{B}_p + \mathbf{B}_{\phi}$ with $\mathbf{B}_p = \frac{1}{r}(-\frac{\partial\psi}{\partial z}\mathbf{e}_r + \frac{\partial\psi}{\partial r}\mathbf{e}_z)$ and $\mathbf{B}_{\phi} = \frac{f}{r}e_{\phi}$.

98

The geometry of the tokamak (see Fig. 2) determines various sub-domains that are used to set the expression of $j_{\phi}$ accordingly.

100

- $\Omega_{Fe} \subset D$ denotes those parts of $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) = \mu_{Fe}(|\nabla\psi|^2 r^{-2}) \ge \mu_0$; elsewhere, that is in $\Omega_{Fe}^c$, $\mu = \mu_0$ the permeability of vacuum;

102
# Contents

<table><tr><td>1</td><td>Introduction</td><td>4</td></tr><tr><td>2</td><td>Continuous model</td><td>5</td></tr><tr><td>2.1</td><td>Equilibrium</td><td>5</td></tr><tr><td>2.2</td><td>Resistive current diffusion</td><td>6</td></tr><tr><td>2.3</td><td>Electron energy transport, resistivity and bootstrap current</td><td>7</td></tr><tr><td>3</td><td>Discrete problem</td><td>10</td></tr><tr><td>3.1</td><td>The equilibrium (ψ, f) system</td><td>10</td></tr><tr><td>3.2</td><td>Matrix problem and Newton method for the (ψ, f) system</td><td>13</td></tr><tr><td>3.3</td><td>The T<sub>e</sub> equation</td><td>14</td></tr><tr><td>4</td><td>Numerical results</td><td>15</td></tr><tr><td>4.1</td><td>Simulation setup</td><td>15</td></tr><tr><td>4.2</td><td>Free displacement of the plasma. Test case 1 and 2.</td><td>19</td></tr><tr><td>4.3</td><td>Controlled evolution</td><td>25</td></tr><tr><td>5</td><td>Conclusion</td><td>28</td></tr><tr><td>A</td><td>Test of the original numerical method from [25]</td><td>29</td></tr><tr><td>B</td><td>Sources, resistivity and bootstrap current</td><td>31</td></tr><tr><td>C</td><td>Derivatives computation</td><td>32</td></tr></table>
104 - $\Omega_{c_i} \subset \mathcal{D}$, $1 \le i \le N_c$, denotes the intersection with the poloidal plane of the $i$th poloidal field coil carrying currents and $N_c$ denotes the number of coils. The $i$th coil has cross section area $|\Omega_{c_i}|$ and carries a total current $I_i$;

108 - $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}$, characterized by an electric conductivity $\sigma_j \ne 0$, $j = 1, ..., N_{ps}$ (note that the tokamak vacuum vessel is a passive structure too);

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

114 - $\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 Grad-Shafranov equation [30, 31] is satisfied in this domain.

The equilibrium of the plasma in a tokamak has to satisfy at each instant the following non-linear initial boundary value problem: starting from an initial condition $\psi_0(r, z)$ find $\psi(r, z, t)$ such that

$$ \mathcal{L}\psi = j_{\phi} = \begin{cases} rp'(\psi, t) + \frac{1}{\mu_0 r} ff'(\psi, t) & \text{in } \Omega_p(\psi), \\ \frac{1}{|\Omega_{c_i}|} I_i\left(\frac{d\psi}{dt}, t\right) & \text{in } \Omega_{c_i}, i = 1, \dots, N_c, \\ -\frac{\sigma_j}{r} \frac{d\psi}{dt} & \text{in } \Omega_{ps_j} \subset S, j = 1, \dots, N_{ps}, \\ 0 & \text{elsewhere,} \end{cases} \\ I_i\left(\frac{d\psi}{dt}, t\right) = \sum_{j=1}^{N_s} S_{ij}V_j(t) + \sum_{j=1}^{N_c} R_{ij}\frac{1}{|\Omega_{c_j}|}\int_{\Omega_{c_j}}\frac{d\psi}{dt}drdz, \quad i = 1, \dots, N_c \tag{3} $$

with boundary conditions $\psi(0, z, t) = 0$ and $\lim_{||(r,z)||_2 \to +\infty} \psi(r, z, t) = 0$. These boundary conditions are taken into account using the boundary integral method presented in [32] and we have used the exact same implementation as in [28, 33, 23].

The quantity $\frac{d\psi}{dt}$ is the time derivative of $\psi$ at fixed $(r, z)$. The current $I_i$ relates to voltages $V_j(t)$ in the power supplies and to self and mutual induction via electric circuit equations (see [28] Appendix A and [29] Appendix B). $N_s$ is the number of power supplies, matrix $S$ has dimension $N_c \times N_s$ and matrix $R$ has dimension $N_c \times N_s$. The plasma domain $\Omega_p(\psi)$ is unknown and
# 1 Introduction

Numerical simulation of tokamak plasma equilibrium evolution at the timescale of current and heat diffusion also known as resistive diffusion timescale, is essential to analyse existing Tokamaks discharges, to prepare new plasma scenarios and to develop and test plasma controllers [34]. The goal of this paper is to present the model and numerical methods recently developed in order to achieve such simulations. We refer to text books [4, 31] for the basics of tokamak plasma modelling. In particular the notion of average over magnetic surfaces which is used in this work can be found in these references.

Since the first studies [22, 23, 28] introducing this technique of averaging over magnetic surfaces allowing to write a set of 1D averaged transport equations simulating the diffusion of the plasma pressure and of the toroidal magnetic field inside the plasma, several papers have been published which describe the numerical simulation, and the associated codes, of the coupling between this set of 1D transport equations and the 2D Grad-Shafranov equilibrium equation [12, G2M], [7, SCED], [48, 6], [30, TSC], [15, CORSICA], [36, DINA-CH-CRONOS], [29, TOPICS], [18, ASTRA-SPIDER], [40, JINTRAC-CREATE-NL], [19], [39, FEEQS-METIS]. From a high perspective such a coupling goes as follows. Given the pressure profile $p$ and the poloidal current function $f$ in the plasma as well as the voltages in poloidal field coils circuits, the 2D equilibrium code evolves the poloidal flux $\psi(r, z, t)$ and the equilibrium configuration. On the other hand, given an equilibrium configuration the 1D transport code evolves the pressure profile and the poloidal current function. This is commonly done by including in the transport system the so-called resistive diffusion equation, a 1D evolution equation for $\psi(\rho, t)$ where $\rho$ is a label for magnetic surfaces. Using the average Grad-Shafranov equation one can deduce the evolution of $f$ from the evolution of $\psi(\rho, t)$. Among the numerous numerical difficulties arising with such a strategy the most challenging is certainly to ensure the consistency between the two representations of the same quantity, $\psi(r, z, t)$ and $\psi(\rho, t)$ which evolve with 2 different equations. Different techniques and algorithms to tackle this issue can be found in the previous references. It is necessary to iterate between the equilibrium equation and the resistive diffusion equation at each time step. Moreover a good choice for the boundary condition to be applied at the plasma boundary in the resistive diffusion equation is not very clear (see for e.g. [19, Section 4.3]) and such algorithms are known not to be very robust.

In a recent paper [25] the author proposes an alternative approach for the simulation of plasma equilibrium evolution by solving together with the equilibrium equation an equation for the evolution of the poloidal current function $f$. The idea of directly solving an equation for $f$ can already be found e.g. in [47, 12] but using a surface averaged 1D equation which makes it very difficult to be solved simultaneously with the 2D equilibrium equation. In [25] the equation for $f$ is still essentially 1D since the unknown is of the form $f(\psi(r, z, t), t)$ but magnetic surface averaging is avoided using a weak formulation. This enables to solve simultaneously and consistently the 2 equations for $\psi(r, z, t)$ and $f(\psi(r, z, t), t)$ as a system of two non-linear equations.

This contribution develops on the method proposed in [25] for the coupling of magnetic equilibrium and current diffusion. The original implementation of [25] is based on a classical $C^0$ piecewise linear P1 finite element method (FEM) for the discretization of the poloidal flux $\psi(r, z, t)$ and Legendre polynomials up to degree 10 for the discretization of the poloidal current function $f$. We introduce a new space discretization for the poloidal flux using coupled $C^0$ and $C^1$ FEM [17, 8]. This, together with the use of cubic spline functions to represent $f$, enables to restrain numerical oscillations which can occur with the original method (see A) Moreover the use of a high order $C^1$ FEM to represent the plasma enables a smooth evolution of the plasma boundary and magnetic axis which is not possible with P1 FEM.
depends non-linearly on the poloidal flux $\psi$ (it is a free-boundary problem). The boundary of $\Omega_p(\psi)$ is an iso-contour of $\psi$ and is either tangent to the boundary of $\Omega_L$ (limiter configuration) or contains one or more saddle points (X-points) of $\psi$ (divertor configuration). The value of $\psi$ at the magnetic axis is denoted $\psi_a$ and on the plasma boundary by $\psi_b$ (the subscript $b$ is used for values at the plasma boundary).

The two functions $p'$ and $f f'$ (where the superscript $'$ denotes here derivation with respect to $\psi$) are the usual terms in the Grad-Shafranov equation. The resistive diffusion equation presented in the next Section provides a way to compute $f$. The pressure term $p$ is computed from the electron temperature as explained further below in Section 2.3.

## 2.2. Resistive current diffusion

In the above equilibrium equation the Grad-Shafranov pressure term $p'$ and diamagnetic term $f f'$ are not known and have to be provided. The pressure term does not lead to significant numerical difficulties. It is not the case for the $f f'$ term. In the literature this term is usually computed from the 1D flux surface average $\psi$ evolution equation, classically called the resistive current diffusion equation, together with the flux surface average Grad-Shafranov equation. With this approach the poloidal flux has two representations, a 2D representation $\psi(r, z)$ from the equilibrium model and a 1D representation $\psi(\rho)$ where $\rho$ is a label for the flux surfaces. This resistive current diffusion equation is usually part of a transport model which is evolved in parallel to the equilibrium model. With such a procedure it is very difficult to ensure numerically the consistency of these two representations of the poloidal flux and to provide a robust simulator. This is the main reason why the idea of using directly an evolution equation for $f$ is attractive. A flux surface average evolution equation for $f$ can be found for example in [20, eqn.(20)] [7, eqn.(24)] but because of the averaging process and of the so-called geometric coefficients which appear it is difficult to solve the equilibrium equation together with the $f$ equation as a single system of equations using Newton's method to resolve the non-linearity. The formulation proposed in [19] solves this issue. From Ohm's law and Faraday's law [19] derives an evolution equation for $f$ given in weak form using 1D test function $g(\psi)$ but integrating over the 2D plasma domain. With this approach no geometric coefficients are needed and the system of equations can be solved in $(\psi, f)$ directly. The evolution problem for $f$ reads as follows: starting from an initial condition $f_0(\psi_0)$ find $f(\psi, t)$ such that
We add to the model an evolution equation for electron temperature in the plasma. This enables us to evolve the pressure term in the simulations and to compute consistently the plasma resistivity and the non-inductive currents terms called bootstrap current needed in the resitive diffusion equation for f.

The model is presented in Section 2 and the discretization methods in Section 3. These methods are implemented in the code NICE [20] and numerical results are provided in Section 4. A vertical displacement event (VDE) for the WEST tokamak [11, 9, 10]is simulated and comparison with experimental results from shot 59695 are used to validate the simulation. Eventually the code is coupled to a magnetic feedback controller [38] enabling the simulation of a prescribed plasma scenario. The results for an X-point formation scenario in the WEST tokamak are presented as an illustration of the efficiency of the developed numerical methods.

## 2 Continuous model

The model for plasma equilibrium evolution presented here is made up of 3 main state equations: the 2D quasi-static equilibrium evolution equation for the poloidal flux $\psi$ described and used in [6, 4, 26, 5, 25, 8], the resistive current diffusion equation given as an evolution equation for the poloidal current function $f(\psi, t)$ provided in [25] and an equation for the electron energy transport inside the plasma. The knowlegde of the electron temperature $T_e$ enables us to compute consistently the pressure term needed in the equilibrium equation as well as the parallel resistivity $\eta^\parallel$ and the bootstrap current which appear in the evolution equation for f.

### 2.1 Equilibrium

Let us introduce a cylindrical $(r, \phi, z)$ coordinate system and assume axisymmetry. Maxwell's equations together with force balance in the plasma reduce to the following equation for the poloidal flux $\psi(r, z)$ in the poloidal section of the tokamak contained in the positive half plane $D = [0, \infty] \times [-\infty, \infty]$:

$$ -\Delta^* \psi = j_{\phi}, \qquad (1) $$

where $j_{\phi}$ is the toroidal component of the current density $\mathbf{j}$, and the second order elliptic differential operator $\Delta^*$ is defined by

$$ \Delta^* \psi := \nabla \cdot \left( \frac{1}{\mu(\psi)r} \nabla \psi \right), \qquad (2) $$

with $\nabla$ the 2D operator in the $(r, z)$-plane and $\mu(\psi)$ is the magnetic permeability. The magnetic field decomposes in a poloidal and toroidal component

$$ \mathbf{B} = \mathbf{B}_p + \mathbf{B}_{\phi} \text{ with } \mathbf{B}_p = \frac{1}{r} \left( -\frac{\partial \psi}{\partial z} \mathbf{e}_r + \frac{\partial \psi}{\partial r} \mathbf{e}_z \right) \text{ and } \mathbf{B}_{\phi} = \frac{f}{r} \mathbf{e}_{\phi}. $$

The geometry of the tokamak (see Fig. 2) determines various sub-domains that are used to set the expression of $j_{\phi}$ accordingly.

- $\Omega_{Fe} \subset D$ denotes those parts of $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) = \mu_{Fe}(|\nabla\psi|^2 r^{-2}) \ge \mu_0$; elsewhere, that is in $\Omega_{Fe}^c$, $\mu = \mu_0$ the permeability of vacuum;

- $\Omega_{c_i} \subset D$, $1 \le i \le N_c$, denotes the intersection with the poloidal plane of the $i$th poloidal field coil carrying currents and $N_c$ denotes the number of coils. The $i$th coil has cross section area $|\Omega_{c_i}|$ and carries a total current $I_i$;
$$ \int_{\Omega_p(\psi)} \frac{df(\psi, t)}{dt} \frac{g(\psi)}{r} dr dz + \int_{\Omega_p(\psi)} \frac{d\psi}{dt} \frac{f(\psi, t)}{r} g'(\psi) dr dz \\ + \int_{\Omega_p(\psi)} \eta^{\parallel}(rp'f(\psi, t) + \frac{f^2(\psi, t) + |\nabla\psi|^2}{\mu_0 r} f'(\psi, t) - r\mathbf{j}_{ni} \cdot \mathbf{B})g'(\psi) dr dz = 0, \quad (4) $$

for any test function $g(\psi)$ with $g(\psi_b(t)) = 0$. The value of $f$ at the plasma boundary is given by $f(\psi_b(t), t) = f_b(t) = B_0 r_0$ with $B_0$ the external toroidal magnetic field at $r_0$ center of the tokamak vacuum vessel. The scalar product term $\mathbf{j}_{ni} \cdot \mathbf{B}$ takes into account the non-inductive bootstrap current density $\mathbf{j}_{ni}$ which differs from Ohmic current and represents the source terms in the plasma current which are not created by variations of the poloidal flux. Since this current density is not well known [19] takes the approximation $\mathbf{j}_{ni} \cdot \mathbf{B} \approx j_{ni}(\psi)|\mathbf{B}|$ and uses an arbitrary analytic expression for $j_{ni}(\psi)$. The same is done for the resistivity $\eta^{\parallel}$ (the parallel component coefficient of the anisotropic resistivity tensor of the plasma). In the next paragraph we complete the model with an evolution equation for the electron temperature which enables us to approximate more precisely $\mathbf{j}_{ni} \cdot \mathbf{B}$ with its average over magnetic surfaces $\langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle$, and to compute the resistivity $\eta^{\parallel}$.

## 2.3. Electron energy transport, resistivity and bootstrap current

We use a simplified transport model including an electron energy transport equation only. Such a simplification is justified by the fact that the most important parameters for the plasma state description are the poloidal flux and the electron temperature [34, 35]. At this stage our goal is not to have a complete transport model such as those found e.g. in [6, 10, 36, 37, 38, 39] but rather to have a minimal model enabling to compute in a realistic way the resistivity $\eta^{\parallel}$ and the non-inductive current term $\langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle$.

In the plasma the conservation of energy of electrons [6, 2] reads

$$ \frac{3}{2} \dot{p}_e + \mathrm{div}(\mathbf{Q}_e + \frac{5}{2} p_e \mathbf{u}_e) = \mathbf{j} \cdot \mathbf{E} - \mathbf{u}_i \cdot \mathrm{grad} \, p_i - S_s + S_a, \quad (5) $$

where $p_e = n_e T_e^J$ is the pressure for electrons, $T_e^J$ the electron temperature in $[J]$, $n_e$ the electron density in $[m^{-3}]$, $p_i = n_i T_i^J$ the pressure for ions. The quantity $\mathbf{Q}_e$ is the heat flux of the electrons and $\mathbf{u}_e$ the electron velocity. The scalar product $\mathbf{j} \cdot \mathbf{E}$ is the Ohmic source term with $\mathbf{E}$ the electric field and $\mathbf{j}$ the current density. Finally, $S_s$ represents sink terms due to the equipartition
- $S \subset D$, denotes the part of domain $S = \cup_j \Omega_{ps_j}$, where are located the $N_{ps}$ passive structures $\Omega_{ps_j}$, characterized by an electric conductivity $\sigma_j \neq 0$, $j = 1, ..., N_{ps}$ (note that the tokamak vacuum vessel is a passive structure too);

- $\Omega_L \subset D$, denotes the domain bounded by the limiter, thus the domain accessible to 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 Grad-Shafranov equation [24, 45] is satisfied in this domain.

The equilibrium of the plasma in a tokamak has to satisfy at each instant the following non-linear initial boundary value problem: starting from an initial condition $\psi_0(r, z)$ find $\psi(r, z, t)$ such that

$$ -\Delta^*\psi = j_\phi = \begin{cases} rp'(\psi, t) + \frac{1}{\mu_0 r} ff'(\psi, t) & \text{in } \Omega_p(\psi), \\ \frac{1}{|\Omega_{c_i}|} I_i\left(\frac{d\psi}{dt}, t\right) & \text{in } \Omega_{c_i}, i = 1, \dots, N_c, \\ -\frac{\sigma_j}{r} \frac{d\psi}{dt} & \text{in } \Omega_{ps_j} \subset S, j = 1, \dots, N_{ps}, \\ 0 & \text{elsewhere,} \end{cases} \quad (3) \\ I_i\left(\frac{d\psi}{dt}, t\right) = \sum_{j=1}^{N_s} S_{ij}V_j(t) + \sum_{j=1}^{N_c} R_{ij} \frac{1}{|\Omega_{c_j}|} \int_{\Omega_{c_j}} \frac{d\psi}{dt} dr dz, \quad i = 1, \dots, N_c $$

with boundary conditions $\psi(0, z, t) = 0$ and $\lim_{||(r,z)||_2 \to +\infty} \psi(r, z, t) = 0$.

$\frac{d\psi}{dt}$ is the time derivative of $\psi$ at fixed $(r, z)$. The current $I_i$ relates to voltages $V_j(t)$ in the power supplies and to self and mutual induction via electric circuit equations (see [26] Appendix A and [5] Appendix B). $N_s$ is the number of power supplies, matrix $S$ has dimension $N_c \times N_s$ and matrix $R$ has dimension $N_c \times N_c$. The plasma domain $\Omega_p(\psi)$ is unknown and depends non-linearly on the poloidal flux $\psi$ (it is a free-boundary problem). The boundary of $\Omega_p(\psi)$ is an iso-contour of $\psi$ and is either tangent the boundary of $\Omega_L$ (limiter configuration) or contains one or more saddle points (X-points) of $\psi$ (divertor configuration). The value of $\psi$ at the magnetic axis is denoted $\psi_a$ and on the plasma boundary by $\psi_b$ (the subscript $b$ is used for values at the plasma boundary).

The two functions $p'$ and $ff'$ (where the superscript $'$ denotes here derivation with respect to $\psi$) are the usual terms in the Grad-Shafranov equation. The resistive diffusion equation presented in the next Section provides a way to compute $f$. The pressure term $p$ is computed from the electron temperature as explained further below in Section 2.3.

## 2.2 Resistive current diffusion

In the above equilibrium equation the Grad-Shafranov pressure term $p'$ and diamagnetic term $ff'$ are not known and have to be provided. The pressure term does not lead to significant numerical difficulties. It is not the case for the $ff'$ term. In the literature this term is usually computed from the 1D flux surface average $\psi$ evolution equation, classically called the resistive current diffusion equation, together with the flux surface average Grad-Shafranov equation. With this approach the poloidal flux has two representations, a 2D representation $\psi(r, z)$ from the equilibrium model and a 1D representation $\psi(\rho)$ where $\rho$ is a label for the flux surfaces. This
between electrons and ions and radiation, $S_a$ represents auxiliary sources and $\mathbf{u}_i$ is the ion velocity. Due to strong anisotropy in the plasma quantities $p_e, n_e, T_e^J$ and $p_i, n_i, T_i^J$ are assumed to be constant on the magnetic surfaces. Moreover in what follows it is assumed that electron and ion temperatures and densities are proportional i.e $n_i = c_{ie}n_e$ and $T_i = r_{ie}T_e$, and that $\mathbf{u}_e = \mathbf{u}_i$ [40].

Let us now introduce the radial coordinate $\rho = \sqrt{\frac{\phi}{\pi B_0}}$ with $\phi$ the toroidal flux. The magnetic field $B_0$ is assumed constant. In order to obtain a 1D equation for electron temperature using $\rho$ as the space variable (5) is averaged over the magnetic surfaces [2] leading to (6)

$$ \frac{3}{2} \frac{\partial}{\partial t} (n_e T_e^J V^{5/3}) + V^{2/3} \frac{\partial}{\partial \rho} (q_e + \frac{5}{2} T_e^J \Gamma_e) = V^{5/3} (\langle \mathbf{j} \cdot \mathbf{E} \rangle - \langle \mathbf{u}_i \cdot \nabla p_i \rangle - \langle \mathbf{u}_\rho \cdot \nabla p_e \rangle - \langle S_s \rangle + \langle S_a \rangle), \quad (6) $$

where partial derivative with respect to time is taken at fixed $\rho$, $\frac{\partial}{\partial t}|_\rho$. The plasma volume enclosed in a flux surface of label $\rho \in [0, \rho_b]$ is $V = V(\rho, t)$ and $V'$ denotes its derivative with respect to $\rho$. The fluxes are defined by $q_e = V' \langle \mathbf{Q}_e \cdot \nabla \rho \rangle$ and $\Gamma_e = V' \langle n_e(\mathbf{u}_e - \mathbf{u}_\rho) \cdot \nabla \rho \rangle$.

The Ohmic source term can be written from [2, eqn. (VI.56) and (VI.57)]

$$ \langle \mathbf{j} \cdot \mathbf{E} \rangle = - \frac{\langle \mathbf{E} \cdot \mathbf{B} \rangle}{\mu_0 f C_3} \frac{\partial}{\partial \rho} (C_2 \frac{\partial \psi}{\partial \rho}) + \langle \mathbf{u}_\rho \cdot \nabla p \rangle \quad (7) $$

where $C_2 = V' \langle |\nabla \rho|^2/r^2 \rangle$ and $C_3 = V' \langle 1/r^2 \rangle$. Let us also define $C_1 = V' \langle |\nabla \rho|^2 \rangle$. This expression involving a second derivative of $\psi(\rho)$ is not very convenient numerically. Contrary to [2] in which $\langle \mathbf{E} \cdot \mathbf{B} \rangle$ is also expressed using derivatives of $\psi(\rho)$ we prefer an expression involving $p', f$ and $f'$. From Ohm's law we have $\mathbf{E} \cdot \mathbf{B} = \eta^{\parallel}(\mathbf{j} - \mathbf{j}_{ni}) \cdot \mathbf{B}$ and since $\mathbf{j} = (rp' + \frac{1}{\mu_0 r} f f')\mathbf{e}_\phi + \frac{f'}{\mu_0} \mathbf{B}_p$ it results that

$$ \langle \mathbf{E} \cdot \mathbf{B} \rangle = \eta^{\parallel} (p'f + \frac{1}{\mu_0} (\langle \frac{1}{r^2} \rangle f^2 + \langle \frac{|\nabla \psi|^2}{r^2} \rangle) f' - \langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle) \quad (8) $$

Using [2, eqn. (VI.42)], and the fact that $\frac{-1}{\mu_0 V'} \frac{\partial}{\partial \rho} (C_2 \frac{\partial \psi}{\partial \rho}) = \langle \frac{j_\phi}{r} \rangle$, we have

$$ \langle \mathbf{j} \cdot \mathbf{E} \rangle = Q_{ohm} + \langle \mathbf{u}_\rho \cdot \nabla p \rangle \quad (9) $$
resistive current diffusion equation is usually part of a transport model which is evolved in parallel to the equilibrium model. With such a procedure it is very difficult to ensure numerically the consistency of these 2 representations of the poloidal flux and to provide a robust simulator. This is the main reason why the idea of using directly an evolution equation for $f$ is attractive. A flux surface average evolution equation for $f$ can be found for example in [47, eqn.(20)] [12, eqn.(24)] but because of the averaging process and of the so-called geometric coefficients which appear it is difficult to solve the equilibrium equation together with the $f$ equation as a single system of equations using Newton's method to resolve the non-linearity. The formulation proposed in [25] solves this issue. From Ohm's law and Faraday's law [25] derives an evolution equation for $f$ given in weak form using 1D test function $g(\psi)$ but integrating over the 2D plasma domain. With this approach no geometric coefficients are needed and the system of equations can be solved in $(\psi, f)$ directly. The evolution problem for $f$ reads as follows: starting from an initial condition $f_0(\psi_0)$ find $f(\psi, t)$ such that

$$ \int_{\Omega_p(\psi)} \frac{df(\psi, t)}{dt} \frac{g(\psi)}{r} drdz + \int_{\Omega_p(\psi)} \frac{d\psi}{dt} \frac{f(\psi, t)}{r} g'(\psi) drdz \\ + \int_{\Omega_p(\psi)} \eta^{\parallel}(rp'f(\psi, t) + \frac{f^2(\psi, t) + |\nabla\psi|^2}{\mu_0 r}) f'(\psi, t) - r\mathbf{j}_{ni} \cdot \mathbf{B})g'(\psi) drdz = 0, \quad (4) $$

for any test function $g(\psi)$ with $g(\psi_b(t)) = 0$. The value of $f$ at the plasma boundary is given by $f(\psi_b(t), t) = f_b(t) = B_0 r_0$ with $B_0$ the external toroidal magnetic field at $r_0$ center of the tokamak vacuum vessel. The scalar product term $\mathbf{j}_{ni} \cdot \mathbf{B}$ takes into account the non-inductive bootstrap current density $\mathbf{j}_{ni}$ which differs from ohmic current and represents the source terms in the plasma current which are not created by variations of the poloidal flux. Since this current density is not well known [25] takes the approximation $\mathbf{j}_{ni} \cdot \mathbf{B} \approx j_{ni}(\psi)|\mathbf{B}|$ and uses an arbitrary analytic expression for $j_{ni}(\psi)$. The same is done for the resistivity $\eta^{\parallel}$ (the parallel component coefficient of the anisotropic resistivity tensor of the plasma). In the next paragraph we complete the model with an evolution equation for the electron temperature which enables us to approximate more precisely $\mathbf{j}_{ni} \cdot \mathbf{B}$ with its average over magnetic surfaces $\langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle$, and to compute the resistivity $\eta^{\parallel}$.

## 2.3 Electron energy transport, resistivity and bootstrap current

We use a simplified transport model including an electron energy transport equation only. Such a simplification is justified by the fact that the most important parameters for the plasma state description are the poloidal flux and the electron temperature [50, 21]. At this stage our goal is not to have a complete transport model such as those found e.g. in [28, 6, 41, 2, 14, 32] but rather to have a minimal model enabling to compute in a realistic way the resistivity $\eta^{\parallel}$ and the non-inductive current term $\langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle$.

In the plasma the conservation of energy of electrons [28, 4] reads

$$ \frac{3}{2} \dot{p}_e + \text{div}(\mathbf{Q}_e + \frac{5}{2} p_e \mathbf{u}_e) = \mathbf{j} \cdot \mathbf{E} - \mathbf{u}_i \cdot \text{grad} p_i - S_s + S_a, \quad (5) $$

where $p_e = n_e T_e^J$ is the pressure for electrons, $T_e^J$ the electron temperature in $[J]$, $n_e$ the electron density in $[m^{-3}]$, $p_i = n_i T_i^J$ the pressure for ions. $\mathbf{Q}_e$ is the heat flux of the electrons and $\mathbf{u}_e$ the electron velocity. $\mathbf{j} \cdot \mathbf{E}$ is the Ohmic source term with $\mathbf{E}$ the electric field and $\mathbf{j}$ the current density. $S_s$ represents sink terms due to the equipartition between electrons and ions and radiation, $S_a$ represents auxiliary sources and $\mathbf{u}_i$ is the ion velocity. Due to strong anisotropy in the plasma quantities $p_e, n_e, T_e^J$ and $p_i, n_i, T_i^J$ are assumed to be constant on the
with

$$Q_{ohm} = \frac{\eta^{\parallel} \langle j_{\phi}/r \rangle}{f \langle 1/r^2 \rangle} (p'f + \frac{1}{\mu_0}(\langle 1/r^2 \rangle f^2 + \langle |\nabla\psi|^2/r^2 \rangle) f' - \langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle) \quad (10)$$

216 Injecting (9) into (6) and noting $Q_s = \langle S_s \rangle$, $Q_a = \langle S_a \rangle$ leads to

$$\begin{aligned} & \frac{3}{2} \frac{\partial}{\partial t} (n_e T_e^J V'^{5/3}) + V'^{2/3} \frac{\partial}{\partial \rho} (q_e + \frac{5}{2} T_e^J \Gamma_e) \\ & = V'^{5/3} (Q_{ohm} - \langle \mathbf{u}_i - \mathbf{u}_\rho \rangle \cdot \nabla p_i) - Q_s + Q_a. \end{aligned} \quad (11)$$

The assumptions concerning the fluxes [40] are $n_e(\mathbf{u}_e - \mathbf{u}_\rho) = -\chi_{ne}\nabla n_e$ and $Q_e = -\chi_{T_e}n_e\nabla T_e$ with $\chi_{T_e}$ and $\chi_{n_e}$ diffusion coefficients profiles. This leads to

$$\frac{3}{2} \frac{\partial}{\partial t} (n_e T_e^J V'^{5/3}) - V'^{2/3} \frac{\partial}{\partial \rho} (\chi_{T_e} C_1 n_e \frac{\partial T_e^J}{\partial \rho} + \frac{5}{2} \chi_{n_e} C_1 \frac{\partial n_e}{\partial \rho} T_e^J) = V'^{5/3} Q \quad (12)$$

220 with $Q = (Q_{ohm} + Q_i - Q_s + Q_a)$, $Q_i = \frac{\chi_{n_e} C_1}{n_e} \frac{\partial n_e}{\partial \rho} \frac{\partial p_i}{\partial \rho}$ and $p_i = \frac{1}{1+\frac{1}{r_{ie}c_{ie}}p}$.

In order to work on the fixed interval $[0, 1]$ we rewrite (12) in the normalized variable $x = \frac{\rho}{\rho_b} = \rho_N$ using $\frac{\partial}{\partial t}|_{\rho} = \frac{\partial}{\partial t}|_x - \frac{x}{\rho_b} \frac{d\rho_b}{dt} \frac{\partial}{\partial x}$. We use the same symbol for a quantity seen as function of $\rho$ or $x$. All left hand side spatial derivative terms in (12) can be gathered using the identity [41, eqn. (0.6)]:

$$\begin{aligned} & \frac{3}{2} V'^{-5/3} \frac{\partial (n_e T_e^J V'^{5/3})}{\partial t} \Big|_{\rho} = \frac{3}{2} V'^{-5/3} \frac{\partial (n_e T_e^J V'^{5/3})}{\partial t} \Big|_x \\ & - V'^{-1} \frac{\partial (V'^{\frac{3}{2}} k_{\phi} \rho n_e T_e^J)}{\partial \rho} + \left( \frac{3}{2} k_{\rho} - k_{\phi} \rho \frac{\partial (\log V')}{\partial \rho} \right) (n_e T_e^J), \end{aligned} \quad (13)$$

with $k_\phi = \frac{1}{2\phi_b} \frac{d\phi_b}{dt}$, $k_\rho = \frac{1}{\rho_b} \frac{d\rho_b}{dt}$ and $k_\phi = 0$ since here $\frac{dB_0}{dt} = 0$. We also use the identity [41, eqn. (0.7)]: $\rho \frac{\partial(\log V')}{\partial \rho} = 1 - \rho \frac{\partial(\log(f < 1/r^2 >))}{\partial \rho}$. Finally we define $T_e^J = eT_e$ where $e$ is the electron charge and $T_e$ in $[eV]$, and $n_e = n_0 N_e$ with $n_0 = 1.0 \times 10^{19} [m^{-3}]$ to obtain

$$\begin{aligned} & \frac{3}{2} \frac{\partial (N_e T_e V'^{5/3})}{\partial t} - V'^{2/3} \frac{1}{\rho_b} \frac{\partial}{\partial x} \left( \left[ \frac{1}{\rho_b} C_1 N_e \chi_{T_e} \right] \frac{\partial T_e}{\partial x} \right) \\ & - V'^{2/3} \frac{1}{\rho_b} \frac{\partial}{\partial x} \left( \left[ C_1 \frac{5}{2} \chi_{n_e} \frac{\partial N_e}{\partial \rho} + \frac{3}{2} V' k_\rho \rho N_e \right] T_e \right) \\ & + V'^{5/3} \left[ k_\rho \rho \left( \frac{\partial(\log(f < 1/r^2 >))}{\partial \rho} + \frac{1}{2} \right) \right] N_e T_e = V'^{5/3} \frac{1}{n_0 e} Q, \end{aligned} \quad (14)$$
magnetic surfaces. Moreover in what follows it is assumed that electron and ion temperatures and densities are proportional i.e $n_i = c_{ie}n_e$ and $T_i = r_{ie}T_e$, and that $\mathbf{u}_e = \mathbf{u}_i$ [35].

Let us now introduce the radial coordinate $\rho = \sqrt{\frac{\phi}{\pi B_0}}$ with $\phi$ the toroidal flux. The magnetic field $B_0$ is assumed constant. In order to obtain a 1D equation for electron temperature using $\rho$ as the space variable (5) is averaged over the magnetic surfaces [4] leading to (6)

$$ \frac{3}{2} \frac{\partial}{\partial t} (n_e T_e^J V'^{5/3}) + V'^{2/3} \frac{\partial}{\partial \rho} (q_e + \frac{5}{2} T_e^J \Gamma_e) = V'^{5/3} (\langle \mathbf{j} \cdot \mathbf{E} \rangle - \langle \mathbf{u}_i \cdot \nabla p_i \rangle - \langle \mathbf{u}_\rho \cdot \nabla p_e \rangle - \langle S_s \rangle + \langle S_a \rangle), \quad (6) $$

where partial derivative with respect to time is taken at fixed $\rho$, $\frac{\partial}{\partial t}|_\rho$. The plasma volume enclosed in a flux surface of label $\rho \in [0, \rho_b]$ is $V = V(\rho, t)$ and $V'$ denotes its derivative with respect to $\rho$. The fluxes are defined by $q_e = V' \langle \mathbf{Q}_e \cdot \nabla \rho \rangle$ and $\Gamma_e = V' \langle n_e(\mathbf{u}_e - \mathbf{u}_\rho) \cdot \nabla \rho \rangle$. The Ohmic source term can be written from [4, eqn. (VI.56) and (VI.57)]

$$ \langle \mathbf{j} \cdot \mathbf{E} \rangle = - \frac{\langle \mathbf{E} \cdot \mathbf{B} \rangle}{\mu_0 f C_3} \frac{\partial}{\partial \rho} (C_2 \frac{\partial \psi}{\partial \rho}) + \langle \mathbf{u}_\rho \cdot \nabla p \rangle \quad (7) $$

where $C_2 = V' \langle |\nabla \rho|^2/r^2 \rangle$ and $C_3 = V' \langle 1/r^2 \rangle$. Let us also define $C_1 = V' \langle |\nabla \rho|^2 \rangle$. This expression involving a second derivative of $\psi(\rho)$ is not very convenient numerically. Contrary to [4] in which $\langle \mathbf{E} \cdot \mathbf{B} \rangle$ is also expressed using derivatives of $\psi(\rho)$ we prefer an expression involving $p'$, $f$ and $f'$. From Ohm's law we have $\mathbf{E} \cdot \mathbf{B} = \eta^{\parallel}(\mathbf{j} - \mathbf{j}_{ni}) \cdot \mathbf{B}$ and since $\mathbf{j} = (rp' + \frac{1}{\mu_0 r} f f')\mathbf{e}_\phi + \frac{f'}{\mu_0} \mathbf{B}_\rho$ it results that

$$ \langle \mathbf{E} \cdot \mathbf{B} \rangle = \eta^{\parallel} (p'f + \frac{1}{\mu_0} (\langle \frac{1}{r^2} \rangle f^2 + \langle \frac{|\nabla \psi|^2}{r^2} \rangle) f' - \langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle) \quad (8) $$

Using [4, eqn. (VI.42)], $\frac{-1}{\mu_0 V'} \frac{\partial}{\partial \rho} (C_2 \frac{\partial \psi}{\partial \rho}) = \langle \frac{j_\phi}{r} \rangle$, we have

$$ \langle \mathbf{j} \cdot \mathbf{E} \rangle = Q_{ohm} + \langle \mathbf{u}_\rho \cdot \nabla p \rangle \quad (9) $$

with

$$ Q_{ohm} = \frac{\eta^{\parallel} \langle j_{\phi}/r \rangle}{f \langle 1/r^2 \rangle} (p'f + \frac{1}{\mu_0} (\langle 1/r^2 \rangle f^2 + \langle |\nabla \psi|^2/r^2 \rangle) f' - \langle j_{ni} \cdot \mathbf{B} \rangle) \quad (10) $$

Injecting (9) into (6) and noting $Q_s = \langle S_s \rangle$, $Q_a = \langle S_a \rangle$ leads to

$$ \frac{3}{2} \frac{\partial}{\partial t} (n_e T_e^J V'^{5/3}) + V'^{2/3} \frac{\partial}{\partial \rho} (q_e + \frac{5}{2} T_e^J \Gamma_e) = V'^{5/3} (Q_{ohm} - \langle \mathbf{u}_i - \mathbf{u}_\rho \rangle \cdot \nabla p_i - Q_s + Q_a) \quad (11) $$

The assumptions concerning the fluxes [35] are $n_e(\mathbf{u}_e - \mathbf{u}_\rho) = -\chi_{ne} \nabla n_e$ and $\mathbf{Q}_e = -\chi_{T_e} n_e \nabla T_e$ with $\chi_{T_e}$ and $\chi_{n_e}$ diffusion coefficients profiles. This leads to

$$ \frac{3}{2} \frac{\partial}{\partial t} (n_e T_e^J V'^{5/3}) - V'^{2/3} \frac{\partial}{\partial \rho} (\chi_{T_e} C_1 n_e \frac{\partial T_e^J}{\partial \rho} + \frac{5}{2} \chi_{n_e} C_1 \frac{\partial n_e}{\partial \rho} T_e^J) = V'^{5/3} Q \quad (12) $$

with $Q = (Q_{ohm} + Q_i - Q_s + Q_a)$, $Q_i = \frac{\chi_{n_e} C_1}{n_e} \frac{\partial n_e}{\partial \rho} \frac{\partial p_i}{\partial \rho}$ and $p_i = \frac{1}{1 + \frac{1}{r_{ie} c_{ie}} p}$.
where partial derivative with respect to time is taken at fixed $x, \frac{\partial}{\partial t}|_x$

The sink term is $Q_s = Q_{ei} + Q_{rad}$. The expressions used for the electron-ion energy exchange $Q_{ei}$ and for the radiation loss term $Q_{rad}$ are given in Appendix B. Other external auxiliary sources $Q_a$ are here modelized with a Gaussian

$$Q_a = \frac{P_a}{\rho_b \int_0^1 \exp(-(x-x_0)^2/\sigma^2)V'dx} \exp(-(x-x_0)^2/\sigma^2) \quad (15)$$

where $P_a$ is the total auxiliary power delivered in the plasma.

In order to lighten notations we write (14) as:

$$\frac{\partial(BT_e)}{\partial t} - C \frac{\partial}{\partial x} (D \frac{\partial T_e}{\partial x} - ET_e) + FT_e = G. \quad (16)$$

Coefficients $B, C, D, E, F, G$ depend on the plasma equilibrium configuration through several quantities such as its volume and its variations, and average geometric coefficients. They also depend on the electron density profile $n_e$ [$m^{-3}$] which is assumed to be given in this study. The source term $G$ also depends on $\eta^{\parallel}$ and $<\mathbf{j}_{ni} \cdot \mathbf{B}>$.

Equation (16) degenerates at the center of the plasma $x = 0$ where $V'(0) = 0$ and no boundary condition is required. On the plasma boundary $x = 1$ the boundary value $T_{e,b}$ is prescribed.

From the knowledge of the electron temperature profile $T_e$ we compute on the one hand the pressure profile needed in the equilibrium equation, $p = n_e T_e (1 + r_{ie} c_{ie}) e$, and on the other hand the resistivity profile $\eta^{\parallel}$ and the bootstrap current profile $<\mathbf{j}_{ni} \cdot \mathbf{B}>$ with the expressions provided in [42, 43]. The resistivity profile is given by $\eta^{\parallel} = c_{\eta} T_e^{-3/2}$ and the bootstrap current profile by $<\mathbf{j}_{ni} \cdot \mathbf{B}> = c_{bs} f p$ where coefficients $c_{\eta}$ and $c_{bs}$ depend on $T_e, n_e$ and on several other equilibrium profiles: pressure, trapped particle fraction $f_t$ profile computed from the equilibrium using Equations (33-35) in [44] and safety factor $q$.

### 3. Discrete problem

We have to solve equation (3) for $\psi(r, z, t)$, equation (4) for $f(\psi, t)$ and equation (16) for $T_e(x, t)$. The equilibrium system consists of the two first ones and the transport system in our model is reduced to the third equation. Note that this differs from what is found in the literature where the transport
In order to work on the fixed interval [0, 1] we rewrite (12) in the normalized variable $x = \frac{\rho}{\rho_b} = \rho_N$ using $\frac{\partial}{\partial t}|_{\rho} = \frac{\partial}{\partial t}|_x - \frac{x}{\rho_b} \frac{d\rho_b}{dt} \frac{\partial}{\partial x}$. We use the same symbol for a quantity seen as function of $\rho$ or $x$. All left hand side spatial derivative terms in (12) can be gathered using the identity [33, eqn. (0.6)]:

$$ \frac{3}{2} V'^{-5/3} \frac{\partial(n_e T_e^J V'^{5/3})}{\partial t} |_{\rho} = \frac{3}{2} V'^{-5/3} \frac{\partial(n_e T_e^J V'^{5/3})}{\partial t} |_x - V'^{-1} \frac{\partial(V'^{3/2} k_{\phi} \rho n_e T_e^J)}{\partial \rho} + \left( \frac{3}{2} k_{\rho} - k_{\phi} \rho \frac{\partial(\log V')}{\partial \rho} \right) (n_e T_e^J), \quad (13) $$

with $k_{\phi} = \frac{1}{2\phi_b} \frac{d\phi_b}{dt}$, $k_{\rho} = \frac{1}{\rho_b} \frac{d\rho_b}{dt}$ and $k_{\phi} = 0$ since here $\frac{dB_0}{dt} = 0$. We also use the identity [33, eqn. (0.7)]: $\rho \frac{\partial(\log V')}{\partial \rho} = 1 - \rho \frac{\partial(\log(f < 1/r^2 >))}{\partial \rho}$. Finally we define $T_e^J = eT_e$ where $e$ is the electron charge and $T_e$ in [eV], and $n_e = n_0 N_e$ with $n_0 = 1.0e19 [m^{-3}]$ to obtain

$$ \begin{aligned} & \frac{3}{2} \frac{\partial(N_e T_e V'^{5/3})}{\partial t} - V'^{2/3} \frac{1}{\rho_b} \frac{\partial}{\partial x} \left( \left[ \frac{1}{\rho_b} C_1 N_e \chi_{T_e} \right] \frac{\partial T_e}{\partial x} \right) \\ & - V'^{2/3} \frac{1}{\rho_b} \frac{\partial}{\partial x} \left( \left[ C_1 \frac{5}{2} \chi_{n_e} \frac{\partial N_e}{\partial \rho} + \frac{3}{2} V' k_{\rho} \rho N_e \right] T_e \right) \\ & + V'^{5/3} \left[ k_{\rho} \rho \left( \frac{\partial(\log(f < 1/r^2>))}{\partial \rho} + \frac{1}{2} \right) \right] N_e T_e = V'^{5/3} \frac{1}{n_0 e} Q, \end{aligned} \quad (14) $$

where partial derivative with respect to time is taken at fixed $x$, $\frac{\partial}{\partial t}|_x$

The sink term is $Q_s = Q_{ei} + Q_{rad}$. The expressions used for the electron-ion energy exchange $Q_{ei}$ and for the radiation loss term $Q_{rad}$ are given in B. Other external auxiliary sources $Q_a$ are here modelized with a Gaussian

$$ Q_a = \frac{P_a}{\rho_b \int_0^1 \exp(-(x-x_0)^2/\sigma^2) V' dx} \exp(-(x-x_0)^2/\sigma^2) \quad (15) $$

where $P_a$ is the total auxiliary power delivered in the plasma.

In order to lighten notations we write (14) as:

$$ \frac{\partial(BT_e)}{\partial t} - C \frac{\partial}{\partial x} (D \frac{\partial T_e}{\partial x} - ET_e) + FT_e = G, \quad (16) $$

Coefficients $B, C, D, E, F, G$ depend on the plasma equilibrium configuration through several quantities such as its volume and its variations, and average geometric coefficients. They also depend on the electron density profile $n_e$ [$m^{-3}$] which is assumed to be given in this study. The source term $G$ also depends on $\eta^{\parallel}$ and $<\mathbf{j}_{ni} \cdot \mathbf{B}>$.

Equation (16) degenerates at the center of the plasma $x = 0$ where $V'(0) = 0$ and no boundary condition is required. On the plasma boundary $x = 1$ the boundary value $T_{e,b}$ is prescribed.

From the knowledge of the electron temperature profile $T_e$ we compute on one hand the pressure profile needed in the equilibrium equation, $p = n_e T_e (1 + r_{ie} c_{ie}) e$, and on the other hand the resistivity profile $\eta^{\parallel}$ and the bootstrap current profile $<\mathbf{j}_{ni} \cdot \mathbf{B}>$ with the expressions provided in [43, 44]. The resistivity profile is given by $\eta^{\parallel} = c_{\eta} T_e^{-3/2}$ and the bootstrap current profile by $<\mathbf{j}_{ni} \cdot \mathbf{B}> = c_{bs} f p$ where coefficients $c_{\eta}$ and $c_{bs}$ depend on $T_e, n_e$ and on several other equilibrium profiles: pressure, trapped particle fraction $f_t$ profile computed from the equilibrium using Equations (33-35) from [42] and safety factor $q$.
system usually includes the 1D resistive equation for $\psi(\rho, t)$ from which the $ff'$ term needed in (3) is computed. The equations for $\psi$ and $f$ are solved together as a non-linear system using Newton's method at each time step as explained below. This equilibrium system requires the knowledge of pressure $p$, resistivity $\eta^{\|}$ and bootstrap current $\langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle$. These quantities are deduced from the $T_e$ equation. On the other hand to solve the $T_e$ equation it is necessary to know a number of geometric coefficients deduced from the equilibrium. Solving the three equations simultaneously is out of reach of this work because Newton's method would require to compute the derivatives of these geometric coefficients with respect to all three variables which seems hardly possible. Therefore we iterate between the equilibrium system and the $T_e$ equation if necessary.

### 3.1. The equilibrium ($\psi, f$) system

The equilibrium ($\psi, f$) system is discretized in space using FEM for $\psi$ and the Galerkin formulation in which consists the equation for $f$. In [19] a classical $C^0$ piecewise linear Lagrange P1 FEM method is proposed for the approximation of $\psi$ and $f$ is decomposed in the space of polynomials of degree $k$ using Legendre polynomials as a basis. We have implemented this discretization at first but numerical issues (see Appendix A) have lead us to propose a new discretization method.

For the approximation of $\psi$ we adopt a non-overlapping mortar element approach introduced in our previous work [21, 22]. It couples $C^0$ P1 FEM in an annular domain that does not contain the plasma and $C^1$ piecewise cubic reduced Hsieh-Clough-Tocher (rHCT) FEM in an internal domain containing the plasma. This enables us to achieve high regularity for the representation of $\psi$ in the plasma without having to pay the numerical cost of using rHCT FEM in the whole domain and in particular outside the plasma domain where it is not needed.

We work in a domain $\Omega \subset D$, known as the ABB domain (Albanese-Blum-Barbieri [32]) delimited by a semi-circle $\gamma$ of radius $\rho_{\gamma} > 0$ including $\Omega_L \cup \Omega_{\text{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$.

In order to formulate model (3) in weak form we introduce a non-overlapping domain decomposition framework by setting $\Omega = \Omega^{\text{in}} \cup \Omega^{\text{ex}}$ where $\Omega^{\text{in}}$ is a bounded domain containing $\Omega_L$ 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}}$ and $\Omega^{\text{ex}}$, on which we impose the continuity of $\psi$, in a weak sense, through a mortar-like $L^2$ projection [45]. Moreover we consider that all coils $\Omega_{c_i}$ are
# 3 Discrete problem

The equations to be solved are equation (3) for $\psi(r, z, t)$, equation (4) for $f(\psi, t)$ and equation (16) for $T_e(x, t)$. The equilibrium system consists of the two first ones and the transport system in our model is reduced to the third equation. Note that this differs from what is found in the literature where the transport system usually includes the 1D resistive equation for $\psi(\rho, t)$ from which the $ff'$ term needed in (3) is computed. The equations for $\psi$ and $f$ are solved together as a non-linear system thanks to Newton's method at each time step as explained below. This equilibrium system requires the knowledge of pressure $p$, resistivity $\eta^{\parallel}$ and bootstrap current $<\mathbf{j}_{ni} \cdot \mathbf{B}>$. These quantities are deduced from the $T_e$ equation. On the other hand to solve the $T_e$ equation it is necessary to know a number of geometric coefficients deduced from the equilibrium. Solving the 3 equations simultaneously is out of reach of this work because Newton's method would require to compute the derivatives of these geometric coefficients with respect to all three variables which seems hardly possible. Therefore we iterate between the equilibrium system and the $T_e$ equation if necessary.

## 3.1 The equilibrium $(\psi, f)$ system

The equilibrium $(\psi, f)$ system is discretized in space using FEM for $\psi$ and the Galerkin formulation in which consists the equation for $f$. In [25] a classical $C^0$ piecewise linear Lagrange P1 FEM method is proposed for the approximation of $\psi$ and $f$ is decomposed in the space of polynomials of degree $k$ using Legendre polynomials as a basis. We have implemented this discretization at first but numerical issues (see A) have lead us to propose a new discretization method.

For the approximation of $\psi$ we adopt a non-overlapping mortar element approach introduced in our previous work [17, 8]. It couples $C^0$ P1 FEM in an annular domain that does not contain the plasma and $C^1$ piecewise cubic reduced Hsieh-Clough-Tocher (rHCT) FEM in an internal domain containing the plasma.

We work in a domain $\Omega \subset \mathcal{D}$, known as the ABB domain (Albanese-Blum-Barbieri [1]) delimited by a semi-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$.

In order to formulate model (3) in weak form we introduce a non-overlapping domain decomposition framework by setting $\Omega = \Omega^{\text{in}} \cup \Omega^{\text{ex}}$ where $\Omega^{\text{in}}$ is a bounded domain containing $\Omega_L$ 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}}$ and $\Omega^{\text{ex}}$, on which we impose the continuity of $\psi$, in a weak sense, through a mortar-like $L^2$ projection [3]. Moreover we consider that all coils $\Omega_{c_i}$ are in $\Omega^{\text{ex}}$. We also have $\Omega_{Fe} \subset \Omega^{\text{ex}}$. The weak formulation of (3) reads: given $T \ge 0$, find the function $\psi: t \in [0, T] \mapsto \psi(t) = (\psi_{\text{ex}}(t), \psi_{\text{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, s) - J_p(\psi, s; f) - j_{ps}(\frac{d\psi}{dt}, s) - j_c(\frac{d\psi}{dt}, s) - \ell_c(V(t), s) = 0, \quad (17)$$

Inria
in $\Omega^{ex}$. We also have $\Omega_{Fe} \subset \Omega^{ex}$. The weak formulation of (3) 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

296

$$a(\psi, s) - J_p(\psi, s; f) - j_{ps}(\frac{d\psi}{dt}, s) - j_c(\frac{d\psi}{dt}, s) - \ell_c(V(t), s) = 0, \quad (17)$$

where

298

$$ 
\begin{align*}
a(\psi, s) &:= a^{\text{ex}}(\psi_{\text{ex}}, v) + a^{\text{in}}(\psi_{\text{in}}, w), \\
a^{\text{ex}}(\psi, v) &:= \int_{\Omega_{Fe}} \frac{1}{\mu(\psi)r} \nabla \psi \cdot \nabla v \, dr \, dz + \int_{\Omega^{\text{ex}} \setminus \Omega_{Fe}} \frac{1}{\mu_0 r} \nabla \psi \cdot \nabla v \, dr \, dz + c(\psi, v), \\
a^{\text{in}}(\psi, w) &:= \int_{\Omega^{\text{in}}} \frac{1}{\mu_0 r} \nabla \psi \cdot \nabla w \, dr \, dz, \\
J_p(\psi, s; f) &:= \int_{\Omega_p(\psi)} \left( r p'(\psi, t) + \frac{1}{\mu_0 r} f f'(\psi, t) \right) w \, dr \, dz, \\
j_{ps}(\frac{d\psi}{dt}, s) &:= - \sum_{j=1}^{N_{ps}} \int_{\Omega_{ps_j}} \frac{\sigma_j}{r} \frac{d\psi}{dt} (\chi_{\Omega^{\text{ex}}} v + \chi_{\Omega^{\text{in}}} w) \, dr \, dz, \\
j_c(\frac{d\psi}{dt}, s) &:= \sum_{i=1}^{N_c} \sum_{j=1}^{N_s} \frac{1}{|\Omega_{c_i}||\Omega_{c_j}|} \mathbf{R}_{ij} \int_{\Omega_{c_j}} \frac{d\psi}{dt} dr \, dz \int_{\Omega_{c_i}} \chi_{\Omega^{\text{ex}}} v \, dr \, dz, \\
\ell_c(V(t), s) &:= \sum_{i=1}^{N_c} \sum_{j=1}^{N_s} \frac{1}{|\Omega_{c_i}|} \mathbf{S}_{ij} V_j(t) \int_{\Omega_{c_i}} \chi_{\Omega^{\text{ex}}} v \, dr \, dz, \tag{18}
\end{align*}
$$

and $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 [28, 33]. The functional spaces are $\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}}\}$ and $\mathcal{V}_{0,\mathcal{I}} = \{(v, w) \in \mathcal{V}, v|_{\mathcal{I}} = w|_{\mathcal{I}} = 0\}$. We have made explicit the dependence of $J_p(\psi, s; f)$ in $f$ because (17) is to be solved together with (4). This formulation of the equilibrium problem in a domain decomposition framework is suitable to use different discretization methods in $\Omega^{\text{ex}}$ and $\Omega^{\text{in}}$ as in [21, 22].

300

302

304

A mortar finite element approach is applied to (17) to obtain the semi-discrete problem. Let $\tau^{\text{ex}}$ (resp. $\tau^{\text{in}}$) be a mesh of triangles that covers $\Omega^{\text{ex}}$ (resp. $\Omega^{\text{in}}$). We assume that $\mathcal{I}$ is a polygon with nodes and edges in $\tau^{\text{ex}}$. We wish to use in $\Omega^{\text{in}} \supset \Omega_L \supset \Omega_p$, 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 piecewise cubic rHCT finite element space, say $\mathcal{V}^{\text{in}}$, on $\tau^{\text{in}}$ (see [46]). This regularity is not necessary in $\Omega^{\text{ex}}$ therefore we couple rHCT finite elements in $\Omega^{\text{in}}$ with continuous piecewise 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}}$ it 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}}\}$.

306

308

310

312

314

316
where

$$
\begin{align*}
a(\psi, s) &:= a^{\text{ex}}(\psi_{\text{ex}}, v) + a^{\text{in}}(\psi_{\text{in}}, w), \\
a^{\text{ex}}(\psi, v) &:= \int_{\Omega_{Fe}} \frac{1}{\mu(\psi)} \nabla \psi \cdot \nabla v \, drdz + \int_{\Omega^{\text{ex}} \setminus \Omega_{Fe}} \frac{1}{\mu_0 r} \nabla \psi \cdot \nabla v \, drdz + c(\psi, v), \\
a^{\text{in}}(\psi, w) &:= \int_{\Omega^{\text{in}}} \frac{1}{\mu_0 r} \nabla \psi \cdot \nabla w \, drdz, \\
J_p(\psi, s; f) &:= \int_{\Omega_p(\psi)} \left( rp'(\psi, t) + \frac{1}{\mu_0 r} ff'(\psi, t) \right) w \, drdz, \tag{18} \\
j_{ps}(\frac{d\psi}{dt}, s) &:= - \sum_{j=1}^{N_{ps}} \int_{\Omega_{ps,j}} \frac{\sigma_j}{r} \frac{d\psi}{dt} (\chi_{\Omega^{\text{ex}}} v + \chi_{\Omega^{\text{in}}} w) \, drdz, \\
j_c(\frac{d\psi}{dt}, s) &:= \sum_{i=1}^{N_c} \sum_{j=1}^{N_s} \frac{1}{|\Omega_{c_i}||\Omega_{c_j}|} \mathbf{R}_{ij} \int_{\Omega_{c_j}} \frac{d\psi}{dt} \, drdz \int_{\Omega_{c_i}} \chi_{\Omega^{\text{ex}}} v \, drdz, \\
\ell_c(V(t), s) &:= \sum_{i=1}^{N_c} \sum_{j=1}^{N_s} \frac{1}{|\Omega_{c_i}|} \mathbf{S}_{ij} V_j(t) \int_{\Omega_{c_i}} \chi_{\Omega^{\text{ex}}} v \, drdz,
\end{align*}
$$

and $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 [26, 27]. The functional spaces are $V = \{(v, w) \in H^1(\Omega^{\text{ex}}) \times H^1(\Omega^{\text{in}}), v|_{\Gamma_0} = 0, v|_{\mathcal{I}} = w|_{\mathcal{I}}\}$ and $V_{0,\mathcal{I}} = \{(v, w) \in V, v|_{\mathcal{I}} = w|_{\mathcal{I}} = 0\}$. We have made explicit the dependence of $J_p(\psi, s; f)$ in $f$ because (17) is to be solved together with (4). This formulation of the equilibrium problem in a domain decomposition framework is suitable to use different discretization methods in $\Omega^{\text{ex}}$ and $\Omega^{\text{in}}$ as in [17, 8].

A mortar finite element approach is applied to (17) to obtain the semi-discrete problem. Let $\tau^{\text{ex}}$ (resp. $\tau^{\text{in}}$) be a mesh of triangles that covers $\Omega^{\text{ex}}$ (resp. $\Omega^{\text{in}}$). We assume that $\mathcal{I}$ is a polygon with nodes and edges in $\tau^{\text{ex}}$. We wish to use in $\Omega^{\text{in}} \supset \Omega_L \supset \Omega_p$, 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 piecewise cubic rHCT finite element space, say $V^{\text{in}}$, on $\tau^{\text{in}}$ (see [13]). This regularity is not necessary in $\Omega^{\text{ex}}$ therefore we couple rHCT finite elements in $\Omega^{\text{in}}$ with continuous piecewise linear finite elements, say $V^{\text{ex}}$, on $\tau^{\text{ex}}$. The finite element space over the mesh $\tau^{\text{ex}}$ is $V^{\text{ex}} = \{v \in C^0(\Omega^{\text{ex}}), v|_{\Gamma_0} = 0, v|_{\mathcal{I}} \in \mathbb{P}_1(\mathcal{I}), \forall \mathcal{I} \in \tau^{\text{ex}}\}$, whereas over $\tau^{\text{in}}$ it is $V^{\text{in}} = \{w \in C^1(\Omega^{\text{in}}), w|_{\mathcal{I}} \in P_{\text{loc}}(\mathcal{I}), \forall \mathcal{I} \in \tau^{\text{in}}\}$. The space $P_{\text{loc}}(\mathcal{I})$ reads

$$
P_{\text{loc}}(T) = \{w \in 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 at $V_m, V_\ell$ with $m, \ell \in \{1, 2, 3\} \setminus \{i\}$ and at the barycenter $G$, for each $i = 1, 2, 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$.

Let us also write $V^{\text{ex}} = V_{\circ}^{\text{ex}} \oplus \mathcal{E}V_{\circ}^{\text{ex}}$ and $V^{\text{in}} = V_{\circ}^{\text{in}} \oplus \mathcal{E}V_{\circ}^{\text{in}}$, where, for example, $V_{\circ}^{\text{ex}}$ (resp. $V_{\circ}^{\text{ex}}$) is the subspace of $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 $V_{\circ}^{\text{ex}}$ and $V_{\circ}^{\text{in}}$ have vanishing Dirichlet trace on $\mathcal{I}$. Let us define

$$
V_h = \{(u_h^{\text{in}}, u_h^{\text{ex}}) \in V^{\text{in}} \times V^{\text{ex}}, u_h^{\text{ex}}|_{\Gamma_0} = 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 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 semi-discrete evolution problem is obtained by applying the mortar finite element approach to (17). The semi-discrete problem 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 V_h$ such that $\psi_h(0) = \psi_{0,h}$ and, $\forall s_h = (v_h, w_h) \in V_{\circ}^{\text{ex}} \times V_{\circ}^{\text{in}}$, it holds

$$
a(\psi_h, s_h) - J_p(\psi_h, s_h; f(\psi_h)) - j_{ps}(\frac{d\psi_h}{dt}, s_h) - j_c(\frac{d\psi_h}{dt}, s_h) - \ell_c(V(t), s_h) = 0, \quad (19)
$$

with $\psi_{0,h}$ a representation of $\psi_0$ in the space $V_h$. The bilinear forms appearing in (19) are defined as for (17) and evaluated for functions in the discrete spaces.
The space $P_{loc}(T)$ reads

$$P_{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 at $V_m, V_\ell$ with $m, \ell \in \{1, 2, 3\} \setminus \{i\}$ and at the barycenter $G$, for each $i = 1, 2, 3$. In the $P_{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$.

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}$. Let us define

$$\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 semi-discrete evolution problem is obtained by applying the mortar finite element approach to (17). The semi-discrete problem 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}_{\circ}^{\text{ex}} \times \mathcal{V}_{\circ}^{\text{in}}$, it holds

$$a(\psi_h, s_h) - J_p(\psi_h, s_h; f(\psi_h)) - j_{ps}(\frac{d\psi_h}{dt}, s_h) - j_c(\frac{d\psi_h}{dt}, s_h) - \ell_c(V(t), s_h) = 0, \quad (19)$$

with $\psi_{0,h}$ a representation of $\psi_0$ in the space $\mathcal{V}_h$. The bilinear forms appearing in (19) are defined as for (17) and evaluated for functions in the discrete spaces.

We rely on an implicit Euler scheme to fully discretize in time the semi-discrete problem. The discrete problem reads: Find $\psi_h^n = (\psi_h^{ex,n}, \psi_h^{in,n}) \in \mathcal{V}$ approximating $\psi_h(t_n)$, for $t_n \in [0, T]$, such that, $\forall s_h = (v_h, w_h) \in \mathcal{V}_{\circ}^{\text{ex}} \times \mathcal{V}_{\circ}^{\text{in}}$, it holds

$$a(\psi_h^n, s_h) - J_p(\psi_h^n, s_h; f(\psi_h^n, t_n)) - \frac{1}{\Delta t}j_{ps}(\psi_h^n, s_h) - \frac{1}{\Delta t}j_c(\psi_h^n, s_h) \\ + \frac{1}{\Delta t}j_{ps}(\psi_h^{n-1}, s_h) + \frac{1}{\Delta t}j_c(\psi_h^{n-1}, s_h) - \ell_c(V^n, s_h) = 0. \quad (20)$$
We rely on an implicit Euler scheme to fully discretize in time the semi-discrete problem. The discrete problem reads: Find $\psi_h^n = (\psi_h^{ex,n}, \psi_h^{in,n}) \in \mathcal{V}$ approximating $\psi_h(t_n)$, for $t_n \in [0, T]$, such that, $\forall s_h = (v_h, w_h) \in \mathcal{V}_o^{\text{ex}} \times \mathcal{V}_o^{\text{in}}$, it holds

$$a(\psi_h^n, s_h) - J_p(\psi_h^n, s_h; f(\psi_h^n, t_n)) - \frac{1}{\Delta t}j_{ps}(\psi_h^n, s_h) - \frac{1}{\Delta t}j_c(\psi_h^n, s_h) \\ + \frac{1}{\Delta t}j_{ps}(\psi_h^{n-1}, s_h) + \frac{1}{\Delta t}j_c(\psi_h^{n-1}, s_h) - \ell_c(V^n, s_h) = 0. \quad (20)$$

For the approximation in space of $f$ we use a decomposition in $k+1$ cubic splines [16], $\hat{g}_0, \dots, \hat{g}_k$, defined on $[0, 1]$ by nodes $x_0 = 0 < \dots < x_i < \dots < x_k = 1$ and $\hat{g}_i(x_j) = \delta_{ij}$. We introduce the normalized variable $\psi_N = \frac{\psi - \psi_a}{\psi_b - \psi_a}$ for $\psi \in [\psi_a, \psi_b]$ and define $g_i(\psi) = \hat{g}_i(\psi_N)$. We choose to approximate

$$f(\psi_h^n, t_n) \approx f_k^n(\psi_h^n) = \sum_{i=0}^{k} f_{i,n} g_i(\psi_h^n), \quad (21)$$

and we therefore have

$$f'(\psi_h^n, t_n) \approx [f']_k^n(\psi_h^n) = \sum_{i=0}^{k} f_{i,n} g'_i(\psi_h^n). \quad (22)$$

All spline functions are 0 at $\psi = \psi_b$ except the last one and therefore in order to impose the boundary condition $f_k^n(\psi_b) = f_b(t_n)$ we have $f_{k,n} = f_b(t_n)$. For the approximation in time of (4) we again rely on an implicit Euler scheme and follow the semi-Lagrangian discretization of $\frac{df(\psi,t)}{dt}$ proposed in [25]

$$\frac{df(\psi_h^n, t_n)}{dt} \approx \frac{f_k^n(\psi_h^n) - f_k^{n-1}(\psi_h^{n-1})}{\Delta t}. \quad (23)$$

The input profiles $p(\psi, t)$, $\eta^{\|}(\psi, t)$ and $\langle j_{ni} \cdot \mathbf{B} \rangle (\psi, t)$ are also decomposed in a function basis. They are evaluated implicitly in $\psi_h^n$ but their shape are only known at time $t_{n-1}$ since they are computed from the $T_e$ equation. We note this in the following way

$$p(\psi_h^n, t_{n-1}) \approx p_k^{n-1}(\psi_h^n) = \sum_{i=0}^{k} p_{i,n-1} g_i(\psi_h^n), \quad (24)$$

$$\eta^{\|}(\psi_h^n, t_{n-1}) \approx \eta_k^{n-1}(\psi_h^n) = \sum_{i=0}^{k} \eta_{i,n-1} g_i(\psi_h^n), \quad (25)$$

$$\langle j_{ni} \cdot \mathbf{B} \rangle (\psi_h^n, t_{n-1}) \approx [jB]_k^{n-1}(\psi_h^n) = \sum_{i=0}^{k} [jB]_{i,n-1} g_i(\psi_h^n). \quad (26)$$

The discrete resistive equation (4) reads for all $i = 0, \dots, k-1$:

$$\int_{\Omega_p(\psi_h^n)} \frac{f_k^n(\psi_h^n) - f_k^{n-1}(\psi_h^{n-1})}{\Delta t} \frac{g_i(\psi_h^n)}{r} dr dz \\ + \int_{\Omega_p(\psi_h^n)} \frac{\psi_h^n - \psi_h^{n-1}}{\Delta t} \frac{f_k^n(\psi_h^n)}{r} g'_i(\psi_h^n) dr dz \quad (27)$$

$$+ \int_{\Omega_p(\psi_h^n)} \eta_k^{n-1}(\psi_h^n) (r[p']_k^{n-1}(\psi_h^n) f_k^n(\psi_h^n) + \frac{[f_k^n]^2(\psi_h^n) + |\nabla \psi_h^n|^2}{\mu_0 r} [f']_k^n(\psi_h^n) \\ -r[jB]_k^{n-1}(\psi_h^n)) g'_i(\psi_h^n) dr dz = 0. \quad \text{Inria}$$
For the approximation in space of $f$ we use a decomposition in $k+1$ cubic splines [47], $\hat{g}_0, \dots, \hat{g}_k$, defined on $[0, 1]$ by nodes $x_0 = 0 < \dots < x_i < \dots < x_k = 1$ and $\hat{g}_i(x_j) = \delta_{ij}$. We introduce the normalized variable $\psi_N = \frac{\psi - \psi_a}{\psi_b - \psi_a}$ for $\psi \in [\psi_a, \psi_b]$ and define $g_i(\psi) = \hat{g}_i(\psi_N)$. We choose to approximate

$$f(\psi_h^n, t_n) \approx f_k^n(\psi_h^n) = \sum_{i=0}^{k} f_{i,n} g_i(\psi_h^n), \quad (21)$$

and we therefore have

$$f'(\psi_h^n, t_n) \approx [f']_k^n(\psi_h^n) = \sum_{i=0}^{k} f_{i,n} g'_i(\psi_h^n). \quad (22)$$

All spline functions are 0 at $\psi = \psi_b$ except the last one and therefore in order to impose the boundary condition $f_k^n(\psi_b) = f_b(t_n)$ we have $f_{k,n} = f_b(t_n)$. For the approximation in time of (4) we again rely on an implicit Euler scheme and follow the semi-Lagrangian discretization of $\frac{df(\psi,t)}{dt}$ proposed in [19]

$$\frac{df(\psi_h^n, t_n)}{dt} \approx \frac{f_k^n(\psi_h^n) - f_k^{n-1}(\psi_h^{n-1})}{\Delta t}. \quad (23)$$

The input profiles $p(\psi, t)$, $\eta^{\|}(\psi, t)$ and $\langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle (\psi, t)$ are also decomposed in a function basis. They are evaluated implicitly in $\psi_h^n$ but their shape are only known at time $t_{n-1}$ since they are computed from the $T_e$ equation. We note this in the following way

$$p(\psi_h^n, t_{n-1}) \approx p_k^{n-1}(\psi_h^n) = \sum_{i=0}^{k} p_{i,n-1} g_i(\psi_h^n), \quad (24)$$

$$\eta^{\|}(\psi_h^n, t_{n-1}) \approx \eta_k^{n-1}(\psi_h^n) = \sum_{i=0}^{k} \eta_{i,n-1} g_i(\psi_h^n), \quad (25)$$

$$\langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle (\psi_h^n, t_{n-1}) \approx [jB]_k^{n-1}(\psi_h^n) = \sum_{i=0}^{k} [jB]_{i,n-1} g_i(\psi_h^n). \quad (26)$$
For the computation of integrals over $\Omega^{\text{ex}}$ which appear in (20) we use a one-point barycentric quadrature rule. For the computation of integrals over $\Omega^{\text{in}}$ which appear in (20) we use a six-point quadrature rule of order 4 (see [37]) on each sub-triangle. This last rule is also used for the computation of integrals over $\Omega_p(\psi_h^n)$ from (27). A difficulty which arises with rHCT finite elements is the computation of the iso-contours and therefore of the plasma boundary $\Gamma_p(\psi_h)$. Contrary to what has been done in [26, 25] with piecewise linear finite elements, we do not compute explicitely the plasma boundary and do not consider a specific quadrature rule for elements $T$ such that the intersection $T \cap \Omega_p(\psi_h)$ is neither empty nor the whole triangle. We do not look for such elements $T$, but rather set the integrand to be 0 for points $(r, z)$ outside the plasma domain. Once the point defining the plama boundary $(r_b, z_b)$ and the magnetic axis $(r_a, z_a)$ are known, the plasma domain can be defined as $\{(r, z) \in \Omega, 0 \le \psi_{N,h}(r, z) \le 1 \text{ and } \nabla\psi_h(r, z) \cdot (r_a - r, z_a - z)^{\top} \ge 0\}$. Integrals over the plasma domain are then approximated with the same quadrature rule as for integrals over domains not depending on $\psi$. Since the plasma domain evolves with time it can happen that a quadrature point inside $\Omega_p(\psi_h^n)$ is not inside $\Omega_p(\psi_h^{n-1})$. In this case we take $f_k^{n-1}(\psi_h^{n-1}) = B_0 r_0$.

## 3.2 Matrix problem and Newton method for the $(\psi, f)$ system

Let us denote by $\{v_i^{\text{ex}}\}_{i=1,N^{\text{ex}}}$ the 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}}}$ those of $\mathcal{V}^{\text{in}}$ for the rHCT dofs at the $N^{\text{in}}$ nodes $V_j \in \tau^{\text{in}}$.

Let the vector $\psi$ of dimension $N^{\text{ex}}+3N^{\text{in}}$ gather all dofs of $\psi_h \in \mathcal{V}_h$. Let $\mathbf{u}^{\text{ex}}$ and $\mathbf{u}^{\text{in}}$ gather the dofs of $\psi_h^{\text{ex}} \in \mathcal{V}^{\text{ex}}$ and $\psi_h^{\text{in}} \in \mathcal{V}^{\text{in}}$. 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}}$) represent dofs in $V_o^{\text{ex}}$ (resp. $V_o^{\text{in}}$) and $V_\partial^{\text{ex}}$ (resp. $V_\partial^{\text{in}}$), $\mathbf{u}_o^{\text{ex}}$ is of dimension $N_\partial$ and $\mathbf{u}_\partial^{\text{in}}$ of dimension $3N_\partial$. The mortar 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$, 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$ and $k = 1, 3N_\partial$. The coupling condition is taken into account by introducing the reduced variable, $\mathbf{X}_\psi$, and matrix $\mathbf{Q}$ 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} & 0 & 0 \\ 0 & 0 & \mathbf{P}^{-1}\mathbf{D} \\ 0 & \mathbf{I} & 0 \\ 0 & 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}_\psi. $$

Let the vector $\mathbf{X}_f$ gather the $k$ free coefficients of $f_k$. The algebraic system of coupled non-linear equations associated with (20) and (27) to be solved at each time step can be written as

$$ \begin{aligned} \mathbf{e}(\mathbf{X}_\psi^n, \mathbf{X}_f^n; \mathbf{X}_\psi^{n-1}) &= 0_{\mathbb{R}^{N^{\text{ex}}+3N^{\text{in}}-N_\partial}} \\ \mathbf{d}(\mathbf{X}_\psi^n, \mathbf{X}_f^n; \mathbf{X}_\psi^{n-1}, \mathbf{X}_f^{n-1}) &= 0_{\mathbb{R}^k} \end{aligned} \quad (28) $$

where $\mathbf{d}$ represents the resistive diffusion equation for $f$ (component $\mathbf{d}_i$ is (27)) and $\mathbf{e}$ represents the equilibrium model. We have

$$ \begin{aligned} \mathbf{e}(\mathbf{X}_\psi^n, \mathbf{X}_f^n; \mathbf{X}_\psi^{n-1}) := \mathbf{Q}^T &[(\mathbf{A}_0 - \frac{1}{\Delta t}(\mathbf{M}_{ps} + \mathbf{L}^T \mathbf{R} \mathbf{L}))\mathbf{Q}\mathbf{X}_\psi^n \\ &+ \frac{1}{\Delta t}(\mathbf{M}_{ps} + \mathbf{L}^T \mathbf{R} \mathbf{L})\mathbf{Q}\mathbf{X}_\psi^{n-1} \\ &+ \mathbf{A}_\mu(\psi^n) - \mathbf{J}(\psi^n, \mathbf{X}_f^n) - \mathbf{LSV}^n] \end{aligned} \quad (29) $$
The discrete resistive equation (4) reads for all $i = 0, \dots, k-1$:

$$ 
\begin{aligned}
& \int_{\Omega_p(\psi_h^n)} \frac{f_k^n(\psi_h^n) - f_k^{n-1}(\psi_h^{n-1})}{\Delta t} \frac{g_i(\psi_h^n)}{r} dr dz \\
& + \int_{\Omega_p(\psi_h^n)} \frac{\psi_h^n - \psi_h^{n-1}}{\Delta t} \frac{f_k^n(\psi_h^n)}{r} g'_i(\psi_h^n) dr dz \\
& + \int_{\Omega_p(\psi_h^n)} \eta_k^{n-1}(\psi_h^n) (r[p']_k^{n-1}(\psi_h^n) f_k^n(\psi_h^n) + \frac{[f_k^n]^2(\psi_h^n) + |\nabla \psi_h^n|^2}{\mu_0 r}) [f']_k^n(\psi_h^n) \\
& - r[jB]_k^{n-1}(\psi_h^n)) g'_i(\psi_h^n) dr dz = 0. \qquad (27)
\end{aligned}
$$

For the computation of integrals over $\Omega^{\text{ex}}$ which appear in (20) we use a one-point barycentric quadrature rule. For the computation of integrals over $\Omega^{\text{in}}$ which appear in (20) we use a six-point quadrature rule of order 4 (see [48]) on each sub-triangle. This last rule is also used for the computation of integrals over $\Omega_p(\psi_h^n)$ from (27). A difficulty which arises with rHCT finite elements is the computation of the iso-contours and therefore of the plasma boundary $\Gamma_p(\psi_h)$. Contrary to what has been done in [28, 19] with piecewise linear finite elements, we do not compute explicitely the plasma boundary and do not consider a specific quadrature rule for elements $T$ such that the intersection $T \cap \Omega_p(\psi_h)$ is neither empty nor the whole triangle. We do not look for such elements $T$, but rather set the integrand to be 0 for points $(r, z)$ outside the plasma domain. Once the point defining the plasma boundary $(r_b, z_b)$ and the magnetic axis $(r_a, z_a)$ are known, the plasma domain can be defined as $\{(r, z) \in \Omega, 0 \le \psi_{N,h}(r, z) \le 1 \text{ and } \nabla\psi_h(r, z) \cdot (r_a - r, z_a - z)^{\top} \ge 0\}$. Integrals over the plasma domain are then approximated with the same quadrature rule as for integrals over domains not depending on $\psi$. Since the plasma domain evolves with time it can happen that a quadrature point inside $\Omega_p(\psi_h^n)$ is not inside $\Omega_p(\psi_h^{n-1})$. In this case we take $f_k^{n-1}(\psi_h^{n-1}) = B_0 r_0$.

### 3.2. Matrix problem and Newton method for the $(\psi, f)$ system

Let us denote by $\{v_i^{\text{ex}}\}_{i=1,N^{\text{ex}}}$ the 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}}}$ those of $\mathcal{V}^{\text{in}}$ for the rHCT dofs at the $N^{\text{in}}$ nodes $V_j \in \tau^{\text{in}}$.

Let the vector $\psi$ of dimension $N^{\text{ex}} + 3N^{\text{in}}$ gather all dofs of $\psi_h \in \mathcal{V}_h$. Let $\mathbf{u}^{\text{ex}}$ and $\mathbf{u}^{\text{in}}$ gather the dofs of $\psi_h^{\text{ex}} \in \mathcal{V}^{\text{ex}}$ and $\psi_h^{\text{in}} \in \mathcal{V}^{\text{in}}$. 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}}$) represent dofs in $V_o^{\text{ex}}$ (resp. $V_o^{\text{in}}$) and $V_\partial^{\text{ex}}$ (resp. $V_\partial^{\text{in}}$), $\mathbf{u}_o^{\text{ex}}$ is of dimension
where block matrix $A_0$ is associated with the linear part of $a(.,.)$,

$$A_0 = \begin{bmatrix} A_0^{\text{ex}} + C & 0 \\ 0 & A_0^{\text{in}} \end{bmatrix}$$

with

$$(A_0^{\text{ex}})_{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}},$$

$$(A_0^{\text{in}})_{ij} = \int_{\Omega^{\text{in}}} \frac{1}{\mu_0 r} \nabla w_i^{\text{in}} \cdot \nabla w_j^{\text{in}} dr dz, \quad i, j = 1, 3N^{\text{in}},$$

and C (resp. $M_{ps}$ and $L^T \text{RL}$) is the matrix associated with $c(.,.)$ (resp. $j_{ps}(.,.)$ and $j_c(.,.)$). Matrix $L^T \text{RL} = \begin{bmatrix} L^{\text{ex},T} \text{RL}^{\text{ex}} & 0 \\ 0 & 0 \end{bmatrix}$ is defined with $(L^{\text{ex}})_{ij} = \frac{1}{|\Omega_{c_i}|} \int_{\Omega_{c_i}} v_j^{\text{ex}} dr dz$ for $i = 1, N_c$ and $j = 1, N^{\text{ex}}$.

Vector $A_\mu(\psi^n) = \begin{bmatrix} A_\mu^{\text{ex}}(\psi^n) \\ 0 \end{bmatrix}$ comes from the non-linear part of $a(.,.)$ with

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

Vector $J(\psi^n, X_f^n) = \begin{bmatrix} 0 \\ J^{\text{in}}(\psi^n, X_f^n) \end{bmatrix}$ corresponds to the plasma current density term with $(J^{\text{in}}(\psi^n, X_f^n))_i = J_p(\psi_h^n, w_i; f_k^n)$ for $i = 1, 3N^{\text{in}}$.

System (28) is solved using Newton iterations which consists in updating iteratively $(X_\psi^{n,l}, X_f^{n,l})$ until convergence to $(X_\psi^n, X_f^n)$, starting from $(X_\psi^{n,0}, X_f^{n,0}) = (X_\psi^{n-1}, X_f^{n-1})$ and solving the linear system:

$$\begin{bmatrix} D_{X_\psi} \mathbf{e}(X_\psi^{n,l-1}, X_f^{n,l-1}) & D_{X_f} \mathbf{e}(X_\psi^{n,l-1}, X_f^{n,l-1}) \\ D_{X_\psi} \mathbf{d}(X_\psi^{n,l-1}, X_f^{n,l-1}) & D_{X_f} \mathbf{d}(X_\psi^{n,l-1}, X_f^{n,l-1}) \end{bmatrix} \begin{bmatrix} X_\psi^{n,l} - X_\psi^{n,l-1} \\ X_f^{n,l} - X_f^{n,l-1} \end{bmatrix} = - \begin{bmatrix} \mathbf{e}(X_\psi^{n,l-1}, X_f^{n,l-1}) \\ \mathbf{d}(X_\psi^{n,l-1}, X_f^{n,l-1}) \end{bmatrix}. \quad (30)$$

Convergence of the Newton method depends on the computation of the derivatives in (30). We have $D_{X_\psi} \mathbf{e}(X_\psi, X_f) = \mathbf{Q}^T [A_0 - \frac{1}{\Delta t} (M_{ps} + L^T \text{RL}) + D_\psi A_\mu(\psi) - D_\psi J(\psi, X_f)] \mathbf{Q}$ and $D_{X_f} \mathbf{e}(X_\psi, X_f) = -\mathbf{Q}^T D_{X_f} J(\psi, X_f)$. The derivative $D_\psi A_\mu(\psi)$ is known analytically [26]. The derivatives with respect to dofs of $f$, $D_{X_f} J(\psi, X_f)$ and $D_{X_f} \mathbf{d}(X_\psi, X_f)$ can be computed without too many difficulties. On the contrary derivatives with respect to $\psi$, $D_\psi J(\psi, X_f)$ and $D_{X_f} \mathbf{d}(X_\psi, X_f)$ require special care because of the dependence on $\psi$ of the plasma domain $\Omega_p(\psi)$ and of the basis function $g_i(\psi)$. Formulas for these derivatives are provided in C.

### 3.3 The $T_e$ equation

The $T_e$ equation is discretized using centered finite differences in space and the implicit Euler scheme in time. Let $T_{e,j}^n$ denote the approximation of $T_e(x_j, t_n)$ with $x_j = (j + 1/2)\Delta x$ for $j = 0, \dots, N - 1$, the scheme can be written for $j = 1, \dots, N - 2$ as

$$B_j^n T_{e,j}^n - C_j^n \frac{\Delta t}{\Delta x^2} (D_{j+1/2}^n T_{e,j+1}^n - (D_{j-1/2}^n + D_{j+1/2}^n) T_{e,j}^n + D_{j-1/2}^n T_{e,j-1}^n) \\ + C_j^n \frac{\Delta t}{2\Delta x} (E_{j+1}^n T_{e,j+1}^n - E_{j-1}^n T_{e,j-1}^n) + \Delta t F_j^n T_{e,j}^n = \Delta t G_j^{n-1,n} + B_j^{n-1} T_{e,j}^{n-1} \quad (31)$$

Inria
$N_\partial$ and $\mathbf{u}_\partial^{\text{in}}$ of dimension $3N_\partial$. The mortar 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$, 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$ and $k = 1, 3N_\partial$. The coupling condition is taken into account by introducing the reduced variable, $\mathbf{X}_\psi$, and matrix $\mathbf{Q}$ such that

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

Let the vector $\mathbf{X}_f$ gather the $k$ free coefficients of $f_k$. The algebraic system of coupled non-linear equations associated with (20) and (27) to be solved at each time step can be written as

$$ \begin{aligned} \mathbf{e}(\mathbf{X}_\psi^n, \mathbf{X}_f^n; \mathbf{X}_\psi^{n-1}) &= 0_{\mathbb{R}^{N^{\text{ex}}+3N^{\text{in}}-N_\partial}} \\ \mathbf{d}(\mathbf{X}_\psi^n, \mathbf{X}_f^n; \mathbf{X}_\psi^{n-1}, \mathbf{X}_f^{n-1}) &= 0_{\mathbb{R}^k} \end{aligned} \qquad (28) $$

where $\mathbf{d}$ represents the resistive diffusion equation for $f$ (component $\mathbf{d}_i$ is (27)) and $\mathbf{e}$ represents the equilibrium model. We have

$$ \begin{aligned} \mathbf{e}(\mathbf{X}_\psi^n, \mathbf{X}_f^n; \mathbf{X}_\psi^{n-1}) :={}& \mathbf{Q}^T \left[ \begin{aligned} & (\mathbf{A}_0 - \frac{1}{\Delta t} (\mathbf{M}_{ps} + \mathbf{L}^T \mathbf{R} \mathbf{L})) \mathbf{Q} \mathbf{X}_\psi^n \\ & + \frac{1}{\Delta t} (\mathbf{M}_{ps} + \mathbf{L}^T \mathbf{R} \mathbf{L}) \mathbf{Q} \mathbf{X}_\psi^{n-1} \\ & + \mathbf{A}_\mu(\psi^n) - \mathbf{J}(\psi^n, \mathbf{X}_f^n) - \mathbf{LSV}^n \end{aligned} \right] \quad (29) \end{aligned} $$

where block matrix $\mathbf{A}_0$ is associated with the linear part of $a(.,.)$,

$$ \mathbf{A}_0 = \begin{bmatrix} \mathbf{A}_0^{\text{ex}} + \mathbf{C} & 0 \\ 0 & \mathbf{A}_0^{\text{in}} \end{bmatrix} $$

with

$$ (\mathbf{A}_0^{\text{ex}})_{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}}, $$

$$ (\mathbf{A}_0^{\text{in}})_{ij} = \int_{\Omega^{\text{in}}} \frac{1}{\mu_0 r} \nabla w_i^{\text{in}} \cdot \nabla w_j^{\text{in}} dr dz, \quad i, j = 1, 3N^{\text{in}}, $$
and is modified for $j = 0$ to take account the fact that coefficients vanish at $x = 0$ since $V'(0) = 0$, and for $j = N - 1$ to take into account the Dirichlet boundary condition. All coefficients depend on the geometry of the plasma and can be taken implicit. Only the source term $G$ which depends on $\eta^{\parallel}$ and $<\mathbf{j}_{ni} \cdot \mathbf{B}>$ cannot be taken fully implicit and we denote this by $G_j^{n-1,n}$. At each time step the profile $T_e(x, t_n)$ is updated and knowing the current equilibrium solution ($\psi_h^n, f_k^n$), profiles $p_k^n(\psi_h^n), \eta_k^n(\psi_h^n), [jB]_k^n(\psi_h^n)$ can be updated.

# 4 Numerical results

In this section we present three test cases for the WEST Tokamak [11, 9, 10]. Case 1 is a VDE using parametric functions for $p', \eta^{\parallel}$ and $<\mathbf{j}_{ni} \cdot \mathbf{B}>$. Case 2 is the same VDE but using the $T_e$ equation in the model. This VDE corresponds to WEST shot 59695 and comparison with experimental results are used to validate the simulation. Case 3 is the formation of an X-point using a magnetic feedback controller. The test case shown in [25] has also been computed and the results are presented in A.

## 4.1 Simulation setup

All numerical methods used in this paper are implemented in the code NICE [20]. Although there is no restriction in the method to deal with any tokamak, we concentrate in this work on the simulation of the equilibrium evolution of plasmas for the WEST tokamak. The starting point for such a simulation is always the computation of a first equilibrium configuration. This equilibrium is computed using the so-called static inverse mode of NICE: given input Grad-Shafranov profiles, here taken as parametric functions $p'(\psi) = \lambda\beta(1 - \psi_N^{\alpha_p})^{\gamma_p}$ and $ff'(\psi) = \lambda(1 - \beta)(1 - \psi_N^{\alpha_f})^{\gamma_f}$, a plasma current $I_p$ and a desired plasma shape (a set of $(r, z)$ points describing the desired plasma boundary), we solve a PDE-constrained optimization problem to find the currents $I_{c_i}$ in the poloidal field coils for which the solution $\psi$ to the direct static equilibrium problem

$$ -\Delta^* \psi = \begin{cases} rp'(\psi, t) + \frac{1}{\mu_0 r} ff'(\psi, t) & \text{in } \Omega_p(\psi), \\ \frac{1}{|\Omega_{c_i}|} I_i & \text{in } \Omega_{c_i}, \quad i = 1, \dots, N_c, \\ 0 & \text{elsewhere,} \end{cases} \qquad (32) $$

shows a plasma boundary close to the desired one.

This inverse equilibrium computation is done using only P1 FEM. We recompute the equilibrium solving (32) using the currents provided by the inverse solution but using the coupled P1-rHCT FEM method. The solution $\psi$ is taken as the initial $\psi_0$ for evolution simulations.

Integrating $ff'(\psi)$ we compute the initial $f_0(\psi_0)$ profile. Integrating $p'(\psi)$ the initial pressure profile is computed and assuming a given electron density profile $n_e$ (see Fig. 1a), the initial electron temperature profile $T_{e,0}$ (see Fig. 1b) is constructed. The electron density profile is represented using the following parametric expression:

$$ n_e(x) = (\gamma - \beta)(1 - \alpha(x - x_0)) \frac{\exp(-\frac{x-x_0}{h})}{\exp(-\frac{x-x_0}{h}) + \exp(-\frac{x-x_0}{h})} + \beta, \qquad (33) $$

with $x_0 = 0.9$, $h = 0.2$, $\alpha = 0.5$, $\beta = 1.1 \times 10^{17}$, $\gamma = 3.0 \times 10^{19}$. Parameters for functions $p'(\psi)$ and $ff'(\psi)$ are: $\beta = 0.54537755$, $\alpha_p = \alpha_f = 1.0$, $\gamma_p = \gamma_f = 1.09075510$. The total plasma current is $I_p = 482000$ [A] and the toroidal magnetic field is $B_0 = 3.524$ [T].
and C (resp. $M_{ps}$ and $L^T R L$) is the matrix associated with $c(.,.)$ (resp. $j_{ps}(.,.)$ and $j_c(.,.)$). Matrix $L^T R L = \begin{bmatrix} L^{\text{ex},T} R L^{\text{ex}} & 0 \\ 0 & 0 \end{bmatrix}$ is defined with $(L^{\text{ex}})_{ij} = \frac{1}{|\Omega_{c_i}|} \int_{\Omega_{c_i}} v_j^{\text{ex}} dr dz$ for $i=1, N_c$ and $j=1, N^{\text{ex}}$.

390

Vector $A_\mu(\psi^n) = \begin{bmatrix} A_\mu^{\text{ex}}(\psi^n) \\ 0 \end{bmatrix}$ comes from the non-linear part of $a(.,.)$ with

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

392

Vector $J(\psi^n, X_f^n) = \begin{bmatrix} 0 \\ J^{\text{in}}(\psi^n, X_f^n) \end{bmatrix}$ corresponds to the plasma current density term with $(J^{\text{in}}(\psi^n, X_f^n))_i = J_p(\psi_h^n, w_i; f_k^n)$ for $i=1, 3N^{\text{in}}$.

394

System (28) is solved using Newton iterations which consists in updating iteratively $X^{n,l} = (X_{\psi}^{n,l}, X_{f}^{n,l})$ until convergence to $(X_{\psi}^{n}, X_{f}^{n})$, starting from $(X_{\psi}^{n,0}, X_{f}^{n,0}) = (X_{\psi}^{n-1}, X_{f}^{n-1})$ and solving the linear system:

396

$$ \begin{bmatrix} D_{X_{\psi}} \mathbf{e}(\mathbf{X}_{\psi}^{n,l-1}, \mathbf{X}_{f}^{n,l-1}) & D_{X_{f}} \mathbf{e}(\mathbf{X}_{\psi}^{n,l-1}, \mathbf{X}_{f}^{n,l-1}) \\ D_{X_{\psi}} \mathbf{d}(\mathbf{X}_{\psi}^{n,l-1}, \mathbf{X}_{f}^{n,l-1}) & D_{X_{f}} \mathbf{d}(\mathbf{X}_{\psi}^{n,l-1}, \mathbf{X}_{f}^{n,l-1}) \end{bmatrix} \begin{bmatrix} \mathbf{X}_{\psi}^{n,l} - \mathbf{X}_{\psi}^{n,l-1} \\ \mathbf{X}_{f}^{n,l} - \mathbf{X}_{f}^{n,l-1} \end{bmatrix} = - \begin{bmatrix} \mathbf{e}(\mathbf{X}_{\psi}^{n,l-1}, \mathbf{X}_{f}^{n,l-1}) \\ \mathbf{d}(\mathbf{X}_{\psi}^{n,l-1}, \mathbf{X}_{f}^{n,l-1}) \end{bmatrix}. \quad (30) $$

398

The stopping criteria is given by a tolerance $\epsilon_{tol}$ on the relative error $||X^{n,l} - X^{n,l-1}||/||X^{n,l-1}|| < \epsilon_{tol}$. Convergence of the Newton method depends on the computation of the derivatives in (30). We have $D_{X_{\psi}} \mathbf{e}(\mathbf{X}_{\psi}, \mathbf{X}_{f}) =$

400

$Q^T[A_0 - \frac{1}{\Delta t}(M_{ps} + L^T R L) + D_{\psi}A_{\mu}(\psi) - D_{\psi}J(\psi, \mathbf{X}_{f})]Q$ and $D_{X_{f}} \mathbf{e}(\mathbf{X}_{\psi}, \mathbf{X}_{f}) =$

402

$-Q^T D_{X_{f}} J(\psi, \mathbf{X}_{f})$. The derivative $D_{\psi}A_{\mu}(\psi)$ is known analytically [28]. The derivatives with respect to dofs of $f$, $D_{X_{f}} J(\psi, \mathbf{X}_{f})$ and $D_{X_{f}} \mathbf{d}(\mathbf{X}_{\psi}, \mathbf{X}_{f})$ can be computed without too many difficulties. On the contrary derivatives with respect to $\psi$, $D_{\psi}J(\psi, \mathbf{X}_{f})$ and $D_{X_{\psi}} \mathbf{d}(\mathbf{X}_{\psi}, \mathbf{X}_{f})$ require special care because of the dependence on $\psi$ of the plasma domain $\Omega_{\psi}(\psi)$ and of the basis function $g_i(\psi)$. Formulas for these derivatives are provided in Appendix C.

404

406

### 3.3. The $T_e$ equation

408

The $T_e$ equation is discretized using centered finite differences in space and the implicit Euler scheme in time. Let $T_{e,j}^n$ denote the approximation
<figure><img src="image_4.png" /><figcaption>Figure 1: Electron density profile and initial temperature profile</figcaption></figure>

The mesh contains 15868 nodes (30815 triangles) in the external domain and 3612 nodes (6683 triangles) in the internal domain. The different domains are represented on Fig. 2 together with the equilibrium computed using the mortar method.

The transport solver is initialized from the $\psi_0$ equilibrium and computes the values of $T_e^0(\rho^0)$, $p_k^0(\psi_h^0)$, $\eta_k^0(\psi_h^0)$, $[jB]^0_k(\psi_h^0)$.

After this initialization phase, the time stepping goes as follows. At a time $n$, the transport model provides profiles $p_k^n(\psi_h^n)$, $\eta_k^n(\psi_h^n)$, $[jB]^n_k(\psi_h^n)$ needed by the equilibrium solver, this one interpolates the profiles on the chosen decomposition basis then computes the next equilibrium. Once the equilibrium is computed, the data are sent to the transport solver to compute the new profiles $p_k^{(n+1)}(\psi_h^{(n+1)})$, $\eta_k^{(n+1)}(\psi_h^{(n+1)})$, $[jB]^{(n+1)}_k(\psi_h^{(n+1)})$.

All the profiles computed from the 2D equilibrium and needed in the $T_e$ equation, such as $<1/r^2>$ for example, depend on the choice of the discrete $\psi_N$ values for which isocontours of $\psi$ are computed and used in the average computations. In order to restrain numerical errors near the magnetic axis that is to say near $x=0$ in the resolution of the $T_e$ equation, a good choice for these discrete $\psi_N$ values is to have them related to regularly spaced discrete $\rho_N$ values as shown on Fig. 3. The choice of these $\psi_N$ values changes at each time step since $\rho(\psi)$ evolves.

Although it has not appeared necessary in our numerical experiments it is possible to iterate between the equilibrium and the transport solver to replay the time step with these new profiles and ensure accurate resolution of the non-linearities. The simulation algorithm is summarized as follows in Algorithm 1.
of $T_e(x_j, t_n)$ with $x_j = (j + 1/2)\Delta x$ for $j = 0, \dots, N - 1$, the scheme can be written for $j = 1, \dots, N - 2$ as

$$B_j^n T_{e,j}^n - C_j^n \frac{\Delta t}{\Delta x^2} (D_{j+1/2}^n T_{e,j+1}^n - (D_{j-1/2}^n + D_{j+1/2}^n) T_{e,j}^n + D_{j-1/2}^n T_{e,j-1}^n) \\ + C_j^n \frac{\Delta t}{2\Delta x} (E_{j+1}^n T_{e,j+1}^n - E_{j-1}^n T_{e,j-1}^n) + \Delta t F_j^n T_{e,j}^n = \Delta t G_j^{n-1,n} + B_j^{n-1} T_{e,j}^{n-1} \quad (31)$$

and is modified for $j = 0$ to take account the fact that coefficients vanish at $x = 0$ since $V'(0) = 0$, and for $j = N - 1$ to take into account the Dirichlet boundary condition. All coefficients depend on the geometry of the plasma and can be taken implicit. Only the source term $G$ which depends on $\eta^{\|}$ and $\langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle$ cannot be taken fully implicit and we denote this by $G_j^{n-1,n}$. At each time step the profile $T_e(x, t_n)$ is updated and knowing the current equilibrium solution $(\psi_h^n, f_k^n)$, profiles $p_k^n(\psi_h^n)$, $\eta_k^n(\psi_h^n)$, $[jB]_k^n(\psi_h^n)$ can be updated.

## 4. Numerical results

In this section we present three test cases for the WEST Tokamak [24, 25, 26]. Case 1 is a free displacement of the plasma using parametric functions for $p', \eta^{\|}$ and $\langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle$. Case 2 is the same free displacement of the plasma as in case 1 but using the $T_e$ equation in the model. Both cases are presented in section 4.2 which ends with results from a similar simulation but in which the resistivity is set to 0. This free displacement of the plasma corresponds to WEST shot 59695 and comparison with experimental results are used to validate the simulation. Case 3 is the formation of an X-point using a magnetic feedback controller. It is presented in section 4.3 as well as results for a VDE obtained by turning off the vertical position magnetic control. The test case shown in [19] has also been computed and the results are presented in Appendix A.

### 4.1. Simulation setup

All numerical methods used in this paper are implemented in the code NICE [23]. Although there is no restriction in the method to deal with any tokamak, we concentrate in this work on the simulation of the equilibrium evolution of plasmas for the WEST tokamak. The starting point for such a simulation is always the computation of a first equilibrium configuration. This equilibrium is computed using the so-called static inverse mode of NICE:
<figure><img src="image_5.png" /><figcaption>Figure 2: Left: WEST computational domain geometry. The ABB domain boundary and the boundary $\mathcal{I}$ of domain $\Omega_{in}$ are shown in red. The poloidal field coils (named A, Bb, Bh, ...) are in blue. The iron domain is in green and passive structures in black. Right: Zoom on an equilibrium computed using the mortar coupling method. The plasma boundary is shown in red and the boundary $\mathcal{I}$ of domain $\Omega_{in}$ is in blue.</figcaption></figure>
given in input Grad-Shafranov profiles (here taken as parametric functions given in details lower in the text), a plasma current $I_p$ and a desired plasma shape (a set of $(r, z)$ points describing the desired plasma boundary), we solve a PDE-constrained optimization problem to find the currents $I_{c_i}$ in the poloidal field coils for which the solution $\psi$ to the direct static equilibrium problem

$$ \mathcal{L}\psi = \begin{cases} rp'(\psi) + \frac{1}{\mu_0 r} f f'(\psi) & \text{in } \Omega_p(\psi), \\ \frac{1}{|\Omega_{c_i}|} I_i & \text{in } \Omega_{c_i}, \quad i = 1, \dots, N_c, \\ 0 & \text{elsewhere,} \end{cases} \qquad (32) $$

shows a plasma boundary close to the desired one.

This inverse equilibrium computation is done using only P1 FEM. We recompute the equilibrium solving (32) using the currents provided by the inverse solution but using the coupled P1-rHCT FEM method. The solution $\psi$ is taken as the initial $\psi_0$ for evolution simulations.

Integrating $f f'(\psi)$ we compute the initial $f_0(\psi_0)$ profile. Integrating $p'(\psi)$ the initial pressure profile is computed and assuming a given electron density profile $n_e$ (see Fig. 1a), the initial electron temperature profile $T_{e,0}$ (see Fig. 1b) is constructed. The electron density profile is represented using the following parametric expression:

$$ n_e(x) = (c-b)(1-a(x-x_0)) \frac{\exp(-\frac{x-x_0}{h})}{\exp(-\frac{x-x_0}{h}) + \exp(\frac{x-x_0}{h})} + b, \qquad (33) $$

where $x = \rho/\rho_b$ is the normalized radial coordinate and $x_0 = 0.9$, $h = 0.2$, $a = 0.5$, $b = 1.1 \times 10^{17} [m^{-3}]$, $c = 3.0 \times 10^{19} [m^{-3}]$.

Parameters for functions $p'(\psi) = \lambda\beta(1 - \psi_N^{\alpha_p})^{\gamma_p}$ and $f f'(\psi) = \lambda(1 - \beta)(1 - \psi_N^{\alpha_f})^{\gamma_f}$ are: $\beta = 0.54537755$, $\alpha_p = \alpha_f = 1.0$, $\gamma_p = \gamma_f = 1.09075510$. The total plasma current is $I_p = 482000 [A]$ and the toroidal magnetic field is $B_0 = 3.524 [T]$.

The mesh contains 15868 nodes (30815 triangles) in the external domain and 3612 nodes (6683 triangles) in the internal domain. The different domains are represented on Fig. 2 together with the equilibrium computed using the mortar method.

The transport solver is initialized from the $\psi_0$ equilibrium and computes the values of $T_e^0(\rho^0)$, $p_k^0(\psi_h^0)$, $\eta_k^0(\psi_h^0)$, $[jB]_k^0(\psi_h^0)$.
<figure><img src="image_6.png" /><figcaption>Figure 3: Typical normalized toroidal flux coordinate $\rho_N$ with respect to $\psi_N$</figcaption></figure>

**Algorithm 1** NICE with resistive diffusion and transport using P1-rHCT FEM. The computational equilibrium solver is noted EquiSolver, the transport solver is TranspSolver, the data structure containing an equilibrium is dataEqui and the one containing the outputs from transport computations is dataTransp.

**Initialization**

* Compute static inverse P1

* Compute static direct P1-rHCT

► dataEqui ← EquiSolver

► TranspSolver ← dataEqui

**Time Stepping**

$n \leftarrow 0$

**while** $n < n_{\text{timestep}}$ **do**

rep $\leftarrow 0$

**for** rep $<$ replay **do**

▷ Either number of replay or convergence criterion

► dataTransp ← TranspSolver

► EquiSolver ← dataTransp

* Compute OneTimeStep EquiSolver P1-rHCT

► dataEqui ← EquiSolver

► TranspSolver ← dataEqui

* Compute OneTimeStep TranspSolver

**end for**

**end while**
<figure><img src="image_2.png" /><figcaption>Figure 1: Electron density profile and initial temperature profile</figcaption></figure>

After this initialization phase, the time stepping goes as follows. At a time $t_n$, the transport model provides profiles $p_k^n(\psi_h^n)$, $\eta_k^n(\psi_h^n)$, $[jB]^n_k(\psi_h^n)$ needed by the equilibrium solver, this one interpolates the profiles on the chosen decomposition basis then computes the next equilibrium. Once the equilibrium is computed, the data are sent to the transport solver to compute the new profiles $p_k^{n+1}(\psi_h^{n+1})$, $\eta_k^{n+1}(\psi_h^{n+1})$, $[jB]^n_k(\psi_h^{n+1})$.

All the profiles computed from the 2D equilibrium and needed in the $T_e$ equation, such as $<1/r^2>$ for example, depend on the choice of the discrete $\psi_N$ values for which isocontours of $\psi$ are computed and used in the average computations. In order to restrain numerical errors near the magnetic axis that is to say near $x=0$ in the resolution of the $T_e$ equation, a good choice for these discrete $\psi_N$ values is to have them related to regularly spaced discrete $\rho_N$ values as shown on Fig. 3. The choice of these $\psi_N$ values changes at each time step since $\rho(\psi)$ evolves.

Although it has not appeared necessary in our numerical experiments it is possible to iterate between the equilibrium and the transport solver to replay the time step with these new profiles and ensure accurate resolution of the non-linearities. The simulation algorithm is summarized as follows in Algorithm 1.
## 4.2 Free displacement of the plasma. Test case 1 and 2.

For the first two cases, the initial condition of the plasma is set up using the currents and plasma boundary from WEST shot 59695 at time 5.502[s] where the magnetic control is turned off. As a first numerical experiment, test case 1 is conducted to check the good behavior of the algorithm. Starting from the initial condition described above we let the plasma evolve freely. Voltages in the poloidal field coils circuits are constant equal to 0 during the time stepping. The time step is set to $\Delta t = 10^{-3}$ [s]. The $T_e$ evolution equation is not included at first and the resistivity and non-inductive current are held constant in time. We take $\langle j_{ni} \cdot B \rangle = 0$ for the non-inductive current term and use the analytic expression $\eta^{\parallel}(\psi_N) = \frac{10^{-8}}{1.1-\psi_N}$ for resistivity. The toroidal magnetic field $B_0$ is constant and so is the boundary condition $f_b(t) = B_0 r_0$. The $p'$ profile is also kept fixed during the simulation.

An important step in the development of the Newton method (30) is to check the accuracy of the computed derivatives presented in C. This is done by verifying that the first order Taylor expansion

$$ \|\mathbf{F}(\mathbf{X} + \varepsilon\mathbf{h}) - \mathbf{F}(\mathbf{X}) - \varepsilon D_X \mathbf{F}(\mathbf{X})\mathbf{h}\| = O(\varepsilon^2) $$

holds for the non-linear operators as illustrated in Fig. 4.

<figure><img src="image_7.png" /><figcaption>Figure 4: Example of a gradient test for derivatives $D_{X_\psi}e$, $D_{X_f}e$, $D_{X_\psi}d$, $D_{X_f}d$</figcaption></figure>

At each time step the Newton method converges in 6 to 12 iterations with a given stopping criterion of $10^{-10}$. Convergence history examples are given in Table 1 for test case 1 and Table 2 for case 2.

Figure 5 shows the results obtained with the developed coupled P1-rHCT FEM method using cubic splines for the discretization of $f$ with 10 degrees of freedom. The plasma moves upwards from the initial X-point configuration to a limiter configuration. Function $f$ evolves smoothly during the time stepping. The simulation fails to converge after 40 time steps while the plasma shrinks against the limiter.
<figure><img src="image_3.png" /><figcaption>Figure 2: Left: WEST computational domain geometry. The ABB domain boundary and the boundary $\mathcal{I}$ of domain $\Omega_{in}$ are shown in red. The poloidal field coils (named A, Bb, Bh, ...) are in blue. The iron domain is in green and passive structures in black. The limiter contour appears as a black dashed line. Right: Zoom on an equilibrium computed using the mortar coupling method. The plasma boundary is shown in red and the boundary $\mathcal{I}$ of domain $\Omega_{in}$ is in blue.</figcaption></figure>
<table><caption>Table 1: Convergence history of Newton iterations for case 1 : iteration number *n* and residual relative error ||**X**<sup>*n*</sup> − **X**<sup>*n*−1</sup>||/||**X**<sup>*n*−1</sup>||.</caption><thead><tr><th>n</th><th>t = 0.01</th><th>t = 0.02</th><th>t = 0.03</th></tr></thead><tbody><tr><td>1</td><td>1.112 × 10<sup>-3</sup></td><td>1.009 × 10<sup>-3</sup></td><td>9.493 × 10<sup>-4</sup></td></tr><tr><td>2</td><td>3.234 × 10<sup>-4</sup></td><td>4.217 × 10<sup>-4</sup></td><td>6.127 × 10<sup>-4</sup></td></tr><tr><td>3</td><td>1.817 × 10<sup>-4</sup></td><td>2.486 × 10<sup>-4</sup></td><td>4.444 × 10<sup>-4</sup></td></tr><tr><td>4</td><td>9.909 × 10<sup>-5</sup></td><td>1.513 × 10<sup>-4</sup></td><td>3.357 × 10<sup>-4</sup></td></tr><tr><td>5</td><td>7.580 × 10<sup>-5</sup></td><td>7.409 × 10<sup>-5</sup></td><td>1.806 × 10<sup>-4</sup></td></tr><tr><td>6</td><td>3.695 × 10<sup>-5</sup></td><td>4.944 × 10<sup>-5</sup></td><td>1.605 × 10<sup>-4</sup></td></tr><tr><td>7</td><td>1.659 × 10<sup>-5</sup></td><td>2.812 × 10<sup>-5</sup></td><td>8.099 × 10<sup>-5</sup></td></tr><tr><td>8</td><td>6.390 × 10<sup>-6</sup></td><td>1.499 × 10<sup>-5</sup></td><td>5.770 × 10<sup>-5</sup></td></tr><tr><td>9</td><td>4.341 × 10<sup>-6</sup></td><td>3.641 × 10<sup>-9</sup></td><td>2.687 × 10<sup>-5</sup></td></tr><tr><td>10</td><td>7.052 × 10<sup>-6</sup></td><td>3.202 × 10<sup>-14</sup></td><td>1.854 × 10<sup>-8</sup></td></tr><tr><td>11</td><td>8.087 × 10<sup>-10</sup></td><td></td><td>1.726 × 10<sup>-13</sup></td></tr><tr><td>12</td><td>2.739 × 10<sup>-14</sup></td><td></td><td></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 5: Case 1. Left : initial (blue) and final (red) plasma boundary. Right : *f* profile at different time steps.</figcaption></figure>

Let us now introduce in the simulation the electron temperature equation in order to have an evolving pressure profile as well as more realistic <**j**<sub>ni</sub> ⋅ **B**> and η|| profiles. The space discretization is chosen to be Δx = 2.0 × 10<sup>-3</sup>.

The source auxiliary $Q_a$ delivers a power $P_a = 0$ [MW] as it is the case in the experiment. Sources $Q_{ei}$ and $Q_{rad}$ are given in B. In the study case, the plasma is modeled with an effective atomic number $Z_{eff} = 3$. Moreover, the plasmas created in WEST have tungsten impurities which are responsible for a great part of the radiation losses. In order to simulate this effect, impurities with $Z_W = 74$ are added as explained in B. Since we don't use any model for the electron density evolution, we choose $\chi_{ne} = 0$ leading to $Q_i = 0$. The coefficient $c_{ie}$ is computed as explained in B. Concerning the remaining parameters we choose $r_{ie} = 1.0$ for the ion-electron temperature ratio (leading to $Q_{ei} = 0$) and choose a diffusion coefficient profile $\chi_{Te}$ proportional
**Algorithm 1** NICE with resistive diffusion and transport using P1-rHCT FEM. The computational equilibrium solver is noted EquiSolver, the transport solver is TranspSolver, the data structure containing an equilibrium is dataEqui and the one containing the outputs from transport computations is dataTransp.

**Initialization**

• Compute static inverse P1

• Compute static direct P1-rHCT

▶ dataEqui ← EquiSolver

▶ TranspSolver ← dataEqui

**Time Stepping**

$n \leftarrow 0$

**while** $n < n_{timestep}$ **do**

rep $\leftarrow 0$

**for** rep < replay **do**

▷ Either number of replay or convergence

criterion

▶ dataTransp ← TranspSolver

▶ EquiSolver ← dataTransp

• Compute OneTimeStep EquiSolver P1-rHCT

▶ dataEqui ← EquiSolver

▶ TranspSolver ← dataEqui

• Compute OneTimeStep TranspSolver

**end for**

**end while**
to the safety factor $q$ ($\chi_{T_e} = kq$ with $k = 0.15$ in all numerical experiments presented in this work).

<table><caption>Table 2: Convergence history of Newton iterations for case 2: iteration number $n$ and residual relative error $||\mathbf{X}^n - \mathbf{X}^{n-1}||/||\mathbf{X}^{n-1}||$.</caption><thead><tr><th>n</th><th>t = 0.01</th><th>t = 0.02</th><th>t = 0.03</th></tr></thead><tbody><tr><td>1</td><td>1.503 × 10<sup>-3</sup></td><td>1.230 × 10<sup>-3</sup></td><td>1.327 × 10<sup>-3</sup></td></tr><tr><td>2</td><td>2.526 × 10<sup>-4</sup></td><td>1.797 × 10<sup>-4</sup></td><td>2.359 × 10<sup>-4</sup></td></tr><tr><td>3</td><td>4.650 × 10<sup>-5</sup></td><td>3.376 × 10<sup>-5</sup></td><td>2.696 × 10<sup>-5</sup></td></tr><tr><td>4</td><td>1.283 × 10<sup>-5</sup></td><td>6.966 × 10<sup>-6</sup></td><td>9.006 × 10<sup>-6</sup></td></tr><tr><td>5</td><td>1.735 × 10<sup>-6</sup></td><td>1.608 × 10<sup>-9</sup></td><td>1.851 × 10<sup>-9</sup></td></tr><tr><td>6</td><td>5.714 × 10<sup>-11</sup></td><td>2.523 × 10<sup>-14</sup></td><td>3.029 × 10<sup>-14</sup></td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Figure 6: Case 2. Left : initial (blue) and final (red) plasma boundary. Right : $f$ profile at different time steps.</figcaption></figure>

Changes in the internal characteristics of the plasma, non-linear interplay between the different quantities and the introduction of more realistic bootstrap current and resistivity profiles, seem to improve convergence of the Newton method (Tab. 2) and to keep the algorithm converging over more time steps. The simulation holds for more than 100 [ms] (Fig. 6).

As in the previous case the plasma moves upwards, switches from an X-point configuration to a limiter configuration and finally shrinks against the limiter. From Figure 7 it can be seen that the experimental plasma shape evolution from WEST shot 59695 is quite well reproduced.

Different simulated profiles at several time steps during the simulation are shown on Fig. 8 and 9. Temperature drops in the plasma center and increases near the edge as time goes on (Fig. 8a). This change in temperature leads to a greater pressure on the edge, a decrease in the resistivity and an increase in the bootstrap current on the edge. Function $f$ globally increases. The Ohmic source is directly linked to the derivatives of $f$, the effect of this is a spike around $\rho_N = 0.9$.

In order to test the good numerical behavior of the code, multiple runs of the simulation are done varying the space and time discretization, but also the Newton stopping criterion or the
<figure><img src="image_4.png" /><figcaption>Figure 3: Typical normalized toroidal flux coordinate $\rho_N$ with respect to $\psi_N$</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 7: Comparison between the simulation (red) and shot 59695 (blue) for the initial equilibrium (left plot), after 100 time steps (middle plot) and the evolution of the geometric axis coordinates (right plots).</figcaption></figure>

even the replay.

Varying the space and time discretization does not change significantly the solution. An example (Fig. 10 and Fig. 11) shows the differences between the simulation above and a simulation using a time step divided by 2. The results are similar if we refine the space discretization in the transport solver or if we refine the FEM mesh.

Moreover the Newton stopping criterion can be set to $10^{-5}$ instead of $10^{-10}$ without changing the result substantially. Adding the replay loop between the equilibrium solver and the transport solver does not change significantly the solution either.
## 4.2. Free displacement of the plasma

Test case 1 and 2. For these test cases, the initial condition of the plasma is set up using the currents and plasma boundary from WEST shot 59695 at time 5.502[s]. This shot is dedicated to experience a free displacement of the plasma. After a controlled phase where the plasma is stabilized in an X-point configuration, the magnetic control is suddenly turned off at time 5.502[s] and the voltages in the poloidal field coils circuits set to 0.

As a first numerical experiment, test case 1 is conducted to check the good behavior of the algorithm. Starting from the initial condition described above we let the plasma evolve freely. Voltages in the poloidal field coils circuits are constant equal to 0 during the time stepping. The time step is set to $\Delta t = 10^{-3}$ [s]. The $T_e$ evolution equation is not included at first and the resistivity and non-inductive current are held constant in time. We take $\langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle = 0$ for the non-inductive current term and use the analytic expression $\eta^{\parallel}(\psi_N) = \frac{10^{-8}}{1.1-\psi_N}$ for resistivity. The toroidal magnetic field $B_0$ is constant and so is the boundary condition $f_b(t) = B_0 r_0$. The $p'$ profile is also kept fixed during the simulation.

An important step in the development of the Newton method (30) is to check the accuracy of the computed derivatives presented in Appendix C. This is done by verifying that the first order Taylor expansion

$$||\mathbf{F}(\mathbf{X} + \varepsilon\mathbf{h}) - \mathbf{F}(\mathbf{X}) - \varepsilon D_X \mathbf{F}(\mathbf{X})\mathbf{h}|| = O(\varepsilon^2)$$

holds for the non-linear operators as illustrated in Fig. 4.

At each time step the Newton method converges in 6 to 12 iterations with a given relative tolerance $\varepsilon_{tol} = 10^{-10}$ for the stopping criteria. Convergence history examples are given in Table 1 for test case 1 and Table 2 for case 2.

Figure 5 shows the results obtained with the developed coupled P1-rHCT FEM method using cubic splines for the discretization of $f$ with 10 degrees of freedom. The plasma moves upwards from the initial X-point configuration to a limiter configuration. Function $f$ evolves smoothly during the time stepping. The simulation fails to converge after 40 time steps while the plasma shrinks against the limiter.

Let us now introduce in the simulation the electron temperature equation in order to have an evolving pressure profile as well as more realistic $\langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle$ and $\eta^{\parallel}$ profiles. The space discretization is chosen to be $\Delta x = 2.0 \times 10^{-3}$.

The source auxiliary $Q_a$ delivers a power $P_a = 0$ [MW] as it is the case in the experiment. Sources $Q_{ei}$ and $Q_{rad}$ are given in Appendix B. In the
<figure><img src="image_11.png" /><figcaption>Figure 8: Case 2. Electron temperature, plasma pressure, resistivity and bootstrap current profiles at several time steps during the simulation</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9: Case 2. Source profiles at several time steps during the simulation</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Example of a gradient test for derivatives $D_{X_\psi}e, D_{X_f}e, D_{X_\psi}d, D_{X_f}d$</figcaption></figure>

<table><caption>Table 1: Convergence history of Newton iterations for case 1 : iteration number *n* and residual relative error ||**X**<sup>*n*</sup> - **X**<sup>*n*−1</sup>||/||**X**<sup>*n*−1</sup>||.</caption><thead><tr><th>n</th><th>t = 0.01</th><th>t = 0.02</th><th>t = 0.03</th></tr></thead><tbody><tr><td>1</td><td>1.112 × 10<sup>-3</sup></td><td>1.009 × 10<sup>-3</sup></td><td>9.493 × 10<sup>-4</sup></td></tr><tr><td>2</td><td>3.234 × 10<sup>-4</sup></td><td>4.217 × 10<sup>-4</sup></td><td>6.127 × 10<sup>-4</sup></td></tr><tr><td>3</td><td>1.817 × 10<sup>-4</sup></td><td>2.486 × 10<sup>-4</sup></td><td>4.444 × 10<sup>-4</sup></td></tr><tr><td>4</td><td>9.909 × 10<sup>-5</sup></td><td>1.513 × 10<sup>-4</sup></td><td>3.357 × 10<sup>-4</sup></td></tr><tr><td>5</td><td>7.580 × 10<sup>-5</sup></td><td>7.409 × 10<sup>-5</sup></td><td>1.806 × 10<sup>-4</sup></td></tr><tr><td>6</td><td>3.695 × 10<sup>-5</sup></td><td>4.944 × 10<sup>-5</sup></td><td>1.605 × 10<sup>-4</sup></td></tr><tr><td>7</td><td>1.659 × 10<sup>-5</sup></td><td>2.812 × 10<sup>-5</sup></td><td>8.099 × 10<sup>-5</sup></td></tr><tr><td>8</td><td>6.390 × 10<sup>-6</sup></td><td>1.499 × 10<sup>-5</sup></td><td>5.770 × 10<sup>-5</sup></td></tr><tr><td>9</td><td>4.341 × 10<sup>-6</sup></td><td>3.641 × 10<sup>-9</sup></td><td>2.687 × 10<sup>-5</sup></td></tr><tr><td>10</td><td>7.052 × 10<sup>-6</sup></td><td>3.202 × 10<sup>-14</sup></td><td>1.854 × 10<sup>-8</sup></td></tr><tr><td>11</td><td>8.087 × 10<sup>-10</sup></td><td></td><td>1.726 × 10<sup>-13</sup></td></tr><tr><td>12</td><td>2.739 × 10<sup>-14</sup></td><td></td><td></td></tr></tbody></table>

study case, the plasma is modeled with an effective atomic number $Z_{eff} = 3$. Moreover, the plasmas created in WEST have tungsten impurities which are
<figure><img src="image_13.png" /><figcaption>Figure 10: Case 2. Comparison of results computed at 100 [ms] with a time step $\Delta t$ (blue) and $\Delta t/2$ (red). Equilibrium plasma boundary (left) and the magnetic axis $r$ and $z$ position (right).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 11: Case 2. Comparison of results computed at 100 [ms] with a time step $\Delta t$ (blue) and $\Delta t/2$ (red). $f$ (left) and $T_e$ (right).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Case 1. Left: initial (blue) and final (red) plasma boundary. Right: $f$ profile at different time steps.</figcaption></figure>

responsible for a great part of the radiation losses. In order to simulate this effect, impurities with $Z_W = 74$ are added as explained in Appendix B. Since we don't use any model for the electron density evolution, we choose $\chi_{n_e} = 0$ leading to $Q_i = 0$. The coefficient $c_{ie}$ is computed as explained in Appendix B. Concerning the remaining parameters we choose $r_{ie} = 1.0$ for the ion-electron temperature ratio (leading to $Q_{ei} = 0$) and choose a diffusion coefficient profile $\chi_{T_e}$ proportional to the safety factor $q$ ($\chi_{T_e} = kq$ with $k = 0.15$ in all numerical experiments presented in this work).

<table><caption>Table 2: Convergence history of Newton iterations for case 2: iteration number $n$ and residual relative error $||\mathbf{X}^n - \mathbf{X}^{n-1}||/||\mathbf{X}^{n-1}||$.</caption><thead><tr><th>n</th><th>t = 0.01</th><th>t = 0.02</th><th>t = 0.03</th></tr></thead><tbody><tr><td>1</td><td>1.503 × 10<sup>-3</sup></td><td>1.230 × 10<sup>-3</sup></td><td>1.327 × 10<sup>-3</sup></td></tr><tr><td>2</td><td>2.526 × 10<sup>-4</sup></td><td>1.797 × 10<sup>-4</sup></td><td>2.359 × 10<sup>-4</sup></td></tr><tr><td>3</td><td>4.650 × 10<sup>-5</sup></td><td>3.376 × 10<sup>-5</sup></td><td>2.696 × 10<sup>-5</sup></td></tr><tr><td>4</td><td>1.283 × 10<sup>-5</sup></td><td>6.966 × 10<sup>-6</sup></td><td>9.006 × 10<sup>-6</sup></td></tr><tr><td>5</td><td>1.735 × 10<sup>-6</sup></td><td>1.608 × 10<sup>-9</sup></td><td>1.851 × 10<sup>-9</sup></td></tr><tr><td>6</td><td>5.714 × 10<sup>-11</sup></td><td>2.523 × 10<sup>-14</sup></td><td>3.029 × 10<sup>-14</sup></td></tr></tbody></table>

Changes in the internal characteristics of the plasma, non-linear interplay between the different quantities and the introduction of more realistic
4.3 Controlled evolution

A VDE is a brutal event. Its simulation can be used at first to test the robustness of the code. However during a tokamak discharge the plasma position is usually controlled by the poloidal field coils.

As a third numerical result (Case 3) we consider the use of a magnetic feedback controller developed for WEST [38] in order to prevent any VDE to occur and to drive the plasma from an initial limiter configuration to an X-point configuration. Since the code NICE is written in C++ and the controller in Matlab Simulink we have chosen to use MUSCLE3 library [49] to link one to each other. At each time step NICE communicates through MUSCLE3 the equilibrium data needed by the controller (plasma geometric axis, plasma current and coils currents) and the controller returns the computed voltages to be applied in the poloidal field circuits in order to correct the feed forward coils currents.

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

We simulate 2 seconds of plasma with a time step of $10^{-3}$ [s]. The stopping criterion for Newton iterations at each time step is set to $10^{-5}$ as it is shown in the previous section 4.2 that this does not significantly impact the result. The simulation runs the 2000 time steps in 2h30 on a laptop with an i7-13800H processor. This simulation is also faster than the previous test cases at each time step as the algorithm takes approximately 3 iterations to converge when the plasma is stabilized.

The controller takes into account the geometric axis as it has a smoother evolution than the magnetic axis if the equilibrium is computed using simple P1 FEM. This could be modified in the future since the magnetic axis evolves smoothly with rHCT FEM. The plasma has the expected behaviour. It follows the path indicated by the controller as shown on Fig. 12 and Fig. 13. The controller also tries to maintain the plasma current at 482000 [A].

<figure><img src="image_16.png" /><figcaption>Figure 12: Plasma current (Ip), coils currents (Icoils). The dashed lines show the controller requests.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Case 2. Left : initial (blue) and final (red) plasma boundary. Right : f profile at different time steps.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Comparaison between the simulation (red) and shot 59695 (blue) for the initial equilibrium (left plot), after 100 time steps (middle plot) and the evolution of the geometric axis coordinates (right plots).</figcaption></figure>

bootstrap current and resistivity profiles, seem to improve convergence of the Newton method (Tab. 2) and to keep the algorithm converging over more
<figure><img src="image_17.png" /><figcaption>Figure 13: Case 3. Plasma boundary at different time steps (top left). Then from top to bottom and left to right: magnetic (red) and geometric (blue) axis with respect to time and plasma current. The dashed lines show the controller requests. Plasma toroidal current, diamagnetic function, pressure, safety factor, electron temperature, $ff'$ and $p'$ profiles at different time steps.</figcaption></figure>

On figure 12 the effects of the voltages appear as the currents increase or decrease in the coils. The command of the currents for the coils is well respected. Coil A is left free in order to control the plasma current. Integrated power of the different sources introduced in section 2.3 over the volume of the plasma are shown on Fig. 14. We keep $Q_a$ constant at 0 [MW] for the X point formation and $Q_i$ is 0 since $\chi_{ne} = 0$. The radiative losses are about 200kW and the ohmic source about 500kW. $Q_{ei}$ is also equal to 0 as we keep $r_{ie} = 1.0$.
<figure><img src="image_9.png" /><figcaption>Figure 8: Case 2. Electron temperature, plasma pressure, resistivity and bootstrap current profiles at several time steps during the simulation</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Case 2. Source profiles at several time steps during the simulation</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 14: Power of the different sources, $Q_{ohm}$ (green), $Q_{rad}$ (red). Sources $Q_i$ (purple), $Q_a$ (blue) and $Q_{ei}$ (yellow) are superimposed with value 0.</figcaption></figure>
time steps. The simulation holds for more than 100 [ms] (Fig. 6).

As in the previous case the plasma moves upwards, switches from an X-point configuration to a limiter configuration and finally shrinks against the limiter. From Figure 7 it can be seen that the experimental plasma shape evolution from WEST shot 59695 is quite well reproduced.

Different simulated profiles at several time steps during the simulation are shown on Fig. 8 and 9. Temperature drops in the plasma center and increases near the edge as time goes on (Fig. 8a). This change in temperature leads to a greater pressure on the edge, a decrease in the resistivity and an increase in the bootstrap current on the edge. Function *f* globally increases. The Ohmic source is directly linked to the derivatives of *f*, the effect of this is a spike around $\rho_N = 0.9$.

In order to test the good numerical behavior of the code, multiple runs of the simulation are done varying the space and time discretization, but also the Newton stopping criterion or the even the replay iterations.

Varying the space and time discretization does not change significantly the solution. An example (Fig. 10 and Fig. 11) shows the differences between the simulation above and a simulation using a time step divided by 2. The results are similar if we refine the space discretization in the transport solver or if we refine the FEM mesh.

Moreover the Newton stopping criterion can be set to $10^{-5}$ instead of $10^{-10}$ without changing the result substantially. Adding the replay loop between the equilibrium solver and the transport solver does not change significantly the solution either.
# 5 Conclusion

We have presented the latest development of numerical methods implemented in the code NICE in order to simulate the evolution of the plasma equilibrium at the resistive diffusion time scale. The proposed model: the evolutive equilibrium equation for the poloidal flux, the evolution equation for the poloidal current function from [25] and an evolution equation for electron temperature, together with a coupled P1-rHCT FEM method have enabled us to obtain simulations which were conducted in quite a generic manner without particularly complicated tuning of the code in order to keep consistency between equilibrium and transport. The algorithm has proven to be robust and a VDE simulation has been validated against a WEST experiment. We have coupled the code to a magnetic feedback controller from the WEST tokamak thanks to the MUSCLE3 library in order to simulate an X-point formation scenario. This is extremeley promising and we expect to carry on this work with more complete and complex plasma discharge simulations and on other tokamaks than WEST.

# Acknowledgment

We would like to thank the MUSCLE3 team and in particular Daan van Vugt for their help in coupling NICE and the WEST magnetic feedback controller.

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 - EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.
<figure><img src="image_11.png" /><figcaption>Figure 10: Case 2. Comparison of results computed at 100 [ms] with a time step $\Delta t$ (blue) and $\Delta t/2$ (red). Equilibrium plasma boundary (left) and the magnetic axis $r$ and $z$ position (right).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Case 2. Comparison of results computed at 100 [ms] with a time step $\Delta t$ (blue) and $\Delta t/2$ (red). $f$ (left) and $T_e$ (right).</figcaption></figure>
A Test of the original numerical method from [25]

As a starting point for this paper we have implemented the same numerical methods and the
same test case as in [25]. This means that we do not use the coupled P1-rHCT FEM method
but a simpler P1 FEM method for the space discretization of (3) and (4). Moreover in (4) *f*
is decomposed in a basis of Legendre polynomials. The $T_e$ equation is not included, we take
$\langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle = 0$ and use the analytic expression $\eta^{\|}(\psi_N) = \frac{10^{-9}}{1.001-\psi_N}$. Voltages in the poloidal
field circuits are kept constant during the time stepping. Using the same time step as in [25],
$\Delta t = 10^{-1}$ [s], we obtain very similar results with a smooth evolution of *f* as shown in Fig. 15.
However if we use a more realistic time step $\Delta t = 10^{-3}$ (which is not done in [25]), numerical
oscillations appear in the approximation of *f* as shown in Fig. 16. In this case the plasma position
is not controlled and the plasma undergoes a VDE and the jumps from node to node with P1 FEM
in the X-point position defining the plasma boundary and in magnetic axis initiate an instability
in the resolution of (4). Another reason for these oscillations seem to be the quadrature method
used to compute the integrals in (4). Indeed the same quadrature rule as for (3) is used, that
is to say a single point barycentric rule for each triangle containing some plasma. This seems
not precise enough to discretize (4) with *f* represented with Legendre polynomials up to degree
10. The simulation shown in Fig. 16 stops because of the non convergence of the Newton
method after 16 iterations. Numerical oscillations on *f* and $f f'$ cause significant difficulties in
the numerical resolution.

In order to overcome these numerical issues we propose in this work to use an rHCT FEM
method to discretize $\psi$ in $\Omega_{in}$ and natural Splines of degree 3 to discretize $f$. This enables us to
represent $f$ with polynomials of much smaller degree while still using 11 (or more) basis functions.
Moreover the use of rHCT FEM has 2 advantages: first it frees the X-point and the magnetic
axis position from mesh nodes and second it necessitates the use of a higher order quadrature rule
since the basis FEM functions are polynomials of degree 3. The same quadrature rule is used to
discretize the equation for $f$.

<figure><img src="image_19.png" /><figcaption>Figure 15: Test case from [25] with $\Delta t = 10^{-1}$ [s]</figcaption></figure>
*Adiabatic evolution.* As a last test of the code in this section we conduct a similar numerical simulation but setting the resistivity to 0, as well as the sources, transport and diffusion coefficients in the $T_e$ equation. In this configuration, the plasma experiences an adiabatic evolution.

In order to perturbate the initial equilibrium, voltages in the poloidal field circuits are set to 0 for the first 7 [ms] of the simulation before being set to their values corresponding to the coils currents of the initial equilibrium. The simulation is stopped after 340 time steps (0.340 [s]) just before the plasma gets into a limiter configuration as this transition perturbates the plasma profiles.

The plasma geometry evolves during the simulation leading to an evolution of $\rho_b$. In order to check the consistency between the computed $\rho_b$ and the geometry we use the fact that $\rho_b = \sqrt{\frac{\phi_b}{\pi B_0}} \approx \sqrt{\frac{S}{\pi}} \approx a\sqrt{\kappa}$ with $S$ the surface of the plasma cross section, $a$ the minor radius and $\kappa$ the elongation. This is shown on Figure 12 where $\rho_b$ is compared to $a\sqrt{\kappa}$ throughout the simulation.

<figure><img src="image_13.png" /><figcaption>Figure 12: Left: initial (blue) and final (red) plasma boundary. Right: $a\sqrt{\kappa}$ versus $\rho_b$ (see text) at different time points of the simulation.</figcaption></figure>

During an adiabatic evolution of the plasma it is expected that quantities $p\frac{\partial V}{\partial \rho}$ and $\frac{\partial \psi}{\partial \rho}$ remain constant in time for a given $\rho$ value [7, 2]. The latter implies that the safety factor $q$ also remains constant in time for a given $\rho$ value. Figures 13 and 14 show that this is satisfied in our simulation for all $\rho$ values present at both initial and final time.

Figure 15 shows the $p, p', f$ and $ff'$ profiles at initial and final time. Even though there are no sources in the $T_e$ equation the pressure slightly evolves
<figure><img src="image_20.png" /><figcaption>Figure 16: Same test case as in [25] but with $\Delta t = 10^{-3}$ [s]</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: Initial (blue) and final (red) $p \frac{\partial V}{\partial \rho}^{\frac{5}{3}}$ profiles. The 2 dots indicate the profile value for the initial and final $\rho_b$ values.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14: Left: initial (blue) and final (red) $\frac{\partial \psi}{\partial \rho}$ profiles. Right: Initial (blue) and final (red) $q$ profiles. The 2 dots indicate the profile value for the initial and final $\rho_b$ values.</figcaption></figure>

because the plasma volume changes and the adiabatic quantity $p \frac{\partial V}{\partial \rho}^{\frac{5}{3}}$ should remain constant in time. The $f f'$ profile is perturbed on the edge of the plasma during the evolution [49]. This should be linked to the fact that the model can take into account the induced edge current observed in tokamak plasmas.
B Sources, resistivity and bootstrap current

Following [2] the electron-ion energy exchange source term $Q_{ei}$ is written as

$$Q_{ei} = \frac{3}{2} \sum_{j} \frac{n_{j} Z_{j}^{2}}{A_{j}} \frac{e(T_{e} - T_{i})}{(m_{p}/2m_{e})\tau_{e}}, \qquad (34)$$

where the energy confinement time is

$$\tau_e = \frac{12\pi^{3/2}\epsilon_0^2}{e^4} \sqrt{\frac{m_e}{2}(eT_e)^{3/2} n_e \Lambda_{ei}}, \qquad (35)$$

and the Coulomb logarithm is

$$\Lambda_{ei} = 15.2 - \frac{1}{2} \log \frac{n_e}{10^{20}} + \log \frac{T_e}{10^3}. \qquad (36)$$

The proton mass is $m_p$, the electron mass $m_e$ and the permittivity of vacuum $\epsilon_0$. For ion j, $n_j$ is its density, $Z_j$ its charge number and $A_j$ its mass number. Each ion density is proportional to $n_e$. If we consider a plasma with one main ion and one impurity, we choose $Z_{eff}$ such that $Z_{main} < Z_{eff} < Z_{imp}$ (in our case $Z_{main} = 1$ and $Z_{imp} = 7$). Then the $Z_{eff}$ equation, $n_e Z_{eff} = n_{main} Z_{main}^2 + n_{imp} Z_{imp}^2$ together with the electroneutrality equation, $n_e = n_{main} Z_{main} + n_{imp} Z_{imp}$, give $n_{main} = c_{main} n_e$, $n_{imp} = c_{imp} n_e$ and $n_i = n_{main} + n_{imp} = c_{ie} n_e$ with

$$\begin{align}
c_{main} &= \frac{Z_{imp} - Z_{eff}}{Z_{main}(Z_{imp} - Z_{main})}, \nonumber \\
c_{imp} &= \frac{Z_{eff} - Z_{main}}{Z_{imp}(Z_{imp} - Z_{main})}, \tag{37} \\
c_{ie} &= \frac{Z_{imp} + Z_{main} - Z_{eff}}{Z_{imp}Z_{main}}. \nonumber
\end{align}$$

The radiative loss term $Q_{rad}$ is computed as

$$Q_{rad} = \sum_{j} n_e n_j L_j(T_e) \qquad (38)$$

where $L_j(T_e)$ is read from a table computed from the Atomic Data and Analysis Structure ADAS database [46].

In order to add the radiative effect of the tungsten we introduce a density $n_W = c_{W}n_{e}$ of tungsten (with $c_{W} = 1.0 \times 10^{-5}$) and revise the equations as follow :

$$\begin{align}
n_e^* &= n_e - n_W Z_W, \nonumber \\
c_{main}^* &= c_{main} - c_{W_e}, \tag{39} \\
c_{ie}^* &= c_{ie} - c_{W_e}. \nonumber
\end{align}$$

Hence electronic neutrality still holds.

The resistivity and bootstrap current profiles are computed following [44]. For the resistivity we have

$$\eta^{\parallel} = \frac{1}{\sigma_{Sptz}} \qquad (40)$$
<figure><img src="image_16.png" /><figcaption>Figure 15: Initial (blue) and final (red) *p*', *ff*', *p* and *f* profiles.</figcaption></figure>
with

$$ \sigma_{Sptz} = \frac{1.9012 \cdot 10^4}{Z_{eff} N(Z_{eff}) \ln \Lambda_e} T_e^{\frac{3}{2}}, \qquad (41) $$

$$ N(Z_{eff}) = 0.58 + \frac{0.74}{0.76 + Z_{eff}}, \qquad (42) $$

and

$$ \ln \Lambda_e = 31.3 - \ln(\frac{\sqrt{n_e}}{T_e}) \qquad (43) $$

For the bootstrap current we have

$$ \langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle = fp \left[ \mathcal{L}_{31} \frac{\partial \ln n_e}{\partial \psi} + \frac{p_e}{p} (\mathcal{L}_{31} + \mathcal{L}_{32}) \frac{\partial \ln T_e}{\partial \psi} + (1 - \frac{p_e}{p}) (\mathcal{L}_{31} + \alpha \mathcal{L}_{34}) \frac{\partial \ln T_i}{\partial \psi} \right], \quad (44) $$

and the expressions for profiles $\mathcal{L}_{31}$, $\mathcal{L}_{32}$, $\mathcal{L}_{34}$ and $\alpha$ are detailed in [44] and highly depend on trapped fraction profile.

## C Derivatives computation

The derivatives of the equilibrium system of equations are computed analytically and read as follows. Their accuracy is verified by the gradient test as shown on Fig. 4.

Derivative of the equation for $\psi$ with respect to $\psi$:

$$ (D_{X_{\psi}} \mathbf{e}(\mathbf{X}_{\psi}^{n}, \mathbf{X}_{f}^{n}))_{i,j} = \int_{\Omega^{\text{in}}} A + C - \frac{\partial J_{i}}{\partial \psi_{j}} drdz \qquad (45) $$

with

$$ \begin{aligned} \frac{\partial J_i}{\partial \psi_j} = & r[p'']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} w_i + r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} w_i \\ & + \frac{1}{\mu_0 r} \left( \left( [f']_k^n(\psi_{N,h}^n) \right)^2 + [f'']_k^n(\psi_{N,h}^n) f_k^n(\psi_{N,h}^n) \right) \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} w_i \\ & + \frac{1}{\mu_0 r} [f']_k^n(\psi_{N,h}^n) f_k^n(\psi_{N,h}^n) \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} w_i \end{aligned} \qquad (46) $$

Derivative of the equation for $f$ with respect to $\psi$:

$$ (D_{X_{\psi}} \mathbf{d}(\mathbf{X}_{\psi}^{n}, \mathbf{X}_{f}^{n}))_{i,j} = \int_{\Omega^{\text{in}}} \frac{\partial R D_i}{\partial \psi_j} dr dz \qquad (47) $$

with
### 4.3. Controlled evolution

Test case 3. Free plasma displacement simulations can be used at first to test the robustness of the code. However during a tokamak discharge the plasma position is controlled by the poloidal field coils.

Hence we consider now the use of a magnetic feedback controller developed for WEST [27] in order to prevent any VDE to occur and to drive the plasma from an initial limiter configuration to an X-point configuration. Since the code NICE is written in C++ and the controller in Matlab Simulink we have chosen to use MUSCLE3 library [50] to link one to each other. At each time step NICE communicates through MUSCLE3 the equilibrium data needed by the controller (plasma geometric axis, plasma current and coils currents) and the controller returns the computed voltages to be applied in the poloidal field circuits in order to correct the feed forward coils currents.

<table><tr><td>Muscle3</td></tr></table>

↑ *dataEqui* ↓ *Voltages*

↑ *Voltages* ↓ *dataEqui*

NICE

Controller

We simulate 2 [s] of plasma with a time step of 10⁻³ [s]. The stopping criterion for Newton iterations at each time step is set to 10⁻⁵ as it is shown in the previous Section 4.2 that this does not significantly impact the result.

The simulation runs the 2000 time steps in 2h30 on a laptop with an i7-13800H processor. This simulation is also faster than the previous test cases at each time step as the algorithm takes approximately 3 iterations to converge when the plasma is stabilized.

The controller takes into account the geometric axis as it has a smoother evolution than the magnetic axis if the equilibrium is computed using simple P1 FEM. This could be modified in the future since the magnetic axis evolves smoothly with rHCT FEM. The plasma has the expected behaviour. It follows the path indicated by the controller as shown on Fig. 16 and Fig. 17. The controller also tries to maintain the plasma current at 482000 [A].

On Fig. 16 the effects of the voltages appear as the currents increase or decrease in the coils. The command of the currents for the coils is well respected. Coil A is left free in order to control the plasma current. Integrated power of the different sources introduced in Section 2.3 over the volume of the plasma are shown on Fig. 18. We keep $Q_a$ constant at 0 [MW] for the X point formation and $Q_i$ is 0 since $\chi_{ne} = 0$. The radiative losses are about
$$
\begin{equation}
\begin{aligned}
\frac{\partial RD_i}{\partial \psi_j} = & \frac{[f']_k^n(\psi_{N,h}^n)}{\Delta t} \frac{\partial \psi_N}{\partial \psi_j} \frac{1}{r} g_i(\psi_{N,h}^n) + \frac{f_k^n(\psi_{N,h}^n) - f_k^{n-1}(\psi_{N,h}^{n-1})}{\Delta t} \frac{1}{r} g'_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \\
& + \frac{1}{\Delta t} \frac{\partial \psi_h^n}{\partial \psi_j} \frac{f_k^n(\psi_{N,h}^n)}{r} g'_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} + \frac{\psi_h^n - \psi_h^{n-1}}{\Delta t} \frac{[f']_k^n(\psi_{N,h}^n)}{r} \frac{\partial \psi_N}{\partial \psi_j} g'_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \\
& + \frac{\psi_h^n - \psi_h^{n-1}}{\Delta t} \frac{f_k^n(\psi_{N,h}^n)}{r} g''_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} + \frac{\psi_h^n - \psi_h^{n-1}}{\Delta t} \frac{f_k^n(\psi_{N,h}^n)}{r} g'_i(\psi_{N,h}^n) \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} \\
& + [\eta']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \left[ r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} f_k^n(\psi_{N,h}^n) \right. \\
& \left. + \frac{[f_k^n]^2(\psi_{N,h}^n) + |\nabla \psi_h^n|^2}{\mu_0 r} [f']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} - r[jB]_k^{n-1}(\psi_{N,h}^n) \right] g'_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \\
& + \eta_k^{n-1}(\psi_{N,h}^n) \left[ r[p'']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} f_k^n(\psi_{N,h}^n) + r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} [f']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \right. \\
& \left. + r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} f_k^n(\psi_{N,h}^n) + 2 \frac{f_k^n(\psi_{N,h}^n)}{\mu_0 r} ([f']_k^n(\psi_{N,h}^n))^2 \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} \right. \\
& \left. + \frac{2\nabla \psi}{\mu_0 r} \cdot \nabla \left( \frac{\partial \psi_h}{\partial \psi_j} \right) [f']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} + \frac{[f_k^n]^2(\psi_{N,h}^n) + |\nabla \psi_h^n|^2}{\mu_0 r} [f''']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} \right. \\
& \left. + \frac{[f_k^n]^2(\psi_{N,h}^n) + |\nabla \psi_h^n|^2}{\mu_0 r} [f']_k^n(\psi_{N,h}^n) \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} - r[jB']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \right] g'_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \\
& + \eta_k^{n-1}(\psi_{N,h}^n) \left[ r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} f_k^n(\psi_{N,h}^n) \right. \\
& \left. + \frac{[f_k^n]^2(\psi_{N,h}^n) + |\nabla \psi_h^n|^2}{\mu_0 r} [f']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} - r[jB]_k^{n-1}(\psi_{N,h}^n) \right] g''_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} \\
& + \eta_k^{n-1}(\psi_{N,h}^n) \left[ r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} f_k^n(\psi_{N,h}^n) \right. \\
& \left. + \frac{[f_k^n]^2(\psi_{N,h}^n) + |\nabla \psi_h^n|^2}{\mu_0 r} [f']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} - r[jB]_k^{n-1}(\psi_{N,h}^n) \right] g'_i(\psi_{N,h}^n) \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j}
\end{aligned}
\tag{48}
\end{equation}
$$

Derivative of the equation for $\psi$ with respect to $f$:

$$ (D_{X_f} \mathbf{e}(\mathbf{X}_{\psi}^n, \mathbf{X}_f^n))_{i,j} = \int_{\Omega^{\text{in}}} \frac{\partial J_i}{\partial f_j} dr dz \quad (49) $$

with

$$ \frac{\partial J_i}{\partial f_j} = \frac{1}{\mu_0 r} [g'_i(\psi_{N,h}^n) f_k^n(\psi_{N,h}^n) + [f']_k^n(\psi_{N,h}^n) g_i(\psi_{N,h}^n)] \frac{\partial \psi_N}{\partial \psi} w_i \quad (50) $$
<figure><img src="image_17.png" /><figcaption>Figure 16: Plasma current (Ip), coils currents (Icoils). The dashed lines show the controller requests.</figcaption></figure>

200kW and the ohmic source about 500kW. The quantity $Q_{ei}$ is also equal to 0 as we keep $r_{ie} = 1.0$.

620
Derivative of the equation for f with respect to f:

$$ (D_{X_f} \mathbf{d}(\mathbf{X}_{\psi}^n, \mathbf{X}_f^n))_{i,j} = \int_{\Omega^{\text{in}}} \frac{\partial RD_i}{\partial f_j} dr dz \quad (51) $$

with

$$ \frac{\partial RD_i}{\partial f_j} = \frac{1}{\Delta t} \frac{g_j(\psi_{N,h}^n)}{r} g_i(\psi_{N,h}^n) + \frac{\psi_h^n - \psi_h^{n-1}}{\Delta t} \frac{g_j(\psi_{N,h}^n)}{r} g_i'(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \\ + \eta_k^{n-1}(\psi_{N,h}^n) \left[ r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} g_j(\psi_{N,h}^n) + \frac{2f_k^n(\psi_{N,h}^n)g_j(\psi_{N,h}^n)}{\mu_0 r} [f']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \right. \\ \left. + \frac{[f_k^n]^2(\psi_{N,h}^n) + |\nabla \psi_h^n|^2}{\mu_0 r} g_j'(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \right] g_i'(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \quad (52) $$

and

$$ \frac{\partial \psi_N}{\partial \psi_j} = \frac{\partial \psi_N}{\partial \psi} \frac{\partial \psi_h}{\partial \psi_j} + \frac{\partial \psi_N}{\partial \psi_a} \frac{\partial \psi_a(\psi_h)}{\partial \psi_j} + \frac{\partial \psi_N}{\partial \psi_b} \frac{\partial \psi_b(\psi_h)}{\partial \psi_j} \\ \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} = \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_a} \frac{\partial \psi_a(\psi_h)}{\partial \psi_j} + \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_b} \frac{\partial \psi_b(\psi_h)}{\partial \psi_j} $$

with

$$ \frac{\partial \psi_N}{\partial \psi} = \frac{1}{\psi_b - \psi_a} \\ \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_a} = \frac{1}{(\psi_b - \psi_a)^2} \\ \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_b} = \frac{-1}{(\psi_b - \psi_a)^2} $$
<figure><img src="image_18.png" /><figcaption>Figure 17: Case 3. Plasma boundary at different time steps (top left). Then from top to bottom and left to right: magnetic (red) and geometric (blue) axis with respect to time and plasma current. The dashed lines show the controller requests. Plasma toroidal current, diamagnetic function, pressure, safety factor, electron temperature, $ff'$ and $p'$ profiles at different time steps.</figcaption></figure>
References

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

[2] J.F. Artaud, V. Basiuk, F. Imbeaux, M. Schneider, J. Garcia, G. Giruzzi, P. Huynh, T. Aniel, F. Albajar, J.M. Ané, A. Bécoulet, C. Bourdelle, A. Casati, L. Colas, J. Decker, R. Dumont, L.G. Eriksson, X. Garbet, R. Guirlet, P. Hertout, G.T. Hoang, W. Houlberg, G. Huysmans, E. Joffrin, S.H. Kim, F. Köchl, J. Lister, X. Litaudon, P. Maget, R. Masset, B. Pégourié, Y. Peysson, P. Thomas, E. Tsitrone, and F. Turco. The CRONOS suite of codes for integrated tokamak modelling. *Nuclear Fusion*, 50(4):043001, 2010.

[3] C. Bernardi, Y. Maday, and A.T. Patera. 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.

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

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

[6] J. Blum and J. Le Foll. Plasma equilibrium evolution at the resistive diffusion timescale. *Computer Physics Reports*, 1(7-8):465–494, 1984.

[7] J. Blum, J. Le Foll, and B. Thooris. The self-consistent equilibrium and diffusion code SCED. *Computer Physics Communications*, 24(3):235–254, 1981.

[8] C. Boulbe, B. Faugeras, G. Gros, and F. Rapetti. Tokamak free-boundary plasma equilibrium computations in presence of non-linear materials. *J. Scientific Computing*, 96:42, 2023. https://hal.sciences/hal-03423469.

[9] C. Bourdelle, J.F. Artaud, V. Basiuk, M. Bécoulet, S. Brémond, J. Bucalossi, H. Bufferand, G. Ciraolo, L. Colas, Y. Corre, X. Courtois, J. Decker, L. Delpech, P. Devynck, G. Dif-Pradalier, R.P. Doerner, D. Douai, R. Dumont, A. Ekedahl, N. Fedorczak, C. Fenzi, M. Fir-daouss, J. Garcia, P. Ghendrih, C. Gil, G. Giruzzi, M. Goniche, C. Grisolia, A. Gros-man, D. Guilhem, R. Guirlet, J. Gunn, P. Hennequin, J. Hillairet, T. Hoang, F. Imbeaux, I. Ivanova-Stanik, E. Joffrin, A. Kallenbach, J. Linke, T. Loarer, P. Lotte, P. Maget, Y. Marandet, M.L. Mayoral, O. Meyer, M. Missirlian, P. Mollard, P. Monier-Garbet, P. Moreau, E. Nardon, B. Pégourié, Y. Peysson, R. Sabot, F. Saint-Laurent, M. Schnei-der, J.M. Travère, E. Tsitrone, S. Vartanian, L. Vermare, M. Yoshida, R. Zagorski, and JET Contributors. WEST Physics Basis. *Nuclear Fusion*, 55(6):063017, June 2015.

[10] J. Bucalossi, J. Achard, O. Agullo, T. Alarcon, L. Allegretti, H. Ancher, G. Antar, S. An-tusch, V. Anzallo, C. Arnas, D. Arranger, J. F. Artaud, M. H. Aumeunier, S. G. Baek, X. Bai, J. Balbin, C. Balorin, T. Barbui, A. Barbuti, J. Barlerin, V. Basiuk, T. Batal, O. Baulaigue, A. Bec, M. Bécoulet, E. Benoit, E. Benard, J. M. Benard, N. Bertelli, E. Bertrand, P. Beyer, J. Bielecki, P. Bienvenu, R. Bisson, V. Bobkov, G. Bodner, C. Bot-tereau, C. Bouchand, F. Bouquey, C. Bourdelle, J. Bourg, S. Brezinsek, F. Brochard, C. Brun, V. Bruno, H. Bufferand, A. Bureau, S. Burles, Y. Camenen, and Antti Hakola. Op-erating a full tungsten actively cooled tokamak: Overview of WEST first phase of operation. *Nuclear Fusion*, 62(4), February 2022.
<figure><img src="image_19.png" /><figcaption>Figure 18: Power of the different sources, $Q_{ohm}$ (green), $Q_{rad}$ (red). Sources $Q_i$ (purple), $Q_a$ (blue) and $Q_{ei}$ (yellow) are superimposed with value 0.</figcaption></figure>
[11] J. Bucalossi, M. Missirlian, P. Moreau, F. Samaille, E. Tsitrone, D. Van Houtte, T. Batal, C. Bourdelle, M. Chantant, Y. Corre, X. Courtois, L. Delpech, L. Doceul, D. Douai, H. Dougnac, F. Faisse, C. Fenzi, F. Ferlay, M. Firdaouss, L. Gargiulo, P. Garin, C. Gil, A. Grosman, D. Guilhem, J. Gunn, C. Hernandez, D. Keller, S. Larroque, F. Leroux, M. Lipa, P. Lotte, A. Martinez, O. Meyer, F. Micolon, P. Mollard, E. Nardon, R. Nouailletas, A. Pilia, M. Richou, S. Salasca, and J.-M. Travère. The WEST project: Testing ITER divertor high heat flux component technology in a steady state tokamak environment. *Fusion Engineering and Design*, 89(7-8):907-912, October 2014.

[12] R.N Byrne and H.H Klein. G2M, a two-dimensional diffusion time scale tokamak code. *Journal of Computational Physics*, 26(3):352-378, 1978.

[13] R.W. Clough and J.L. Tocher. 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, October 1965.

[14] D.P. Coster, V. Basiuk, G. Pereverzev, D. Kalupin, R. Zagorksi, R. Stankiewicz, P. Huynh, and F. Imbeaux. The european transport solver. *Plasma Science, IEEE Transactions on*, 38(9):2085-2092, Sep. 2010.

[15] J.A. Crotinger. CORSICA; a comprehensive simulation of toroidal magnetic-fusion devices. Technical report uclr-id-126284, Lawrence Livermore National Laboratory, 1997.

[16] C. De Boor. *A Practical Guide To Splines*. Springer-Verlag, 1978.

[17] A. Elarif, B. Faugeras, and F. Rapetti. Tokamak free-boundary plasma equilibrium computation using finite elements of class C0 and C1 within a mortar element approach. *J. Computational Physics*, 439:110388, 2021.

[18] E. Fable, C. Angioni, F. J. Casson, D. Told, A. A. Ivanov, F. Jenko, R. M. McDermott, S. Yu. Medvedev, G. V. Pereverzev, F. Ryter, W. Treutterer, E. Viezzer, and the ASDEX Upgrade Team. Novel free-boundary equilibrium and transport solver with theory-based models and its validation against ASDEX Upgrade current ramp scenarios. *Plasma Physics and Controlled Fusion*, 55(12):124028, nov 2013.

[19] G. L. Falchetto, D. Coster, R. Coelho, B.D. Scott, L. Figini, D. Kalupin, E. Nardon, L.L. Alves, J.F. Artaud, V. Basiuk, J. Bizarro, C. Boulbe, A. Dinklage, D. Farina, B. Faugeras, J. Ferreira, A. Figueiredo, P. Huynh, F. Imbeaux, I. Ivanova-Stanik, T. Jonsson, H.-J. Klingshirn, C. Konz, A. Kus, N.B. Marushchenko, E. Nardon, S. Nowak, G. Pereverzev, M. Owsiak, E. Poli, Y. Peysson, R. Reimer, J. Signoret, O. Sauter, R. Stankiewicz, P. Strand, I. Voitsekhovitch, E. Westerhof, T. Zok, W. Zwingmann, ITM-TF contributors, the ASDEX Upgrade Team, and JET-EFDA Contributors. The European Integrated Tokamak Modelling (ITM) effort: achievements and first physics results. *Nuclear Fusion*, 54(4):043018, 2014.

[20] B. Faugeras. An overview of the numerical methods for tokamak plasma equilibrium computation implemented in the NICE code. *Fusion Eng. Design*, 160:112020, 2020. https://hal.archives-ouvertes.fr/hal-02955053.

[21] F. Felici, O. Sauter, S. Coda, B.P. Duval, T.P. Goodman, J-M. Moret, J.I. Paley, and the TCV Team. Real-time physics-model-based simulation of the current density profile in tokamak plasmas. *Nuclear Fusion*, 51(8):083052, aug 2011.

[22] H. Grad and J. Hogan. Classical diffusion in a tokamak. *Phys. Rev. Lett.*, 24:1337-1340, Jun 1970.
*Vertical Displacement Event.* Another test to verify the good behaviour of the plasma with our proposed approach is a VDE simulation. This instability appears in tokamaks when the plasma position is not controlled vertically. In such a situation the plasma moves up or down exponentially fast. In our simulation we keep a controller on the radial PF coils (Fh and Fb in WEST) while the other voltages are set to 0 in order to maintain the radial position of the plasma so that it does not hit the low field side limiter. Thus after a plasma stabilization period of 0.7 [s] with the controller fully on, almost all the generators are shut down, only the radial control is kept active. The same parameters as for the previous simulations are used. Figure 19 shows the evolution of the plasma boundary and the magnetic axis, as well as the exponential decrease of the z coordinate of the magnetic axis after $t = 0.7$ [s].

<figure><img src="image_20.png" /><figcaption>Figure 19: Left: initial (blue) and final (red) plasma shape. Right: evolution of the z coordinate of the plasma magnetic axis. The vertical control is turned off at $t = 0.7$ [s]</figcaption></figure>
[23] H. Grad, P.N. Hu, and D.C. Stevens. Adiabatic evolution of plasma equilibrium. Proc. Nat. Acad. Sci. USA, 72:3789-3793, Oct. 1975.

[24] H. Grad and H. Rubin. Hydromagnetic equilibria and force-free fields. Proceedings of the 2nd UN Conf. on the Peaceful Uses of Atomic Energy, 31:190, 1958.

[25] H. Heumann. A Galerkin method for the weak formulation of current diffusion and force balance in tokamak plasmas. Journal of Computational Physics, 442:110483, 2021.

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

[27] H. Heumann and F. Rapetti. A finite element method with overlapping meshes for free-boundary axisymmetric plasma equilibria in realistic geometries. J. Comput. Phys., 334:522-540, 2017.

[28] F.L. Hinton and R.D. Hazeltine. Theory of plasma transport in toroidal confinement systems. Rev. Mod. Phys., 48:239-308, Apr 1976.

[29] M. Honda. Simulation technique of free-boundary equilibrium evolution in plasma ramp-up phase. Computer Physics Communications, 181(9):1490-1500, SEP 2010.

[30] S. C. Jardin, N. Pomphrey, and J. DeLucia. Dynamic modeling of transport and positional control of tokamaks. J. Comput. Phys., 66:481-507, October 1986.

[31] S.C. Jardin. Computational methods in plasma physics. Boca Raton, FL: CRC Press/Taylor & Francis, 2010.

[32] D. Kalupin, I. Ivanova-Stanik, I. Voitsekhovitch, J. Ferreira, D. Coster, L.L. Alves, Th. Aniel, J.F. Artaud, V. Basiuk, Joao P.S. Bizarro, R. Coelho, A. Czarnecka, Ph. Huynh, A. Figueiredo, J. Garcia, L. Garzotti, F. Imbeaux, F. Kochl, M.F. Nave, G. Pereverzev, O. Sauter, B.D. Scott, R. Stankiewicz, P. Strand, ITM-TF contributors, and JET-EFDA Contributors. Numerical analysis of JET discharges with the European Transport Simulator. Nuclear Fusion, 53(12):123007, nov 2013.

[33] D. Kalupin, G. Pereverzev, and R. Stankiewicz. European Transport Simulator: transport equations, variables and Fortran implementation. Technical documentation, Eurofusion, (revision 01.10.2015) 2015.

[34] S. Kerboua-Benlarbi, R. Nouailletas, B. Faugeras, E. Nardon, and P. Moreau. Magnetic control of WEST plasmas through deep reinforcement learning. In SOFE 2023, 30th IEEE Symposium On Fusion Engineering, Oxford, UK, July 2023.

[35] R.R. Khayrutdinov and V.E. Lukash. Studies of plasma equilibrium and transport in a tokamak fusion device with the inverse-variable technique. Journal of Computational Physics, 109(2):193-201, 1993.

[36] S. H. Kim, J. F. Artaud, V. Basiuk, V. Dokuka, R. R. Khayrutdinov, J. B. Lister, and V. E. Lukash. Full tokamak discharge simulation of ITER by combining DINA-CH and CRONOS. Plasma Physics and Controlled Fusion, 51(10):105007, 2009.

[37] J.N. Lyness and D. Jespersen. Moderate degree symmetric quadrature rules for the triangle. J. Inst. Maths. Applics., 15:19-32, 1975.
## 5. Conclusion

We have presented the latest development of numerical methods implemented in the code NICE in order to simulate the evolution of the plasma equilibrium at the resistive diffusion time scale. The proposed model: the evolutive equilibrium equation for the poloidal flux, the evolution equation for the poloidal current function from [19] and an evolution equation for electron temperature, together with a coupled P1-rHCT FEM method have enabled us to obtain simulations which were conducted in quite a generic manner without particularly complicated tuning of the code in order to keep consistency between equilibrium and transport. The algorithm has proven to be robust and a free plasma displacement simulation has been validated against a WEST experiment. We have coupled the code to a magnetic feedback controller from the WEST tokamak thanks to the MUSCLE3 library in order to simulate an X-point formation scenario. The goal of the present work is to present and test the robustness of the proposed new numerical approach for the simulation of plasma equilibrium evolution. This is extremely promising and we expect to carry on this work with more complete and complex plasma discharge simulations and on other tokamaks than WEST. In order to produce real physics results the interfacing of NICE with MUSCLE3 will be of great help to couple the code with a code for the evolution of the electron density and a code computing more realistic transport coefficients.

We would like to thank the MUSCLE3 team and in particular Daan van Vugt for their help in coupling NICE and the WEST magnetic feedback controller. We would also like to thank the two anonymous reviewers for their constructive remarks from which the paper significantly benefited.

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 - EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.
[38] R. Nouailletas, E. Nardon, P. Moreau, C. Reux, and T.-P.-H. Truong. WEST magnetic control. In *2019 IEEE 58th Conference on Decision and Control (CDC)*, pages 3214–3219, 2019.

[39] V. Ostuni, J.F. Artaud, G. Giruzzi, E. Joffrin, H. Heumann, and H. Urano. Tokamak discharge simulation coupling free-boundary equilibrium and plasma model with application to JT-60SA. *Nuclear Fusion*, 61(2):026021, jan 2021.

[40] V. Parail, R. Albanese, R. Ambrosino, J.-F. Artaud, K. Besseghir, M. Cavinato, G. Corrigan, J. Garcia, L. Garzotti, Y. Gribov, F. Imbeaux, F. Koechl, C.V. Labate, J. Lister, X. Litaudon, A. Loarte, P. Maget, M. Mattei, D. McDonald, E. Nardon, G. Saibene, R. Sartori, and J. Urban. Self-consistent simulation of plasma scenarios for ITER using a combination of 1.5d transport codes and free-boundary equilibrium codes. *Nuclear Fusion*, 53(11):113002, sep 2013.

[41] G.V. Pereverzev and P.N. Yushmanov. ASTRA Automated System for TRansport Analysis. Ipp report 5/98, Max-Planck-Institut fur Plasmaphysik, Feb. 2002.

[42] O. Sauter. Geometric formulas for system codes including the effect of negative triangularity. *Fusion Engineering and Design*, 112:633-645, 2016.

[43] O. Sauter, C. Angioni, and Y. R. Lin-Liu. Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime. *Physics of Plasmas*, 6(7):2834-2839, 07 1999.

[44] O. Sauter, C. Angioni, and Y. R. Lin-Liu. Erratum: "neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime" [phys. plasmas 6, 2834 (1999)]. *Physics of Plasmas*, 9(12):5140-5140, 11 2002.

[45] V.D. Shafranov. On magnetohydrodynamical equilibrium configurations. *Soviet Journal of Experimental and Theoretical Physics*, 6:545, 1958.

[46] H.P. Summers. The ADAS user manual, version 2.6. Technical report, http://www.adas.ac.uk, 2004.

[47] J.B. Taylor. Evolution of resistive plasma in an axisymmetric torus. Ppn-10/75, Theory division, Culham Laboratory, 1975.

[48] A.D. Turnbull and R.G. Storer. A plasma resistive diffusion model. *Journal of Computational Physics*, 50(3):409-435, 1983.

[49] L. E. Veen and A. G. Hoekstra. Easing Multiscale Model Design and Coupling with MUSCLE 3. In Valeria V. Krzhizhanovskaya, Gábor Závodszky, Michael H. Lees, Jack J. Dongarra, Peter M. A. Sloot, Sérgio Brissos, and João Teixeira, editors, *Computational Science - ICCS 2020*, pages 425-438, Cham, 2020. Springer International Publishing.

[50] E. Witrant, E. Joffrin, S. Brémond, G. Giruzzi, D. Mazon, O. Barana, and P. Moreau. A control-oriented model of the current profile in tokamak plasma. *Plasma Physics and Controlled Fusion*, 49(7):1075, 2007.
Appendix A. Test of the original numerical method from [19]

As a starting point for this paper we have implemented the same numerical methods and the same test case as in [19]. This means that we do not use the coupled P1-rHCT FEM method but a simpler P1 FEM method for the space discretization of (3) and (4). Moreover in (4) *f* is decomposed in a basis of Legendre polynomials. The *T*<sub>e</sub> equation is not included, we take <**j**<sub>ni</sub>·**B**> = 0 and use the analytic expression η<sup>||</sup>(ψ<sub>N</sub>) = <10<sup>-9</sup> / (1.001-ψ<sub>N</sub>).* Voltages in the poloidal field circuits are kept constant during the time stepping. Using the same time step as in [19], Δ*t* = 10<sup>-1</sup> [s], we obtain very similar results with a smooth evolution of *f* as shown in Fig. A.20. However if we use a more realistic time step Δ*t* = 10<sup>-3</sup> (which is not done in [19]), numerical oscillations appear in the approximation of *f* as shown in Fig. A.21. In this case the plasma position is not controlled and the plasma undergoes a VDE and the jumps from node to node with P1 FEM in the X-point position defining the plasma boundary and in magnetic axis initiate an instability in the resolution of (4). Another reason for these oscillations seem to be the quadrature method used to compute the integrals in (4). Indeed the same quadrature rule as for (3) is used, that is to say a single point barycentric rule for each triangle containing some plasma. This seems not precise enough to discretize (4) with *f* represented with Legendre polynomials up to degree 10. The simulation shown in Fig. A.21 stops because of the non convergence of the Newton method after 16 iterations. Numerical oscillations on *f* and *ff'* cause significant difficulties in the numerical resolution.

In order to overcome these numerical issues we propose in this work to use an rHCT FEM method to discretize $\psi$ in $\Omega_{in}$ and natural Splines of degree 3 to discretize $f$. This enables us to represent $f$ with polynomials of much smaller degree while still using 11 (or more) basis functions. Moreover the use of rHCT FEM has 2 advantages: first it frees the X-point and the magnetic axis position from mesh nodes and second it necessitates the use of a higher order quadrature rule since the basis FEM functions are polynomials of degree 3. The same quadrature rule is used to discretize the equation for $f$.
informatics mathematics
Inria

RESEARCH CENTRE
SOPHIA ANTIPOLIS - MÉDITERRANÉE

2004 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
<figure><img src="image_21.png" /><figcaption>Figure A.20: Test case from [19] with $\Delta t = 10^{-1}$ [s]</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure A.21: Same test case as in [19] but with $\Delta t = 10^{-3}$ [s]</figcaption></figure>
Appendix B. Sources, resistivity and bootstrap current

Following [37] the electron-ion energy exchange source term $Q_{ei}$ is written as

$$Q_{ei} = \frac{3}{2} \sum_{j} \frac{n_{j} Z_{j}^{2}}{A_{j}} \frac{e(T_{e} - T_{i})}{(m_{p}/2m_{e})\tau_{e}}, \qquad (B.1)$$

where the energy confinement time is

$$\tau_e = \frac{12\pi^{\frac{3}{2}}\epsilon_0^2}{e^4} \sqrt{\frac{m_e}{2} \frac{(eT_e)^{\frac{3}{2}}}{n_e \Lambda_{ei}}}, \qquad (B.2)$$

and the Coulomb logarithm is

$$\Lambda_{ei} = 15.2 - \frac{1}{2} \log \frac{n_e}{10^{20}} + \log \frac{T_e}{10^3}. \qquad (B.3)$$

The proton mass is $m_p$, the electron mass $m_e$ and the permittivity of vacuum $\epsilon_0$. For ion $j$, $n_j$ is its density, $Z_j$ its charge number and $A_j$ its mass number.

Each ion density is proportional to $n_e$. If we consider a plasma with one main ion and one impurity, we choose $Z_{eff}$ such that $Z_{main} < Z_{eff} < Z_{imp}$

(in our case $Z_{main} = 1$ and $Z_{imp} = 7$). Then the $Z_{eff}$ equation, $n_e Z_{eff} = n_{main} Z_{main}^2 + n_{imp} Z_{imp}^2$ together with the electroneutrality equation, $n_e = n_{main} Z_{main} + n_{imp} Z_{imp}$, give $n_{main} = c_{main} n_e$, $n_{imp} = c_{imp} n_e$ and $n_i = n_{main} + n_{imp} = c_{ie} n_e$ with

$$\begin{align}
c_{main} &= \frac{Z_{imp} - Z_{eff}}{Z_{main}(Z_{imp} - Z_{main})}, \nonumber \\
c_{imp} &= \frac{Z_{eff} - Z_{main}}{Z_{imp}(Z_{imp} - Z_{main})}, \tag{B.4} \\
c_{ie} &= \frac{Z_{imp} + Z_{main} - Z_{eff}}{Z_{imp}Z_{main}}. \nonumber
\end{align}$$

The radiative loss term $Q_{rad}$ is computed as

$$Q_{rad} = \sum_{j} n_{e} n_{j} L_{j}(T_{e}) \qquad (B.5)$$

where $L_j(T_e)$ is read from a table computed from the Atomic Data and Analysis Structure ADAS database [51, 52, 53].
In order to add the radiative effect of the tungsten we introduce a density 716 $n_W = c_{We}n_e$ of tungsten (with $c_{We} = 1.0 \times 10^{-5}$) and revise the equations as follow :

$$ 
\begin{aligned}
n_e^* &= n_e - n_W Z_W, \\
c_{main}^* &= c_{main} - c_{We}, \\
c_{ie}^* &= c_{ie} - c_{We}.
\end{aligned} 
\qquad (B.6) $$

Hence electronic neutrality still holds. 718

The resistivity and bootstrap current profiles are computed following [43].

For the resistivity we have 720

$$ \eta^{\parallel} = \frac{1}{\sigma_{Sptz}} \qquad (B.7) $$

with

$$ \sigma_{Sptz} = \frac{1.9012 \cdot 10^4}{Z_{eff} N(Z_{eff}) \ln \Lambda_e} T_e^{\frac{3}{2}}, \qquad (B.8) $$

$$ N(Z_{eff}) = 0.58 + \frac{0.74}{0.76 + Z_{eff}}, \qquad (B.9) $$

and 722

$$ \ln \Lambda_e = 31.3 - \ln(\frac{\sqrt{n_e}}{T_e}) \qquad (B.10) $$

For the bootstrap current we have

$$ \langle \mathbf{j}_{ni} \cdot \mathbf{B} \rangle = fp \left[ \mathcal{L}_{31} \frac{\partial \ln n_e}{\partial \psi} + \frac{p_e}{p} (\mathcal{L}_{31} + \mathcal{L}_{32}) \frac{\partial \ln T_e}{\partial \psi} + (1 - \frac{p_e}{p}) (\mathcal{L}_{31} + \alpha \mathcal{L}_{34}) \frac{\partial \ln T_i}{\partial \psi} \right], \qquad (B.11) $$

and the expressions for profiles $\mathcal{L}_{31}$, $\mathcal{L}_{32}$, $\mathcal{L}_{34}$ and $\alpha$ are detailed in [43] and highly depend on trapped fraction profile. 724

## 726 Appendix C. Derivatives computation

The derivatives of the equilibrium system of equations are computed analytically and read as follows. Their accuracy is verified by the gradient test as shown on Fig. 4. 728

Derivative of the equation for $\psi$ with respect to $\psi$: 730

$$ (D_{X_{\psi}} \mathbf{e}(\mathbf{X}_{\psi}^n, \mathbf{X}_f^n))_{i,j} = \int_{\Omega^{\text{in}}} A + C - \frac{\partial J_i}{\partial \psi_j} drdz \qquad (C.1) $$
with

$$ \begin{aligned} \frac{\partial J_i}{\partial \psi_j} ={}& r[p^{\prime\prime}]_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} w_i + r[p^{\prime}]_k^{n-1}(\psi_{N,h}^n) \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} w_i \\ & + \frac{1}{\mu_0 r} \left( ([f^{\prime}]_k^n(\psi_{N,h}^n))^2 + [f^{\prime\prime}]_k^n(\psi_{N,h}^n) f_k^n(\psi_{N,h}^n) \right) \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} w_i \quad (C.2) \\ & + \frac{1}{\mu_0 r} [f^{\prime}]_k^n(\psi_{N,h}^n) f_k^n(\psi_{N,h}^n) \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} w_i \end{aligned} $$

732

Derivative of the equation for $f$ with respect to $\psi$:

$$ (D_{X_{\psi}} \mathbf{d}(\mathbf{X}_{\psi}^n, \mathbf{X}_f^n))_{i,j} = \int_{\Omega_{\text{in}}} \frac{\partial R D_i}{\partial \psi_j} dr dz \quad (C.3) $$

with
$$
\begin{align}
\frac{\partial RD_i}{\partial \psi_j} = & \frac{[f']_k^n(\psi_{N,h}^n)}{\Delta t} \frac{\partial \psi_N}{\partial \psi_j} \frac{1}{r} g_i(\psi_{N,h}^n) + \frac{f_k^n(\psi_{N,h}^n) - f_k^{n-1}(\psi_{N,h}^{n-1})}{\Delta t} \frac{1}{r} g'_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \nonumber \\
& + \frac{1}{\Delta t} \frac{\partial \psi_h^n}{\partial \psi_j} \frac{f_k^n(\psi_{N,h}^n)}{r} g'_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} + \frac{\psi_h^n - \psi_h^{n-1}}{\Delta t} \frac{[f']_k^n(\psi_{N,h}^n)}{r} \frac{\partial \psi_N}{\partial \psi_j} g'_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \nonumber \\
& + \frac{\psi_h^n - \psi_h^{n-1}}{\Delta t} \frac{f_k^n(\psi_{N,h}^n)}{r} g''_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} + \frac{\psi_h^n - \psi_h^{n-1}}{\Delta t} \frac{f_k^n(\psi_{N,h}^n)}{r} g'_i(\psi_{N,h}^n) \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} \nonumber \\
& + [\eta']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \left[ r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} f_k^n(\psi_{N,h}^n) \right. \nonumber \\
& \left. + \frac{[f_k^n]^2(\psi_{N,h}^n) + |\nabla \psi_h^n|^2}{\mu_0 r} [f']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} - r[jB]_k^{n-1}(\psi_{N,h}^n) \right] g'_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \nonumber \\
& + \eta_k^{n-1}(\psi_{N,h}^n) \left[ r[p'']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} f_k^n(\psi_{N,h}^n) + r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} [f']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \right. \nonumber \\
& \left. + r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} f_k^n(\psi_{N,h}^n) + 2 \frac{f_k^n(\psi_{N,h}^n)}{\mu_0 r} ([f']_k^n(\psi_{N,h}^n))^2 \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} \right. \nonumber \\
& \left. + \frac{2 \nabla \psi}{\mu_0 r} \cdot \nabla \left( \frac{\partial \psi_h}{\partial \psi_j} \right) [f']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} + \frac{[f_k^n]^2(\psi_{N,h}^n) + |\nabla \psi_h^n|^2}{\mu_0 r} [f'']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} \right. \nonumber \\
& \left. + \frac{[f_k^n]^2(\psi_{N,h}^n) + |\nabla \psi_h^n|^2}{\mu_0 r} [f']_k^n(\psi_{N,h}^n) \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} - r[jB]'_{k}^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \right] g'_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \nonumber \\
& + \eta_k^{n-1}(\psi_{N,h}^n) \left[ r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} f_k^n(\psi_{N,h}^n) \right. \nonumber \\
& \left. + \frac{[f_k^n]^2(\psi_{N,h}^n) + |\nabla \psi_h^n|^2}{\mu_0 r} [f']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} - r[jB]_k^{n-1}(\psi_{N,h}^n) \right] g''_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi_j} \frac{\partial \psi_N}{\partial \psi} \nonumber \\
& + \eta_k^{n-1}(\psi_{N,h}^n) \left[ r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} f_k^n(\psi_{N,h}^n) \right. \nonumber \\
& \left. + \frac{[f_k^n]^2(\psi_{N,h}^n) + |\nabla \psi_h^n|^2}{\mu_0 r} [f']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} - r[jB]_k^{n-1}(\psi_{N,h}^n) \right] g'_i(\psi_{N,h}^n) \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} \tag{C.4}
\end{align}
$$
Derivative of the equation for $\psi$ with respect to $f$:

$$ (D_{X_f} \mathbf{e}(\mathbf{X}_{\psi}^n, \mathbf{X}_f^n))_{i,j} = \int_{\Omega^{\text{in}}} \frac{\partial J_i}{\partial f_j} dr dz \quad (C.5) $$

with

$$ \frac{\partial J_i}{\partial f_j} = \frac{1}{\mu_0 r} [g'_i(\psi_{N,h}^n) f_k^n(\psi_{N,h}^n) + [f']_k^n(\psi_{N,h}^n) g_i(\psi_{N,h}^n)] \frac{\partial \psi_N}{\partial \psi} w_i \quad (C.6) $$

Derivative of the equation for $f$ with respect to $f$:

$$ (D_{X_f} \mathbf{d}(\mathbf{X}_{\psi}^n, \mathbf{X}_f^n))_{i,j} = \int_{\Omega^{\text{in}}} \frac{\partial RD_i}{\partial f_j} dr dz \quad (C.7) $$

with

$$ \begin{aligned} \frac{\partial RD_i}{\partial f_j} = & \frac{1}{\Delta t} \frac{g_j(\psi_{N,h}^n)}{r} g_i(\psi_{N,h}^n) + \frac{\psi_h^n - \psi_h^{n-1}}{\Delta t} \frac{g_j(\psi_{N,h}^n)}{r} g'_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \\ & + \eta_k^{n-1}(\psi_{N,h}^n) \left[ r[p']_k^{n-1}(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} g_j(\psi_{N,h}^n) + \frac{2 f_k^n(\psi_{N,h}^n) g_j(\psi_{N,h}^n)}{\mu_0 r} [f']_k^n(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \right. \\ & \left. + \frac{[f_k^n]^2(\psi_{N,h}^n) + |\nabla \psi_h^n|^2}{\mu_0 r} g'_j(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \right] g'_i(\psi_{N,h}^n) \frac{\partial \psi_N}{\partial \psi} \end{aligned} \quad (C.8) $$

and

$$ \frac{\partial \psi_N}{\partial \psi} = \frac{\partial \psi_N}{\partial \psi} \frac{\partial \psi_h}{\partial \psi_j} + \frac{\partial \psi_N}{\partial \psi_a} \frac{\partial \psi_a(\psi_h)}{\partial \psi_j} + \frac{\partial \psi_N}{\partial \psi_b} \frac{\partial \psi_b(\psi_h)}{\partial \psi_j} $$

$$ \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_j} = \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_a} \frac{\partial \psi_a(\psi_h)}{\partial \psi_j} + \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_b} \frac{\partial \psi_b(\psi_h)}{\partial \psi_j} $$

with

$$ \frac{\partial \psi_N}{\partial \psi} = \frac{1}{\psi_b - \psi_a} $$

$$ \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_a} = \frac{1}{(\psi_b - \psi_a)^2} $$

$$ \frac{\partial^2 \psi_N}{\partial \psi \partial \psi_b} = \frac{-1}{(\psi_b - \psi_a)^2} $$
[1] S. Kerboua-Benlarbi, R. Nouailletas, B. Faugeras, E. Nardon, P. Moreau, Magnetic control of WEST plasmas through deep reinforcement learning, in: SOFE 2023, 30th IEEE Symposium On Fusion Engineering, Oxford, UK, 2023.

[2] J. Blum, Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks, Series in Modern Applied Mathematics, Wiley Gauthier-Villars, Paris, 1989.

[3] S. Jardin, Computational methods in plasma physics, Boca Raton, FL: CRC Press/Taylor & Francis, 2010.

[4] H. Grad, J. Hogan, Classical diffusion in a tokamak, Phys. Rev. Lett. 24 (1970) 1337-1340. doi:10.1103/PhysRevLett.24.1337.
URL http://link.aps.org/doi/10.1103/PhysRevLett.24.1337

[5] H. Grad, P. Hu, D. Stevens, Adiabatic evolution of plasma equilibrium, Proc. Nat. Acad. Sci. USA 72 (1975) 3789-3793.

[6] F. Hinton, R. Hazeltine, Theory of plasma transport in toroidal confinement systems, Rev. Mod. Phys. 48 (1976) 239-308. doi:10.1103/RevModPhys.48.239.
URL http://link.aps.org/doi/10.1103/RevModPhys.48.239

[7] R. Byrne, H. Klein, G2M, a two-dimensional diffusion time scale tokamak code, Journal of Computational Physics 26 (3) (1978) 352-378. doi:10.1016/0021-9991(78)90075-X.
URL http://www.sciencedirect.com/science/article/pii/S002199917890075X

[8] J. Blum, J. Le Foll, B. Thooris, The self-consistent equilibrium and diffusion code SCED, Computer Physics Communications 24 (3) (1981) 235-254. doi:https://doi.org/10.1016/0010-4655(81)90149-1.
URL https://www.sciencedirect.com/science/article/pii/S0010465581901491

[9] A. Turnbull, R. Storer, A plasma resistive diffusion model, Journal of Computational Physics 50 (3) (1983) 409-435. doi:https://doi.org/10.1016/0021-9991(83)90104-3.
772 URL https://www.sciencedirect.com/science/article/pii/
0021999183901043

774 [10] J. Blum, J. Le Foll, Plasma equilibrium evolution at the resistive diffusion timescale, Computer Physics Reports 1 (7-8) (1984) 465-494.
doi:10.1016/0167-7977(84)90013-3.

776 URL http://www.sciencedirect.com/science/article/pii/
0167797784900133

780 [11] S. C. Jardin, N. Pomphrey, J. DeLucia, Dynamic modeling of transport and positional control of tokamaks, J. Comput. Phys. 66 (1986) 481-507.
doi:10.1016/0021-9991(86)90077-X.

782 URL http://dl.acm.org/citation.cfm?id=9255.9268

784 [12] J. Crotinger, CORSICA; a comprehensive simulation of toroidal magnetic-fusion devices, Technical report ucll-id-126284, Lawrence Livermore National Laboratory (1997).

786 URL http://wormhole.ucllnl.org/software/corsica/finalreport.pdf

788 [13] S. H. Kim, J. F. Artaud, V. Basiuk, V. Dokuka, R. R. Khayrutdinov, J. B. Lister, V. E. Lukash, Full tokamak discharge simulation of ITER by combining DINA-CH and CRONOS, Plasma Physics and Controlled Fusion 51 (10) (2009) 105007.

792 URL http://stacks.iop.org/0741-3335/51/i=10/a=105007

794 [14] M. Honda, Simulation technique of free-boundary equilibrium evolution in plasma ramp-up phase, Computer Physics Communications 181 (9) (2010) 1490-1500.
doi:http://10.1016/j.cpc.2010.04.014.

796 [15] E. Fable, C. Angioni, F. J. Casson, D. Told, A. A. Ivanov, F. Jenko, R. M. McDermott, S. Y. Medvedev, G. V. Pereverzev, F. Ryter, W. Treutterer, E. Viezzer, the ASDEX Upgrade Team, Novel free-boundary equilibrium and transport solver with theory-based models and its validation against ASDEX Upgrade current ramp scenarios, Plasma Physics and Controlled Fusion 55 (12) (2013) 124028.

800 doi:10.1088/0741-3335/55/12/124028.
802 URL https://dx.doi.org/10.1088/0741-3335/55/12/124028
[16] V. Parail, R. Albanese, R. Ambrosino, J.-F. Artaud, K. Besseghir, M. Cavinato, G. Corrigan, J. Garcia, L. Garzotti, Y. Gribov, F. Imbeaux, F. Koechl, C. Labate, J. Lister, X. Litaudon, A. Loarte, P. Maget, M. Mattei, D. McDonald, E. Nardon, G. Saibene, R. Sartori, J. Urban, Self-consistent simulation of plasma scenarios for ITER using a combination of 1.5d transport codes and free-boundary equilibrium codes, Nuclear Fusion 53 (11) (2013) 113002. doi:10.1088/0029-5515/53/11/113002.
URL https://dx.doi.org/10.1088/0029-5515/53/11/113002

[17] G. L. Falchetto, D. Coster, R. Coelho, B. Scott, L. Figini, D. Kalupin, E. Nardon, L. Alves, J. Artaud, V. Basiuk, J. Bizarro, C. Boulbe, A. Dinklage, D. Farina, B. Faugeras, J. Ferreira, A. Figueiredo, P. Huynh, F. Imbeaux, I. Ivanova-Stanik, T. Jonsson, H.-J. Klingshirn, C. Konz, A. Kus, N. Marushchenko, E. Nardon, S. Nowak, G. Pereverzev, M. Owsiak, E. Poli, Y. Peysson, R. Reimer, J. Signoret, O. Sauter, R. Stankiewicz, P. Strand, I. Voitsekhovitch, E. Westerhof, T. Zok, W. Zwingmann, I.-T. contributors, the ASDEX Upgrade Team, J.-E. Contributors, The European Integrated Tokamak Modelling (ITM) effort: achievements and first physics results, Nuclear Fusion 54 (4) (2014) 043018. doi:http://dx.doi.org/10.1088/0029-5515/54/4/043018.

[18] V. Ostuni, J. Artaud, G. Giruzzi, E. Joffrin, H. Heumann, H. Urano, Tokamak discharge simulation coupling free-boundary equilibrium and plasma model with application to JT-60SA, Nuclear Fusion 61 (2) (2021) 026021. doi:10.1088/1741-4326/abcdb7.
URL https://dx.doi.org/10.1088/1741-4326/abcdb7

[19] H. Heumann, A Galerkin method for the weak formulation of current diffusion and force balance in tokamak plasmas, Journal of Computational Physics 442 (2021) 110483. doi:https://doi.org/10.1016/j.jcp.2021.110483.
URL https://www.sciencedirect.com/science/article/pii/S0021999121003788

[20] J. Taylor, Evolution of resistive plasma in an axisymmetric torus, Ppn-10/75, Theory division, Culham Laboratory (1975).

[21] A. Elarif, B. Faugeras, F. Rapetti, Tokamak free-boundary plasma equilibrium computation using finite elements of class C0 and C1 within a
mortar element approach, J. Computational Physics 439 (2021) 110388,
doi:https://doi.org/10.1016/j.jcp.2021.110388.

840

[22] C. Boulbe, B. Faugeras, G. Gros, F. Rapetti, Tokamak free-boundary plasma equilibrium computations in presence of non-linear materials, J. Scientific Computing 96:42, https://hal.sciences/hal-03423469 (2023). doi:10.1007/s10915-023-02265-8.

842

844

[23] B. Faugeras, An overview of the numerical methods for tokamak plasma equilibrium computation implemented in the NICE code, Fusion Eng. Design 160 (2020) 112020, https://hal.archives-ouvertes.fr/hal-02955053. doi:10.1016/j.fusengdes.2020.112020.

846

848

[24] J. Bucalossi, M. Missirlian, P. Moreau, F. Samaille, E. Tsitrone, D. Van Houtte, T. Batal, C. Bourdelle, M. Chantant, Y. Corre, X. Courtois, L. Delpech, L. Doceul, D. Douai, H. Dougnac, F. Faïsse, C. Fenzi, F. Ferlay, M. Firdaouss, L. Gargiulo, P. Garin, C. Gil, A. Grosman, D. Guilhem, J. Gunn, C. Hernandez, D. Keller, S. Laroque, F. Leroux, M. Lipa, P. Lotte, A. Martinez, O. Meyer, F. Micolon, P. Mollard, E. Nardon, R. Nouailletas, A. Pilia, M. Richou, S. Salasca, J.-M. Travère, The WEST project: Testing ITER divertor high heat flux component technology in a steady state tokamak environment, Fusion Engineering and Design 89 (7-8) (2014) 907-912. doi:10.1016/j.fusengdes.2014.01.062.
URL https://linkinghub.elsevier.com/retrieve/pii/S0920379614000635

850

852

854

856

858

860

[25] C. Bourdelle, J. Artaud, V. Basiuk, M. Bécoulet, S. Brémond, J. Bucalossi, H. Bufferand, G. Ciraolo, L. Colas, Y. Corre, X. Courtois, J. Decker, L. Delpech, P. Devynck, G. Dif-Pradalier, R. Doerner, D. Douai, R. Dumont, A. Ekedahl, N. Fedorczak, C. Fenzi, M. Firdaouss, J. Garcia, P. Ghendrih, C. Gil, G. Giruzzi, M. Goniche, C. Grisolia, A. Grosman, D. Guilhem, R. Guirlet, J. Gunn, P. Hennequin, J. Hillairet, T. Hoang, F. Imbeaux, I. Ivanova-Stanik, E. Joffrin, A. Kallenbach, J. Linke, T. Loarer, P. Lotte, P. Maget, Y. Marandet, M. Mayoral, O. Meyer, M. Missirlian, P. Mollard, P. Monier-Garbet, P. Moreau, E. Nardon, B. Pégourié, Y. Peysson, R. Sabot, F. Saint-Laurent, M. Schneider, J. Travère, E. Tsitrone,

862

864

866

868

870

872
S. Vartanian, L. Vermare, M. Yoshida, R. Zagorski, JET Contributors, WEST Physics Basis, Nuclear Fusion 55 (6) (2015) 063017. doi:10.1088/0029-5515/55/6/063017.

URL https://iopscience.iop.org/article/10.1088/0029-5515/55/6/063017

[26] J. Bucalossi, J. Achard, O. Agullo, T. Alarcon, L. Allegretti, H. Ancher, G. Antar, S. Antusch, V. Anzallo, C. Arnas, D. Arranger, J. Artaud, M. Aumeunier, S. Baek, X. Bai, J. Balbin, C. Balorin, T. Barbui, A. Barbuti, J. Barlerin, V. Basiuk, T. Batal, O. Baulaigue, A. Bec, M. Bécoulet, E. Benoit, E. Benard, J. Benard, N. Bertelli, E. Bertrand, P. Beyer, J. Bielecki, P. Bienvenu, R. Bisson, V. Bobkov, G. Bodner, C. Bottereau, C. Bouchand, F. Bouquey, C. Bourdelle, J. Bourg, S. Brezinsek, F. Brochard, C. Brun, V. Bruno, H. Bufferand, A. Bureau, S. Burles, Y. Camenen, A. Hakola, Operating a full tungsten actively cooled tokamak: Overview of WEST first phase of operation, Nuclear Fusion 62 (4) (Feb. 2022). doi:10.1088/1741-4326/ac2525.

[27] R. Nouailletas, E. Nardon, P. Moreau, C. Reux, T.-P.-H. Truong, WEST magnetic control, in: 2019 IEEE 58th Conference on Decision and Control (CDC), 2019, pp. 3214-3219. doi:10.1109/CDC40024.2019.9029207.

[28] H. Heumann, J. Blum, C. Boulbe, B. Faugeras, G. Selig, J.-M. Ané, S. Brémond, V. Grangirard, P. Hertout, E. Nardon, Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES++: computational methods and applications, J. Plasma Physics 81 (3) (2015) 905810301. doi:https://doi.org/10.1017/S0022377814001251.

[29] J. Blum, H. Heumann, E. Nardon, X. Song, Automating the design of tokamak experiment scenarios, J. Comp. Phys. 394 (2019) 594-614.

[30] H. Grad, H. Rubin, Hydromagnetic equilibria and force-free fields, Proceedings of the 2nd UN Conf. on the Peaceful Uses of Atomic Energy 31 (1958) 190.

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