Blaise Faugeras
email: blaise.faugeras@univ-cotedazur.fr

An overview of the numerical methods for tokamak plasma equilibrium computation implemented in the NICE code

Keywords: tokamak, plasma equilibrium, equilibrium reconstruction, toroidal harmonics, finite element method, Newton method, inverse problem, PDE-constrained optimization PACS: 02.03.Zz, 02.60.-x, 52.55.-s, 52.55.Fa, 52.65.-y

The code NICE (Newton direct and Inverse Computation for Equilibrium) enables to solve numerically several problems of plasma free-boundary equilibrium computations in a tokamak: plasma free-boundary only reconstruction and magnetic measurements interpolation, full free-boundary equilibrium reconstruction from magnetic measurements and possibly internal measurements (interferometry, classical linear approximation polarimetry or Stokes model polarimety, Motional Stark Effect and pressure), direct and inverse, static and quasi-static free-boundary equilibrium computations.

. The strength of NICE is to gather in a single finite element framework different equilibrium computation modes. It makes intensive use of Newton method and Sequential Quadratic Programming method to solve non linear problems.

NICE is used routinely for WEST tokamak operation. It is also adapted to the IMAS (ITER Modelling and Analysis Suite) format which makes it usable on many different fusion tokamak reactors.

In this document we give a general overview of the numerical methods implemented in NICE as well as a number of computation examples.

Introduction

The aim of this paper is to provide a reference document on the different operation modes and numerical methods implemented in the code NICE. NICE stands for Newton direct and Inverse Computation for Equilibrium. Its development started in 2017 with the original goal of developing an equilibrium reconstruction code able to use polarimetry Stokes vector measurements. It demanded the computation of derivatives of certain quantities which did not easily fit in the equilibrium reconstruction code EQUINOX [START_REF] Blum | Reconstruction of the equilibrium of the plasma in a tokamak and identification of the current density profile in real time[END_REF] developed by the author and others. This led to the development and implementation of the optimal control algorithm, presented in this document, with a first step presented in [START_REF] Faugeras | Optimal control of a coupled partial and ordinary differential equations system for the assimilation of polarimetry stokes vector measurements in tokamak free-boundary equilibrium reconstruction with application to ITER[END_REF], which can naturally deal with the non-linear Stokes equations. Moreover it soon appeared appealing to unify in a single performant C++ code the different functionalities from 3 older codes, VacTH [START_REF] Faugeras | 2D interpolation and extrapolation of discrete magnetic measurements with toroidal harmonics for equilibrium reconstruction in a Tokamak[END_REF], EQUINOX [START_REF] Blum | Reconstruction of the equilibrium of the plasma in a tokamak and identification of the current density profile in real time[END_REF] and CEDRES++ [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: computational methods and applications[END_REF] which share many common features. These goals are now achieved and on the way, improvements to the original methods and new functionalities were added.

Several problems of equilibrium computation can be addressed by NICE and are detailed in this document: plasma boundary only reconstruction from magnetic measurements (Section 3). The reference paper for this mode is [START_REF] Faugeras | 2D interpolation and extrapolation of discrete magnetic measurements with toroidal harmonics for equilibrium reconstruction in a Tokamak[END_REF]. In addition to the method proposed in this reference the use of an original regularization term which proved to be efficient is proposed in this document. full equilibrium reconstruction from magnetic measurements and possibly internal measurements (Section 4). NICE uses Sequential Quadratic Programming (SQP) algorithms to solve partial differential equation (PDE) constrained optimization problems and in particular the inverse problem of equilibrium reconstruction. The default algorithm used in NICE is a quasi SQP with reduced Hessian algorithm (which we denote QSQP in this document). Its implementation demands the computation of derivatives of quantities with respect to the state and control variables and enables to use polarimetry Stokes vector measurements [START_REF] Faugeras | Optimal control of a coupled partial and ordinary differential equations system for the assimilation of polarimetry stokes vector measurements in tokamak free-boundary equilibrium reconstruction with application to ITER[END_REF][START_REF] Faugeras | On the identification of the electron temperature profile from polarimetry Stokes vector measurements in tokamak free-boundary equilibrium reconstruction[END_REF]. The code also incorporates the original optimization algorithm from EQUINOX [START_REF] Blum | Reconstruction of the equilibrium of the plasma in a tokamak and identification of the current density profile in real time[END_REF] which appears to be an approximation of the default QSQP algorithm used in NICE. A sensitivity method for error bar computations on every reconstructed equilibrium quantities is also implemented and detailed in Section 4.3.2. direct and inverse, static and quasi-static free-boundary equilibrium computations (Section 5). The reference papers for these implemented modes are [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: computational methods and applications[END_REF] and [START_REF] Blum | Automating the design of tokamak experiment scenarios[END_REF]. The term direct refers to the resolution of the equilibrium equation whereas the term inverse refers to the inverse problem consisting in finding the currents in the poloidal field (PF) coils which give a desired plasma shape in the static case or the voltages in the suppliers which give a desired evolution of the plasma shape in the quasi-static case. The difference between static and quasi-static is explained in Section 2.1

The strengh of NICE is to gather in a single unified framework different functionalities or equilibrium computation modes. It makes intensive use of Newton method and SQP method to solve non linear problems.

NICE is used routinely for the WEST tokamak operation [START_REF] Moreau | The new magnetic diagnostics in the west tokamak[END_REF]. NICE is adapted to the IMAS (ITER Modelling and Analysis Suite [START_REF] Imbeaux | Design and first applications of the ITER integrated modelling and analysis suite[END_REF][START_REF] Pinches | Progress in the ITER Integrated Modelling Programme and the Use and Validation of IMAS within the ITER Members[END_REF]) format which makes it usable on many different tokamaks [START_REF] Coelho | Plasma equilibrium reconstruction of JET discharges using the IMAS modelling infrastructure[END_REF]. Many of the numerical examples provided in this document are obtained using data in the IMAS format.

The document is organized as follows. Section 2 provides a brief description of the equilibrium equations used (Section 2.1), of their finite element discretization (Section 2.2) and of the iterative Newton and QSQP algorithm used (Section 2.3). This section enables to set matrix notations which are very close to what is used in the code itself. Section 3 describes the plasma boundary only reconstruction mode of NICE. This first mode does not rely on finite element discretization but uses a decomposition of the poloidal flux in toroidal harmonics. All other computation modes of NICE use finite elements. Section 4 presents the equilibrium reconstruction mode and Section 5 the direct and inverse equilibrium computation modes. For each mode numerical examples are provided.

Tokamak free-boundary plasma equilibrium

This section presents the general equations which are dealt with in NICE. The finite element discretization is also presented which enables us to set the notations for different matrix and vector operators which are used in the different modes of NICE. Then we present in a generic framework the QSQP algorithm which will then be applied to the inverse problems of the next sections.

Modelization

The equations which govern the equilibrium of a plasma. in the presence of a magnetic field in a tokamak are on the one hand Maxwell's equations satisfied in the whole of space (including the plasma):

∇ • B = 0, ∇ × ( B µ ) = j, (1)
and on the other hand the equilibrium equation for the plasma itself

∇p = j × B, (2)
where B is the magnetic field, µ is the magnetic permeability, p is the kinetic pressure and j is the current density. These equations are sufficient for the static modelization of the plasma equilibrium. In this modelization the current densities in the PF coils are given. In case of quasi-static plasma evolution modelization these equations are augmented with Faraday's law and Ohm's law in the PF coils and passive structures. The current density in these structures is then computed as a function of the input voltages in the PF circuits and/or of the time derivative of the poloidal flux (defined below) as shown in Eqs. [START_REF] Freidberg | Ideal Magnetohydrodynamics[END_REF] and [START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks[END_REF] below. We refer to standard text books (e.g. [START_REF] Freidberg | Ideal Magnetohydrodynamics[END_REF][START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks[END_REF][START_REF] Wesson | Tokamaks[END_REF][START_REF] Goedbloed | Principles of magnetohydrodynamics: with applications to laboratory and astrophysical plasmas[END_REF][START_REF] Jardin | Computational methods in plasma physics[END_REF]) and to [START_REF] Blum | Reconstruction of the equilibrium of the plasma in a tokamak and identification of the current density profile in real time[END_REF][START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: computational methods and applications[END_REF] for details of the derivation and only state the needed equations in what follows.

Introducing a cylindrical coordinate system (e r , e φ , e z ) (r = 0 is the major axis of the tokamak torus) and assuming axial symmetry equations (1) and (2) reduce to the following equation for the poloidal flux ψ(r, z) in the poloidal plane

Ω ∞ = (0, ∞) × (-∞, ∞): -∆ * ψ = j φ , (3)
where j φ is the toroidal component of j, and the second order elliptic differential operator ∆ * is defined by

∆ * ψ := ∂ r 1 µ(ψ)r ∂ r ψ + ∂ z 1 µ(ψ)r ∂ z ψ := ∇ • 1 µ(ψ)r ∇ψ . ( 4
)
Here ∇ is the 2D operator in the (r, z)-plane and µ(ψ) is the magnetic permeability. It is equal to the constant permeability of vacuum µ 0 everywhere except in the possibly existing iron parts of the tokamak (see Fig. 1) where it is a given function of ψ, µ(ψ) = µ f (|∇ψ| 2 r -2 ).

The magnetic field can be decomposed in poloidal and toroidal components

B = B p + B φ , B p = 1 r [∇ψ × e φ ], B φ = B φ e φ = f r e φ , (5)
where f is the diamagnetic function. Equation [START_REF] Faugeras | On the identification of the electron temperature profile from polarimetry Stokes vector measurements in tokamak free-boundary equilibrium reconstruction[END_REF] shows that the magnetic surfaces are generated by the rotation of the iso-flux lines around the axis of the torus.

The toroidal component of the current density. j φ is zero everywhere outside the plasma domain and the poloidal field coils (and possibly the passive structures). The different sub-domains of the poloidal plane of a schematic tokamak (see Figure 1) as well as the corresponding expression for j φ are described below:

-Ω f is the domain of ferromagnetic iron structures where the permeability µ is not constant.

-Ω L is the domain accessible to the plasma. Its boundary is the limiter Γ L .

-Ω p is the plasma domain where equations ( 2) and (1) imply that p and f are constant on each magnetic surface i.e. p = p(ψ) and f = f (ψ). One then deduces the so-called Grad-Shafranov equilibrium equation in the plasma [16, [START_REF] Shafranov | On magnetohydrodynamical equilibrium configurations[END_REF][START_REF] Lüst | Axialsymmetrische magnetohydrodynamische Gleichgewichtskonfigurationen[END_REF]]

-∆ * ψ = rp (ψ) + 1 µ 0 r (f f )(ψ). ( 6
)
The right-hand side of ( 6) is the toroidal component j φ of the plasma current density.

The plasma domain depends on ψ and is unknown, Ω p = Ω p (ψ). We have to deal with a free-boundary problem. This domain is defined by its boundary which is the outermost closed ψ iso-contour contained within the limiter region Ω L . The plasma can either be limited if this iso-contour is tangent to the limiter Γ L or defined by the presence of an X-point (as in Figure 1). In the latter case the plasma domain is bounded by the magnetic separatrix. Functions p and f f are zero outside Ω p .

The current density is non-linear in ψ due to the non-linear functions p and f f and the definition of the plasma domain Ω p (ψ). While Ω p (ψ) is fully determined for a given ψ, the two functions p and f f are not determined in this modelization. It is the goal of the inverse equilibrium reconstruction problem to determine them. For now let us consider that we are given two functions A(x) and B(x) defined on [0, 1] such that in the plasma domain Ω p (ψ)

j φ = λ( r r 0 A(ψ N ) + r 0 r B(ψ N )). ( 7
)
Here r 0 is the characteristic major radius of the tokamak vacuum chamber and λ is a scaling coefficient. The normalized poloidal flux ψ N (r, z) is

ψ N (r, z) = ψ(r, z) -ψ a (ψ) ψ b (ψ) -ψ a (ψ) . ( 8
)
with ψ a and ψ b being the flux values at the magnetic axis and at the boundary of the plasma:

ψ a (ψ) := ψ(r a (ψ), z a (ψ)), ψ b (ψ) := ψ(r b (ψ), z b (ψ)), (9)
with (r a (ψ), z a (ψ)) the magnetic axis where ∇ψ = 0 and ψ has a local extremum in Ω L , and (r b (ψ), z b (ψ)) the coordinates of the point that determines the plasma boundary. The point (r b , z b ) is either an X-point of ψ or the contact point with the limiter Γ L .

-Domains Ω c i represent the poloidal field coils carrying currents

j φ = I i S i , (10)
where S i is the section area of the coil and I i is the current flowing in the coil. In the static modelization I i is a constant whereas in the quasi-static case it relates to voltages v j (t) in the suppliers and to self and mutual induction

I i (t) = Ns j=1 S ij v j (t) + Nc j=1 R ij 1 S j Ωc j ∂ψ ∂t drdz, 1 ≤ i ≤ N c (11)
via electric circuit equations (see [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: computational methods and applications[END_REF] Appendix A and [START_REF] Blum | Automating the design of tokamak experiment scenarios[END_REF] Appendix B). Here N s is the number of suppliers, N c is the number of coils, matrix S has dimension N c × N s and matrix R has dimension N c × N c . In the simplest case where each circuit is composed of one coil and one voltage supplier, S and R are diagonal with

S ii = n i R i and R ii = - 2πn 2 i R i
, n i being the number of turns and R i the resistance.

-Domain Ω ps represents passive structures. The expression for j φ in these domains depends on the modelization considered. For static equilibrium j φ = 0 is usually considered whereas for quasi-static evolution of the equilibrium

j φ = - σ r ∂ψ ∂t ( 12
)
where σ is the conductivity of the passive structure.

Equilibrium equation [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: computational methods and applications[END_REF]. can either be considered in the whole poloidal plane Ω ∞ with ψ = 0 as boundary condition at infinity and on the r = 0 axis, or on a restricted bounded domain (whose boundary can be viewed as a measurement contour Γ and is illustrated by the green dashed line on Figure 1). Both type of problems are addressed in NICE, the first one for the direct and inverse static and quasi-static modes and the second one for the equilibrium reconstruction mode.

FEM discretized operators and derivatives

For the equilibrium reconstruction mode as well as the direct and inverse modes of NICE equation ( 3) is discretized using a P1 finite element method based on a triangular mesh [START_REF] Blum | Plasma equilibrium evolution at the resistive diffusion timescale[END_REF][START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks[END_REF][START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: computational methods and applications[END_REF]. In this Section we provide some details on the discretized operators common to different modes of NICE as well as on the computation of their derivatives which are used in the Newton and SQP algorithms implemented for the resolution of non-linear direct and inverse problems. More operators specific to a given mode are presented later in the document.

There are two possible choices for the spatial domain to be triangulated. It can either be a restricted domain enclosed in the measurement contour Γ or a sufficiently large semi-circle centered at the origin and containing the whole geometry of the tokamak. The first case is called the bounded domain case (see Figure 2) and the second the ABB domain case (see Figures 3 and4). ABB stands for Albanase, Blum, de Barbieri [START_REF] Albanese | On the solution of the magnetic flux equation in an infinite domain[END_REF] who first introduced the boundary integral method on the semi-circle used to take into account conditions at infinity. We refer to [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: computational methods and applications[END_REF] for its precise expression and to [START_REF] Grandgirard | Modélisation de l'équilibre d'un plasma de tokamak[END_REF]Chapter 2.4] for details on its derivation. Alternative approaches for the incorporation of boundary conditions at infinity were more recently presented in [START_REF] Faugeras | FEM-BEM coupling methods for tokamak plasma axisymmetric free-boundary equilibrium computations in unbounded domains[END_REF].

Let us denote by λ i (r, z) the N Lagrangian basis functions associated to the inner vertices of the mesh and by λ d,i the N b ones associated to boundary vertices where Dirichlet conditions are applied (r = 0 in the ABB case or Γ in the bounded domain case).

Approximating the poloidal flux by

ψ h = N i=1 ψ i λ i (r, z) + N b i=1 ψ d,i λ d,i (r, z)
and testing equation (3) against every inner basis function leads to the definition of the following matrices and vectors: state variable. ψ denotes the vector of N finite elements coefficients ψ i elliptic operator. split into a linear and a non linear part (if iron is present), Aψ + A µ (ψ). Matrix A is of size N × N and coefficients

A ij = Ω-Ω f 1 µ 0 r ∇λ i ∇λ j drdz
involving the scalar product of the gradients of basis functions i and j. Vector A µ (ψ) of dimension N has coefficients

A µ,i (ψ) = Ω f 1 µ(ψ h )r ∇λ i ∇ψ h drdz
The integrals are computed as sums of integrals over the triangles T of the mesh on each of which the barycentric quadrature rule is used for the integral approximation. The Jacobian A µ,ψ (ψ) is computed as

[A µ,ψ ] ij = Ω f ∇λ i ∇λ j µ(ψ h )r - 2µ f (|∇ψ h | 2 r -2 ) µ 2 f (|∇ψ h | 2 r -2 )
A d,ij = Ω 1 µ 0 r ∇λ i ∇λ d,j drdz
This matrix is only used in the bounded domain case since for the ABB domain case, ψ d = 0.

In the ABB domain case the boundary integral method on the semi-circle leads to the definition of a matrix C of size N × N (see [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: computational methods and applications[END_REF]). plasma current density. J(ψ, λ, u A , u B ) vector of size N and coefficient

J(ψ, λ, u A , u B ) i = Ωp(ψ h ) λ( r r 0 A(ψ N ) + r 0 r B(ψ N ))λ i drdz
Functions A and B are decomposed in a basis of functions φ i defined on [0, 1].

In NICE piecewise linear and cubic splines functions are implemented.

A(x) = N A i=1 u Ai φ i (x), B(x) = N B i=1 u Bi φ i (x),
and u A , u B denote the vectors of degrees of freedom of A and B in the chosen decomposition basis. In the direct and inverse modes of NICE functions A and B can also have a parametric representation of the type

A(x) = β(1 -x α ) γ , B(x) = (1 -β)(1 -x α ) γ
The mesh does not resolve the boundary of plasma domain Ω p (ψ h ) and a quadrature rule needs to be specified for integrals over intersections of mesh triangles with the plasma domain T ∩ Ω p (ψ h ). Barycenter quadrature is used but the quadrature point and weight depend non linearly on ψ h . This needs to be taken into account in the computation of the Jacobian J ψ (ψ, λ, u A , u B ). This differentiation is not straightforward and technical details are given in [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: computational methods and applications[END_REF]. The computation of Jacobians J u A (ψ, λ, u A , u B ) and J u B (ψ, λ, u A , u B ) which are needed for equilibrium reconstruction on the other hand does not rise any particular difficulty.

Finally we will use the notation

J p (ψ, u A , u B ) = Ωp(ψ h ) ( r r 0 A(ψ N ) + r 0 r B(ψ N ))drdz (13)
such that the total plasma current is computed as

I p = λJ p (ψ, u A , u B )
PF coils. matrix L of size N × N c and vector u I of size N c holding the currents I i of the N c coils contained in the domain under consideration (all coils in the ABB domain case but not in the bounded domain case), with coefficients

L ij = 1 S j Ωc j λ i drdz
such that the PF coils term in the static case is Lu I and in the quasi-static evolution case is LSv + LRL T dψ dt with v the vector of the N s voltages passive structures. matrix J ps of size N × N , associated to the inductive terms in the passive structures with coefficients

[J ps ] ij = - Nps k=1 Ωps k σ k r λ i λ j drdz
where N ps is the number of passive structure and σ k their conductivity.

Newton and PDE-constrained optimization algorithms

With the matrix and vector operators defined in the previous section different direct and inverse problems can be formulated in NICE and are detailed in the remaining sections of this document. Here we present the main algorithms using a generic formulation. We will define different model equations e, state y and control u variables as well as cost functions J(y, u). They will differ from one mode to the other in the code but the resolution algorithms are the same.

The direct problem formulation. is: find y such that e(y, u) = 0, where u is here a constant parameter.

The inverse problem formulation. is: min y,u J(y, u) under the constraint e(y, u) = 0, Throughout this document several explicit formulations are given for e and J. As a starter and to make it more concrete let us consider the static equilibrium problem in the ABB domain detailed in Section 5. The sate variable is y := ψ, the control variable is u := u I the currents in the coils and p and f f are given (thus u A , u B and λ are fixed parameters). The direct model in ABB domain is

e(y, u) := (A + C)ψ + A µ (ψ) -J(ψ, λ, u A , u B ) -Lu I = 0
and the cost function is

J(y, u) = 1 2 ||Ky|| 2 + 1 2 ||Ru|| 2
with a quadratic misfit term, imposing a levelset of the poloidal flux to go through a set of desired points, and a regularization term (see Section 5).

In NICE direct problems are solved thanks to Newton's method. If we denote the Newton increment by ∆y = -e -1 y (y n , u)e(y n , u) this iterative algorithm reads:

y n+1 = y n + ∆y
Inverse problems are solved thanks to a quasi-SQP method with reduced Hessian (QSQP). SQP methods are well documented [START_REF] Nocedal | Numerical optimization[END_REF][START_REF] Hinze | Mathematical Modelling: Theory and Applications[END_REF] and a clear summary is found in [6, Appendix A]. An SQP method can be seen as a Newton method to solve the non-linear system given by the fisrt order optimality condition for the Lagrangian of the PDE-constrained optimization problem L(y, u, p) = J(y, u) + (p, e(y, u))

where vector p of size N is the adjoint variable or Lagrange multiplier for the constraint of the model. This first method also known as Newton-Lagrange method is implemented in NICE and was used in [START_REF] Faugeras | Optimal control of a coupled partial and ordinary differential equations system for the assimilation of polarimetry stokes vector measurements in tokamak free-boundary equilibrium reconstruction with application to ITER[END_REF] for equilibrium reconstruction. It was also already used in [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: computational methods and applications[END_REF] for the static inverse problem (presented below in Section 5). The drawback of this method is that it demands the resolution of linear systems of size (2N + N u ), where N is the dimension of the state and adjoint variables and N u the dimension of the control. When N is large this can become time consumming and an excellent alternative is the SQP formulation with reduced Hessian which is intensively used in NICE and was already used in [START_REF] Blum | Automating the design of tokamak experiment scenarios[END_REF] for the quasi-static inverse problem (see Section 5.2). It is the following two-step iterative algorithm.

quasi SQP with reduced Hessian algorithm (QSQP):.

control variable update step

M (u n+1 -u n ) = -m (15)
2. state variable update step

y n+1 = y n + ∆y + S(u n+1 -u n ) (16)
where ∆y = -e -1 y (y n , u n )e(y n , u n ), ( 17)

S = -e -1 y (y n , u n )e u (y n , u n ), (18)
M = J uu (y n , u n ) + S T J yu (y n , u n ) + J uy (y n , u n )S + S T J yy (y n , u n )S (19)
and

m = J T u (y n , u n ) + S T J T y (y n , u n ) + J uy (y n , u n )∆y + S T J yy (y n , u n )∆y (20)
At each iteration this algorithm demands the resolution of N u + 1 linear systems ( 17)-( 18) of size N involving the same matrix e y (y n , u n ) with different righ-hand sides which can be done very efficiently and of one smaller linear system of size N u in [START_REF] Jardin | Computational methods in plasma physics[END_REF].

One can notice that second order derivatives appear only for the cost function and not for the model equation in [START_REF] Blum | Plasma equilibrium evolution at the resistive diffusion timescale[END_REF] and [START_REF] Albanese | On the solution of the magnetic flux equation in an infinite domain[END_REF]. This is due to the fact that on the one hand the model equation in our applications is affine in the control variable, u, therefore has 0 second order derivative, and on the other hand we want to avoid the expensive computation of other second order derivatives of J(ψ, λ, u A , u B ) and neglect them. Neglecting these second order derivatives of the model equation also allows to avoid having to compute the adjoint state variable. Moreover in some cases where in the cost function the observation operator is non-linear we also neglect some of the second order derivatives terms in J yu , J uy and J yy (see equilibrium reconstruction with internal measurements in Section 4.2). NICE uses a quasi SQP algorithm and not an exact SQP algorithm. The excellent convergence properties of the code allow us to think that these modifications are not an issue for the applications specific to tokamaks dealt with in NICE.

It should also be noted that if one neglects the first order derivative J ψ term in the computation of model derivative e y this algorithm applied to the equilibrium reconstruction problem of Section 4 is exactly the original algorithm from the code EQUINOX. In the original EQUINOX paper [START_REF] Blum | Reconstruction of the equilibrium of the plasma in a tokamak and identification of the current density profile in real time[END_REF] it is described as a fixed-point method to solve the model equation in ψ for which the control variable u is updated at each iteration. It appears clearly that this method is another (more approximated) quasi SQP with reduced Hessian method. In our experience however this method does not always demonstrates the same quality and speed of convergence as the one used in NICE which involves more derivative terms.

Cost functions notation

Throughout this document cost functions of the following type appear

J(ψ, u) = 1 2 ||M(ψ, u) -m|| 2 + 1 2 ||Ru|| 2 (21)
The first term is the data misfit term and the second the regularization term. Each measurement data m i is a given quantity which goes together with a given absolute measurement error, abserr, and an relative measurement error, relerr, and we form the quantity σ i = max(abserr, relerr|m i |). All vectorial quantities related to measurements which appear in cost functions are normalized quantities. The vector m of measurements has components m i = m i /σ i and the observation operator M(ψ, u) computes from the model the equivalent of these normalized measured quantities. In NICE absolute and relative errors can either be provided together with the used data set (for example these quantities should be provided in IMAS) or set by the user through input code parameters.

Regularization terms are needed due to the ill-posedness of the inverse problems under consideration and they guarantee smoothness of the identified quantities. Regularization and penalization parameters are included in R matrices and are provided by the user through input code parameters.

NICE magnetic measurements interpolation and plasma boundary reconstruction mode

The goal of this section is to present the magnetic measurements interpolation and plasma boundary reconstruction mode of NICE.

Method

Toroidal harmonics are used to represent the poloidal flux in an annular domain surrounding the plasma. In such a domain ψ satisfies -∆ * ψ = 0 and can be expanded in a series of toroidal harmonics (see [START_REF] Faugeras | 2D interpolation and extrapolation of discrete magnetic measurements with toroidal harmonics for equilibrium reconstruction in a Tokamak[END_REF] and references therein).

This method allows a fast fit to magnetic measurements and the poloidal field function can then be used to compute the plasma boundary. It allows also to compute precisely the poloidal flux and its gradient on any contour surrounding the plasma and thus on Γ the measurements contour in NICE bounded domain case (green dashed line in Figure 1). These Cauchy conditions on Γ, (g = ψ, h = ∂ n ψ), can then be used in a second step to reconstruct the full non-linear free-boundary Grad-Shafranov equilibrium problem inside the domain Ω limited by Γ. This is presented in the following Section 4.

Below we propose to use a simple additional regularization term which was not originally present in [START_REF] Faugeras | 2D interpolation and extrapolation of discrete magnetic measurements with toroidal harmonics for equilibrium reconstruction in a Tokamak[END_REF]. This regularization operates directly on the toroidal harmonics expansion in order to obtain a smooth representation of flux even faraway from the measurements contour Γ inside the domain Ω.

The toroidal coordinates system [START_REF]Methods of theoritical physics[END_REF][START_REF] Lebedev | Special Functions and their Applications[END_REF] or bipolar coordinates system (if we ignore the angular toroidal variable) (ζ, η) ∈ IR + * × [0, 2π] about the pole F p = (r p , z p ) is related to the cylindrical coordinates system (r, z) by

r = r p sinh ζ cosh ζ -cos η and z -z p = r p sin η cosh ζ -cos η
We assume that F p lies inside the plasma domain where the homogeneous equation, -∆ * ψ = 0, is not satisfied.

Following [START_REF] Faugeras | 2D interpolation and extrapolation of discrete magnetic measurements with toroidal harmonics for equilibrium reconstruction in a Tokamak[END_REF] the poloidal flux in the vacuum surrounding the plasma is written as ψ = ψ C + ψ th . The term ψ C is the given contribution of the poloidal field coils C k with current density j C k = I k S k evaluated thanks to the free space Green function:

ψ C (x) = k C k j C k G(y, x)dy
The toroidal harmonics expansion is represented by ψ th , sum of two terms ψ ex and ψ in , corresponding respectively to external and internal harmonics:

                               ψ th = ψ ex + ψ in , ψ ex = r p sinh ζ √ cosh ζ -cos η × [ n e a n=0 a e n Q 1 n-1/2 (cosh ζ) cos(nη) + n e b n=1 b e n Q 1 n-1/2 (cosh ζ) sin(nη)], ψ in = r p sinh ζ √ cosh ζ -cos η × [ n i a n=0 a i n P 1 n-1/2 (cosh ζ) cos(nη) + n i b n=1 b i n P 1 n-1/2 (cosh ζ) sin(nη)], (22)
where P 1 n-1/2 and Q 1 n-1/2 are the associated Legendre functions of first and second kind, of degree one and half integer order [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF], also called toroidal harmonics when evaluated at point cosh ζ. Functions P 1 n-1/2 have a singularity when ζ → ∞ that is to say at point F p . On the contrary functions Q 1 n-1/2 are singular when ζ → 0 that is to say on the axis r = 0. We denote by u the vector of the unknown expansion coefficients ).

The optimal set of coefficients u opt can be computed from the minimization of the least-square cost function

J 0 (u) = 1 2 ||M u -m|| 2 ( 23
)
where m is the vector of (normalized) flux loops and Bprobes measurements corrected from their known ψ C contribution, and M u the equivalent quantities computed from the toroidal harmonics expansion of Eq. ( 22). Here the rectangular matrix M depends on the choice of the pole of the coordinates system. J 0 is quadratic in u and is minimized simply by solving the associated normal equation, M T M u -M T m = 0, to find the unique optimal set of coefficients u opt .

Assuming enough toroidal harmonics are used the resulting representation of the poloidal field is always precise on Γ the measurements contour but might become highly disturbed when moving towards the pole F p which is a singular point for the internal harmonics. Hence let us introduce the regularized functional

J(u) = J 0 (u) + R(u) (24)
where the regularization term R involves the second order tangential derivative of ψ along a particular circle C surrounding the pole F p ,

R(u) = ε 2 C | d 2 ψ th (u) ds 2 | 2 ds
and is quadratic in u.

The following remark leads to an obvious choice for the circle C. In toroidal coordinates the curves of constant ζ are non-intersecting circles of different radii (r -r p coth ζ) 2 + (z -z p ) 2 = r 2 p sinh 2 ζ that surround the pole but are not concentric. The ζ = 0 curve corresponds to the z-axis (r = 0). As the magnitude of ζ increases, the radius of the circles decreases and their centers approach the pole F p . We choose the circle C to be the circle of a constant ζ 0 value. This choice considerably simplifies the computation of the second order tangential derivative since it is only required to compute derivatives with respect to η to obtain the expression of the regularization term.

The curvilinear abscissa along the constant-ζ 0 circle is given by

s(η) = η -π r p cosh ζ 0 -cos t dt
The second order tangential derivative along the circle of ψ seen as a function of s is

d 2 ψ th ds 2 = ∂ 2 ψ th ∂η 2 (η, ζ 0 )( dη ds ) 2 + ∂ψ th ∂η (η, ζ 0 ) d 2 η ds 2 with dη ds = cosh ζ 0 -cos η r p and d 2 η ds 2 = (cosh ζ 0 -cos η) sin η r 2 p
and we obtain the expression of the regularization term

C ζ 0 | d 2 ψ th ds 2 | 2 ds = π -π [ ∂ 2 ψ th ∂η 2 (η, ζ 0 )( cosh ζ 0 -cos η r p ) 2 + ∂ψ th ∂η (η, ζ 0 ) (cosh ζ 0 -cos η) sin η r 2 p ] 2 [ r p cosh ζ 0 -cos η ]dη (25)
The ψ th derivatives with respect to η can be easily explicitly computed from Eq. ( 22) and the integral computed using a standard quadrature method.

Let us conclude this section with the complete algorithm implemented in NICE.

Algorithm.

1. Initialization: prepare filaments description of PF coils, provide code parameters for number of toroidal harmonics, regularization parameters, test parameter d (see point 4 below), points defining contour Γ, number N p of required plasma boundary points 2. Choose the pole of the coordinates system F p . Default initial guess is F p = (r 0 , 0). Assemble and minimize cost function [START_REF] Hinze | Mathematical Modelling: Theory and Applications[END_REF] to find u opt . 3. Compute the current center F c = (r c , z c ) defined as moments of the plasma current density [START_REF] Zakharov | Equilibrium of a toroidal plasma with noncircular cross-section[END_REF][START_REF] Braams | The interpretation of tokamak magnetic diagnostics[END_REF] which can be precisely computed as integrals on the limiter contour Γ L at every point of which the flux and the poloidal magnetic field can be evaluated thanks to the representation ψ = ψ th (u opt ) + ψ C ,

I p = Γ L 1 µ 0 B s ds (26)
z c I p = Γ L 1 µ 0 (-r log rB n + zB s )ds (27)
r 2 c I p = Γ L 1 µ 0 (2rzB n + r 2 B s )ds ( 28
)
where B s (resp. B n ) denotes the poloidal field tangent (resp. normal) to the integration contour Γ L .

Numerical examples

This mode of NICE also known as VacTH is routinely used at WEST for plasma shape identification and control during a shot [START_REF] Moreau | The new magnetic diagnostics in the west tokamak[END_REF]. Here is provided a typical example of plasma boundary reconstruction for a WEST-like equilibrium. The poloidal flux is computed using 24 toroidal harmonics for ψ th (15 external harmonics with n e a = n e b = 7 and 9 internal harmonics with n i a = n i b = 4) and 24 filaments of currents modelizing the given current density in the PF coils for ψ C . Figure 5 shows the fit to Bprobes measurements with a resulting root mean square error of 2mT.

The ζ value defining the regularization circle is chosen such that the radius of this circle is a(r limiterM ax -r limiterM in ) where a = 0.1 is a code parameter. Another code parameter provides the value of regularization parameter which is here set to ε = 0.001. A 60 points plasma boundary is computed by the method described above and is shown on Figure 6.

The algorithm is not machine dependent and was successfully applied to other tokamaks such as JET, TCV, AUG or also ISTTOK [START_REF] Santos | Plasma boundary reconstruction in ISTTOK using magnetic diagnostic data[END_REF]. Figure 7 shows a JET-like equilibrium example and Figure 8 a TCV-like equilibrium example.

NICE equilibrium reconstruction mode

The goal of this section is to introduce the algorithm implemented in NICE for equilibrium reconstruction. Basic measurements used for equilibrium reconstruction are magnetic measurements but NICE can also use interferometry, polarimetry, motional Stark effect (MSE) and pressure measurements. These are introduced in Section 4.2 and we start by using magnetics only for simplicity.

Problem formulation on a bounded domain with magnetic measurements

Here magnetic measurements are in fact Cauchy data (ψ, ∂ψ ∂n ) provided on the measurement contour Γ. These data are usually computed from a first call to NICE magnetics interpolation mode described in Section 3 but can also be provided by another code specific to a Tokamak which computes the poloidal flux outside the plasma such as XLOC at JET [START_REF] O'brien | Local expansion method for fast plasma boundary identification in JET[END_REF][START_REF] Sartori | JET real-time object-oriented code for plasma boundary reconstruction[END_REF].

Using the matrix notations of Sections 2.2, 2.3 and 2.4 Dirichlet data provide the ψ d vector used in the resolution of the boundary value problem whereas Neumann data are considered as measurements and used in a cost function to be minimized to reconstruct the equilibrium. NICE uses N b values of ∂ψ ∂n given at the middle point of each segment of the mesh forming the Γ contour.

We denote by m the vector of these (normalized) measurements and the equivalent quantities, Hψ + H d ψ d computed from the FEM representation The cost function, which depends on y := ψ and u := (u A , u B ) the degrees of freedom of functions A and B, we consider is the following

J(y, u) := 1 2 ||Hψ + H d ψ d -m|| 2 + 1 2 ||R A u A || 2 + 1 2 ||R B u B || 2 (29)
The first term on the right hand-side is the misfit term and the two others are regularization terms. They are the discrete expressions for

ε A 2 1 0 (A (x)) 2 dx + α A 2 |A(1)| 2 and ε B 2 1 0 (B (x)) 2 dx + α B 2 |B(1)| 2
Parameters ε enable to tune the smoothness of the identified functions whereas parameters α tune the penalization to zero of their value on the plasma boundary.

The free-boundary equilibrium model equation is e(y, u)

:= Aψ + A d ψ d -J(ψ, λ, u A , u B ) -Lu I = 0 ( 30
)
and the equilibrium reconstruction problem is formulated as the folowing PDE-constrained optimization problem: min y,u J(y, u) such that e(y, u) = 0 [START_REF] O'brien | Local expansion method for fast plasma boundary identification in JET[END_REF] This problem is solved thanks to the QSQP algorithm given in Section 2.3. A first initialization step is used to set λ in the current density term given the total plasma current I p which is either given with magnetic measurements or comes out from the toroidal harmonics procedure. Given initial guess ψ 0 , u 0

A and u 0 B , λ is chosen to satisfy

I p -λJ p (ψ 0 , u 0 A , u 0 B ) = 0 (32)
Using the notation of Section 2.3, the ingredients needed for the QSQP algorithm in equations ( 15), ( 16), ( 17), ( 18), [START_REF] Blum | Plasma equilibrium evolution at the resistive diffusion timescale[END_REF] and [START_REF] Albanese | On the solution of the magnetic flux equation in an infinite domain[END_REF]

are ∆y = -(A -J ψ (ψ n , λ, u n A , u n B )) -1 (Aψ n + A d ψ d -J(ψ n , λ, u n A , u n B ) -Lu I ). ( 33
) S = -(A -J ψ (ψ n , λ, u n A , u n B )) -1 (J u (ψ n , λ, u n A , u n B )). ( 34
)
J T y = H T (Hψ + H d ψ d -m), J T u = R T A R A u A R T B R B u B (35)
J yy = H T H, J uu = R T A R A 0 0 R T B R B
, J yu = 0, J uy = 0 (36)

Internal measurements

In addition to magnetics NICE can use internal measurements for equilibrium reconstruction.

Interferometry and polarimetry

Classical polarimetry. Polarimetry measurements in first approximation give the Faraday rotation of the angle of an infrared radiation crossing the section of the plasma along different lines of sight, L i :

α i obs ≈ L i cN e B p .tdl
where c is a constant, N e is the electron density in the plasma and t is the unit vector tangent to the line of sight. In order to be able to use polarimetric measurements the electron density N e has to be known. It is therefore also being identified thanks to interferometric measurements which give the density line integrals over each of the N L line of sights

N i e,obs ≈ L i

N edl

In NICE the electron density, assumed constant on magnetic surfaces, is represented as

N e (ψ N ) = λ N e C(ψ N , u C )
where λ N e is a scaling parameter and function C is either decomposed in a basis of piecewise linear or cubic spline functions

C(x, u C ) = N C i=1 u C,i φ i (x)
or given as a parametric function [START_REF] Groebner | Progress in quantifying the edge physics of the H mode regime in DIII-D[END_REF] C

(x, u C ) = (M -m)(1 -αhz)E(z) + m
where E(z) = e -z e -z + e z , z =

x -x 0 h and the vector of degrees of freedom is

u C = (x 0 , h, α, m, M ).
The interferometry observation operator W(ψ, u C ) is a vector of size N L computed from the FEM discretiezation as

W(ψ, u C ) i = L i λ Ne C(ψ N (ψ h ), u C )dl
and for polarimetry we define Z(ψ, u C ) in a similar way by

Z(ψ, u C ) i = L i cλ Ne C(ψ N (ψ h ), u C )B p (ψ h ).tdl
Here again we have omitted normalization by the assumed observation error. Integrals along the lines of sight are approximated using the trapezoidal quadrature rule.

The vector of state variable is still y := ψ whereas the one of all control variables is now u = (u A , u B , u C ) and we can now define the cost function

J(y, u) = 1 2 ||Hψ + H d ψ d -m|| 2 + 1 2 ||W(ψ, u C ) -m N e || 2 + 1 2 ||Z(ψ, u C ) -m F ar || 2 + 1 2 ||R A u A || 2 + 1 2 ||R B u B || 2 + 1 2 ||R C u C || 2 (37)
which has to be minimized under the constraint of the equilibrium model [START_REF] Santos | Plasma boundary reconstruction in ISTTOK using magnetic diagnostic data[END_REF].

In order to use the SQP algorithm we need to compute derivatives

W ψ (ψ, u C ) W u C (ψ, u C ) Z ψ (ψ, u C ) and Z u C (ψ, u C ).
In NICE we compute the exact derivatives of the discrete operators. This is done on the fly in the same loop as the one used for the assembly of vectors W(ψ, u C ) and Z(ψ, u C ).

The computation of ∆y and S in ( 33) and ( 34) are unchanged. Concerning the gradients of the cost function we have

J T y (y, u) = H T (Hψ + H d ψ d -m) + W ψ (ψ, u C ) T (W(ψ, u C ) -m Ne ) + Z ψ (ψ, u C ) T (Z(ψ, u C ) -m F ar ) ( 38
)
J T u =   R T A R A u A R T B R B u B R T C R C u C + W T u C (W -m Ne ) + Z T u C (Z -m F ar )   ( 39
)
Concerning the second order derivative terms we avoid a complete computation by only considering terms involving products of first order derivatives and form the following quantities

J yy = H T H + W T ψ W ψ + Z T ψ Z ψ ( 40
)
J uu =   R T A R A 0 0 0 R T B R B 0 0 0 R T C R C + W T u C W u C + Z T u C Z u C   ( 41
)
J yu = 0 0 W T ψ W u C + Z T ψ Z u C , J uy =   0 0 W T u C W ψ + Z T u C Z ψ   (42)
Stokes model for polarimetry. The project of performing equilibrium reconstruction using interferometry and Stokes model for polarimetry is at the origin of the development of NICE. It was first presented in [START_REF] Faugeras | Optimal control of a coupled partial and ordinary differential equations system for the assimilation of polarimetry stokes vector measurements in tokamak free-boundary equilibrium reconstruction with application to ITER[END_REF] where a Newton-Lagrange SQP algorithm is used in the ABB domain case for ITER.

In this section we briefly restate the discretization method implemented for Stokes model and describe the more recent QSQP method used for the inverse equilibrium reconstruction problem. The method presented here was used in [START_REF] Faugeras | Equilibrium reconstruction at JET using Stokes model for polarimetry[END_REF] with real measurements at JET. Stokes model consists in a system of differential equations for the the Stokes vector s of dimension 3 to be solved along each line of sight. The system is linear in s but coefficients depend on the electron density and on the magnetic field. Each of the N L line of sight L i is discretized in N i points and a Crank-Nicolson scheme is implemented for the integration of the Stokes model. This can be written

M i (ψ, u B , u C )S i -S i 0 = 0
where M i is a 3N i × 3N i band diagonal matrix, S i = (s i,1 , . . . , s i,N i ) is the vector of all Stokes vector states along the line of sight and S i 0 = (s i 0 , 0, . . . , 0) represents the initial conditions.

The measured quantity can be assumed to be the Stokes vector at the end of the line of sight after integration through the plasma and the observation operator is given a matrix E i of size 3 × 3N i such that E i S i = s i,N i .

Defining the block diagonal matrices

M = diag(M 1 , . . . , M N L ) and E = diag(E 1 , . . . , E N L ),
and collecting the Stokes vector states for all line of sights in a vector

S = (S 1 , . . . , S N L )
one can define a direct model e(y, u) = 0 for the state variable y = (ψ, S)

and control parameters u = (u A , u B , u C ) by e(y, u) := Aψ + A d ψ d -J(ψ, λ; u A , u B ) -Lu I M (ψ, u B , u C )S -S 0 = 0 (43)
as well as the following cost function for equilibrium reconstruction using magnetics, interferometry and Stokes model polarimetry,

J(y, u) = 1 2 ||Hψ + H d ψ d -m|| 2 + 1 2 ||W(ψ, u C ) -m N e || 2 + 1 2 ||ES -m s || 2 + 1 2 ||R A u A || 2 + 1 2 ||R B u B || 2 + 1 2 ||R C u C || 2 (44)
which is minimized under the constraint of the model (equation 43). The Jacobians needed in Eqn. ( 17) and ( 18) of the QSQP algorithm are

e y (y, u) = A -J ψ (ψ, λ; u A , u B ) 0 D ψ [M (ψ, u B , u C )S] M (ψ, u B , u C ) (45)
and

e u (y, u) = -J u A (ψ, λ; u A , u B ) -J u B (ψ, λ; u A , u B ) 0 0 D u B [M (ψ, u B , u C )S] D u C [M (ψ, u B , u C )S] (46
) Concerning the gradients of the cost function we have

J T y (y, u) = H T (Hψ + H d ψ d -m) + W ψ (ψ, u C ) T (W(ψ, u C ) -m Ne ) E T (ES -m s ) (47)
J T u =   R T A R A u A R T B R B u B R T C R C u C + W T u C (W -m Ne )   (48)
Concerning the second order derivative terms we avoid a complete computation by only considering terms involving products of first order derivatives and form the following quantities

J yy = H T H + W T ψ W ψ 0 0 E T E (49)
J uu =   R T A R A 0 0 0 R T B R B 0 0 0 R T C R C + W T u C W u C   (50)
J yu = 0 0 W T ψ W u C 0 0 0 , J uy =   0 0 0 0 W T u C W ψ 0   (51)
Finally let us mention that in [START_REF] Faugeras | On the identification of the electron temperature profile from polarimetry Stokes vector measurements in tokamak free-boundary equilibrium reconstruction[END_REF] dependence on electron temperature is considered in the Stokes model and the numerical experiments are performed for the identification of the T e profile together with p , f f and the N e profiles. As for other profiles the T e profile is decomposed in a function basis with control parameters u D . With the notations of this Section this mainly consists in augmenting the control vector to u = (u A , u B , u C , u D ) and adding a dependence of the Stokes model to T e that is to say using a matrix M (ψ, u B , u C , u D ).

Motional Start effect

Motional Stark effect (MSE) measurements are angular measurements depending on the magnetic field inside the plasma domain. They are given at N M SE points p i = (r i , z i ) as

m i M SE = tan γ i = a i 0 B z (p i ) + a i 1 B r (p i ) + a i 2 B φ (p i ) a i 3 B z (p i ) + a i 4 B r (p i ) + a i 5 B φ (p i )
where the a coefficients are given constants. In NICE the MSE observation operator M(ψ, u B ) is the vector of components

M(ψ, u B ) i = a i 0 B z (ψ h (p i )) + a i 1 B r (ψ h (p i )) + a i 2 S f (ψ N (ψ h (p i )))/r i a i 3 B z (ψ h (p i )) + a i 4 B r (ψ h (p i )) + a i 5 S f (ψ N (ψ h (p i )))/r i J yu = P T ψ P u A 0 , J uy = P T u A P ψ 0 (64) 4.3.
Outputs and uncertainty quantification 4.3.1. Equilibrium outputs, flux surface averages and geometric coefficients Numerous outputs can be extracted from the equilibrium poloidal flux map computed. These include purely geometric information on the plasma shape (plasma boundary, geometric axis, elongation . . . ), global parameters (such as total plasma current I p , poloidal beta β p , internal inductance li, . . . ), 1D profiles of quantities constant on flux isolines in the plasma and 2D maps (ψ itself but also B r , B z , j p , . . . ). All these outputs are standardized and follow the IMAS conventions. We are not going to detail all of them in this paper. Let us however give some details on the computation of some of the important 1D profiles in the plasma.

For

ψ N ∈ [0, 1], S f (ψ N ) = f (ψ)
is computed by integration of B as shown in (52). Let us define a discretization of the unit interval [0, 1] by S +1 values ψ 0 N = 0, . . . , ψ S N = 1. These points are taken as abscissa for all computed 1D profiles. For each ψ s N the contour line Γ ψ s N is extracted from the finite element representation of the solution as a list of N s segments between m l s,1 = (r l s,1 , z l s,1 ) and m l s,2 = (r l s,2 , z l s,2 ) with length |L l s |, for l = 1 to N s . The toroidal flux coordinate is defined as ρ(ψ N ) = φ(ψ N )/πB 0 where φ(ψ

N ) = Ω ψ N f (ψ h (r,z)) r
drdz and Ω ψ N is the domain bounded by the line of flux Γ ψ N . This contour is not resolved by the mesh and the quantities φ s and ρ s are computed from the discrete ψ h for all ψ s N using the same barycentric quadrature rule as for the integrated plasma current density J in Section 2.2:

φ s = T S f (ψ N (b T (ψ h ))) r T (ψ h ) |T ∩ Ω ψ s N |. (65)
The profiles ψ s and ρ s being known one can compute ( ∂ψ ∂ρ ) s = ψ s using finite differences.

In the same way the volume profile is computed as

V ol s = 2π T r T (ψ h )|T ∩ Ω ψ s N | (66)
and ( ∂V ol ∂ρ ) s = V ol s using finite differences.

Following [START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks[END_REF] the average of a quantity A over magnetic surfaces can be computed as

A s = ( Γ ψ s N A r |∇ψ h | dl)/( Γ ψ s N r |∇ψ h | dl). (67)
A number of 1D profiles, also called geometric coefficients, are computed as such averages, e.g. 1/r 2 or |∇ρ| 2 /r 2 . The integrals over flux contour lines involved are approximated as follows:

Γ ψ s N A r |∇ψ h | dl ≈ Ns l=1 1 2 r l s,1 A(m l s,1 ) |∇ψ h|T l s | + r l s,2 A(m l s,2 ) |∇ψ h|T l s | |L l s |. ( 68
)
where T l s is the triangle which is intersected by the segment between m l s,1 and m l s,2 and m l s,• = (r l s,• , z l s,• ) . |∇ψ h|T l s | is constant in the triangle and computed from the 3 values at the nodes of T l s .

Sensitivity method for error bars computation

At convergence of the SQP iterations an optimal u is found and the constraint given by the model is satisfied. Hence y = y(u) and one can define the reduced cost function Ĵ(u) = J(y(u), u). A good approximation to the Hessian of the reduced cost function is SQP algorithm matrix, Ĵuu ≈ M . This matrix would be the exact Hessian if we hadn't neglected some of the second order derivatives terms.

Its diagonal elements are the inverse of the squared a posteriori standard deviations on the components of u also called error bars. Moreover the error bar of any scalar quantity g depending on u (such as the reconstructed functions A(ψ N ) or B(ψ N ) at a given ψ N for example) is computed as

σ 2 g = ∇ u g T M -1 ∇ u g
This enables the computation of the error bar on any output quantity at low numerical cost. If g(u) is in the form g(u) = G(y(u), u) we use

∇ u g = y T u ∇ y G + ∇ u G
where the Jacobian y u = S is already known from SQP algorithm.

Final remarks on equilibrium reconstruction

The equilibrium reconstruction problem can be formulated on the ABB domain. With this formulation magnetic Bprobes and flux loops measurements are directly used without having to perform a first interpolation step on the measurement contour as with the formulation on a bounded domain. However there are two drawbacks to this method. First of all it is computationally more demanding since the mesh is necessarily much larger, and secondly it does not work well for tokamaks with ferromagnetic structures since the modelization of the magnetic permeability is certainly not very precise. Nevertheless the method works well for iron-free tokamaks and a numerical example is provided in Section 5.3. This method is closer to the algorithms implemented in codes like EFIT [START_REF] Lao | Equilibrium analysis of current profiles in Tokamaks[END_REF], EQUAL [START_REF] Zwingmann | Equilibrium analysis of steady state tokamak discharges[END_REF], CLISTE [START_REF] Mc Carthy | The CLISTE Interpretive Equilibrium Code[END_REF] or LIUQE [START_REF] Moret | Tokamak equilibrium reconstruction code LIUQE and its real time implementation[END_REF], which however do not make use of derivative with respect to the state variable as is done here in the QSQP algorithm.

The cost function to be minimized is very similar to ( 29)

J(ψ, u) = 1 2 ||H B ψ -m B || 2 + 1 2 ||H F ψ -m F || 2 + 1 2 ||R A u A || 2 + 1 2 ||R B u B || 2 (69)
Here m B (resp. m F ) are the Bprobes measurements (resp. the flux loops measurements) and H B (resp. H F ) the associated linear observation operators mapping the FEM representation of the poloidal flux to the measurements.

For an iron-free tokamak the model equation constraint reads

(A + C)ψ -J(ψ, λ, u A , u B ) -Lu I = 0 (70)
The QSQP algorithm follows as in the preceding sections.

Clearly the different measurements of the previous sections (magnetics, interferometry, polarimetry, MSE, pressure) can be used all together. In NICE this is controled by code parameters and the user can decide to use any combination of measurements. Finally one can also perform a magnetics only reconstruction and then use the obtained flux map to reconstruct in a second step the electron density profile using interferometry measurements only. rapidly converges to the tolerance value set to 10 -10 . Figure 10 shows the identified p and f f profiles as well as the obtained averaged current density and safety factor. Computed error bars on these profiles are also shown.

WEST-like equilibrium reconstruction with magnetics and interfero-polarimetry. This second example shows results for a WEST-like equilibrium reconstrucion using magnetics and interfero-polarimetry measurements. Figure 11 shows . Figure 12 shows the identified p and f f profiles as well as the obtained averaged current density and safety factor. Computed error bars on these profiles are also shown. Figure 13 shows the identified electron density profile with error bars as well as the interferometry and polarimetry data. Because of invalid measurements 8 chords out of 10 are used for interferometry while 7 are used for polarimetry. JET-like equilibrium reconstruction with magnetics and pressure. This third example shows results for a JET-like equilibrium reconstrucion using magnetics and pressure measurements. Pressure measurements are calculated from High Resolution Thomson Scattering diagnostic with the assumption n i = n e and T i = T e , and data mapped onto flux surface coordinates. Figure 14 shows the computation domain and the obtained flux map. Cauchy measurements on the domain boundary are computed with the first mode of NICE described in Section 3 using external harmonics of order 5 and inner of order 3, then the QSQP algorithm of section 4.2.3 is run. Function A is decomposed in a basis of 14 cubic splines and the regularization parameter ε A varies radially from 10 -1 on the magnetic axis to 10 -3 near the plasma boundary where the function is left free. Function B is decomposed in a basis of 12 cubic splines and the regularization parameter is ε A = 10 -1 and B( 1) is forced to 0 with a large penalization parameter. The relative residual rapidly converges to the tolerance value set to 10 -10 . Figure 15 shows the identified p and f f profiles as well as the obtained averaged current density and safety factor. Computed error bars on these profiles are also shown.

Figure 16 shows the identified pressure profile with error bars as well as the pressure data.

NICE direct and inverse equilibrium computation modes

Static equilibrium computation

Two types of direct static equilibrium computation can be performed with NICE, I p -free (in which the total plasma current is not imposed) or I p -fixed computations (in which the total plasma current is imposed to a given value). Although it is possible to perform this type of computation in the bounded domain case it is most usually performed in the ABB domain case where all poloidal field coils are included. In both types of computation the currents in these coils are given, that is u := u I the control variable for the inverse problem to come, is a given fixed parameter.

In I p -free computations the current density functions p and f f are given, that is λ, u A and u B are given and fixed. The state variable is y := ψ and the model equation is e(y, u) :

= (A + C)ψ + A µ (ψ) -J(ψ, λ; u A , u B ) -Lu I = 0. (71)
In I p -fixed computations the current density functions are given up to the Similar to Section 5.1, in I p -free computations the current density functions p and f f are given that is λ k , u k A and u k B are given at each time step. In I p -fixed computations the current density functions are given up to the scaling factor λ which is computed at each time given a prescribed value of the plasma current I k p . The voltage v i (t) of each of the N v power supplies is a given function of time decomposed in a basis of N bv piecewise linear or spline functions

v i (t) = N bv j=1 u i,j φ j (t)
Noting u V the vector of size N v N bv of all decomposition coefficients u i,j for all power supllies, the

N v × N v N bv matrix B(t) is such that B(t)u V = v(t)
the vector of all voltages.

As for the static case inverse problems can be defined for the quasi-static evolution case. It consists in finding the voltages, hence controling u = u V , such that the plasma shape follows a prescribed evolution. In order to stick with the notations of this document one can define a state variable concatenating the states at each time step y := (y 1 , . . . y N T ) and a model e(y, u) :

=   e 1 (y 1 , y 0 , u V ) . . . e N T (y N T , y N T -1 , u V )  
The simplest cost function which can be defined is the following

J(y, u) = 1 2 N T k=1 N d i=1 w k (ψ k h (x k i ) -ψ k h (x k 0 )) 2 + 1 2 w N T k=1 Nv i=1 ([B k u V ] i ) 2
with a misfit term to a desired plasma boundary at each time step and a penalization term on voltages. This can be rewritten

J(y, u) = 1 2 ||Kψ|| 2 + 1 2 ||Ru V || 2
and the QSQP algorithm to solve the PDE-constrained optimization problem follows. Other types of penalization terms are also implemented in NICE: penalization on the maximum and minimum voltage values and penalization on induced currents in passive strutures. we refer to [START_REF] Blum | Automating the design of tokamak experiment scenarios[END_REF] for a study of this type of inverse problems for the design of tokamak scenarios.

Numerical examples

JT60-SA equilibrium reconstruction in ABB domain with synthetic magnetics generated from a static inverse. In this first example we perform an equilibrium reconstruction with synthetic magnetic data for JT60-SA tokamak in ABB domain. In a first step, the Ip-fixed static inverse mode of NICE is run to compute the currents in the PF coils giving a desired plasma boundary shape provided as a list of input points. The 12 reference currents in (73) are u I 0 = 0 and the penalization weights are chosen as σ i = 10 8 . A parametric representation for A and B is given. The relative residual rapidly converges to 10 -10 . This provides a reference equilibrium with known p and f f functions. From the computed flux map we compute the equivalent of magnetics measurements. A 1% noise is added to these measurements. Then in a second step we run the equilibrium reconstruction mode of NICE in ABB domain that is to say without using the toroidal harmonics step. Functions A and B are to be identified in a basis of 11 cubic splines each. The regularization parameters are ε A = ε B = 10 -2 . The computed poloidal flux map is shown on Figure 17. Figure 18 shows the reference and identified p and f f profiles as well as the averaged current density and safety factor. Computed error bars on these profiles are also shown.

ITER vertical displacement event. In this second example an initial ITER plasma is computed thanks to the Ip-fixed static inverse mode of NICE in a first step. It is shown on Figure 19. The desired plasma shape is provided thanks to a parametric representation.

Then in a second step the Ip-fixed quasi-static evolution mode is run. All voltages are set to zero and the time step is ∆t = 10 -3 s. At each time step the Newton relative residual converges in few iterations to the tolerance value of 10 -10 . A vertical displacement event occurs after approximately 1s. The evolution of the plasma is shown on Figure 20 and the final poloidal flux map at t = 1.5s is shown on Figure 21.

ITER quasi-static evolution inverse mode. In this third example we use the quasi-static evolution inverse mode of NICE to compute the voltages which drive an ITER limiter plasma upwards. The different desired plasma boundaries as time increases are shown on Figure 22 together with the initial poloidal flux map. The comutation is in Ip-fixed mode and function A and B are given by a parametric representation. Each of the 12 voltages time dependent functions to be identified are decomposed in a basis of 11 piecewise linear functions. The weights in the regularization term are chosen as σ i = 10 -5 , 11 time steps of 0.1s are considered. The size of the control vector u is 132. The mesh size is 12574 giving a state variable vector y of size 138314. The QSQP algorithm converges to the tolerance value of 10 -10 in 10 iterations. Figure 23 shows the poloidal flux map at final time with a good match between desired and computed plasma boundaries. Figure 24 shows the voltages computed by the QSQP algorithm.

Conclusion

This document provides a general overview of the numerical methods implemented in the different effective modes of the code NICE. The code is available on the svn repository of the EUROFUSION gateway. It is still evolving and new features a regularly added. Currently the possibility of using higher order C 1 finite elements is under development. Other future developments will also deal with the use of the 1D resistive diffusion equation either in a coupling with the quasi-static evolution mode or as an additional constraint in successive equilibrium reconstructions.

Aknowledgments

The author would like to thank Cédric Boulbe, Hervé Guillard and two anonymous reviewers for their constructive comments on the text, as well as Holger Heumann and Jacques Blum in particular for their inputs on SQP methods.

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Figure 1 :

 1 Figure1: Schematic representation of the poloidal plane of a tokamak. Ferromagnetic iron structures (Ω f ) are represented in gray, PF coils (Ω ci ) in orange and a passive structure (Ω ps ) in blue. The limiter contour (Γ L , purple contour) inside which the plasma (Ω p ) lies is also shown. Cross-circles represent magnetic probes and a fictitious measurement contour Γ is drawn with a green dashed line.

Figure 2 :

 2 Figure 2: Example mesh for WEST tokamak in the bounded domain case. In red is the region inside the limiter.

Figure 3 :

 3 Figure 3: The ABB domain for WEST tokamak.

Figure 4 :

 4 Figure 4: Zoom on an example mesh in the ABB domain case for WEST tokamak. Different parts of the machine appear in different colors (iron structure, PF coils, vacuum vessel, limiter).

4 .

 4 If ||F c -F p || > d set F p := F c , re-assemble and minimize cost function (24) again, since matrix M in (23) depends on the pole of the coordinates system, to find a new u opt . 5. Output Cauchy conditions (ψ, ∂ n ψ) on measurement contour Γ and/or compute plasma boundary with the following algorithm: find P b the point defining the plasma boundary ψ b value. Either a limiter point or an X-point found thanks to a quasi Newton method for ∇ψ = 0 starting from the [F c , P b ) ray define N p -1 rays [F c , P k ) with a rotation of ∆θ = 2π/(N p -1) between each of them. Use Newton method along each ray to find the plasma boundary point where ψ(r, z) -ψ b = 0

Figure 5 :

 5 Figure 5: Measured and computed Bprobes values for a WEST-like plasma boundary reconstruction.

Figure 6 :

 6 Figure 6: Poloidal flux ψ = ψ th (u opt ) + ψ C computed by a fit of toroidal harmonics to magnetic measurements for a WEST-like equilibrium. The measurement contour Γ is shown in black, Bprobes are black circles and flux loops black squares. The limiter contour is the dashed black line. The plasma boundary with a low Xpoint is the thick red contour whereas the regularization circle is the thin red line surrounding the pole of the coordinate system.

Figure 7 :

 7 Figure 7: Poloidal flux ψ = ψ th (u opt ) + ψ C computed by a fit of toroidal harmonics to magnetic measurements for a JET-like equilibrium. The measurement contour Γ is shown in black, Bprobes are black circles and flux loops black squares. The limiter contour is the dashed black line. The plasma boundary with a low Xpoint is the thick red contour whereas the regularization circle is the thin red line surrounding the pole of the coordinate system.

Figure 8 :

 8 Figure 8: Poloidal flux ψ = ψ th (u opt ) + ψ C computed by a fit of toroidal harmonics to magnetic measurements for a TCV-like equilibrium. The measurement contour Γ is shown in black, Bprobes are black circles and flux loops black squares. The limiter contour is the dashed black line. The plasma boundary with an upper Xpoint is the thick red contour whereas the regularization circle is the thin red line surrounding the pole of the coordinate system.

Figure 9 :

 9 Figure 9: Poloidal flux map for a TCV-like equilibrium reconstruction with magnetic measurements. The plasma boundary is shown in red. The limiter contour is the doted line.

4. 5 .

 5 Numerical examples TCV-like equilibrium reconstruction with magnetics only. This first example shows results from a TCV-like equilibrium reconstrucion using magnetics only. Figure9shows the computation domain and the obtained flux map. Cauchy measurements on the domain boundary are computed with the first mode of NICE described in Section 3 using toroidal harmonics of order 5, then the QSQP algorithm of section 4.1 is run. Functions A and B are decomposed in a basis of 6 cubic splines each. The regularization parameters are ε A = ε B = 0.1. The relative residual, relresid = ||(y, u) n+1 -(y, u) n || ||(y, u) n || ,

Figure 10 :

 10 Figure 10: TCV-like equilibrium reconstruction with magnetics. Reconstructed p , f f , j tor = (p + < 1/r 2 > (1/µ 0 )f f )/ < 1/r > and safety factor q profiles, and associated computed error bars.

Figure 11 :

 11 Figure 11: Poloidal flux map for a WEST-like equilibrium reconstruction with magnetics and interfero-polarimetry. The plasma boundary is shown in red. The limiter contour is the doted line and the dashed lines are the interfero-polarimetry chords

Figure 12 :

 12 Figure12: WEST-like equilibrium reconstruction with magnetics and interferopolarimetry. Reconstructed p , f f , j tor = (p + < 1/r 2 > (1/µ 0 )f f )/ < 1/r > and safety factor q profiles, and associated computed error bars.

Figure 13 :

 13 Figure 13: WEST-like equilibrium reconstruction with magnetics and interferopolarimetry. Reconstructed electron density N e and associated computed error bars. Computed and measured interferometry and Faraday angle data.

Figure 14 :

 14 Figure 14: Poloidal flux map for a JET-like equilibrium reconstruction with magnetics and pressure measurements. The plasma boundary is shown in red. The limiter contour is the doted line.

Figure 15 :

 15 Figure15: JET-like equilibrium reconstruction with magnetics and pressure. Reconstructed p , f f , j tor = (p + < 1/r 2 > (1/µ 0 )f f )/ < 1/r > and safety factor q profiles, and associated computed error bars.

Figure 16 :

 16 Figure 16: JET-like equilibrium reconstruction with magnetics and pressure. Reconstructed pressure p and associated computed error bars. Computed and measured pressure data.

Figure 17 :

 17 Figure 17: JT60-SA equilibrium reconstruction with synthetic magnetics. Poloidal flux map in ABB domain.

Figure 18 :

 18 Figure18: JT60-SA equilibrium reconstruction with synthetic magnetics. Reconstructed p , f f , j tor = (p + < 1/r 2 > (1/µ 0 )f f )/ < 1/r > and safety factor q profiles, and associated computed error bars. The reference equilibrium are in black and the reconstructed ones in blue.

Figure 19 :

 19 Figure 19: Initial poloidal flux map in ITER geometry.

Figure 20 :

 20 Figure 20: Evolution of the plasma boundary during a VDE in ITER geometry for t = 1, 1.1, 1.2, 1.3, 1.4 and 1.5s

Figure 21 :

 21 Figure 21: Final poloidal flux map in ITER geometry at t = 1.5s.

Figure 22 :

 22 Figure 22: Initial poloidal flux map in ITER geometry and successive desired plasma boundaries for the time evolution.

Figure 23 :

 23 Figure 23: Final poloidal flux map in ITER geometry. Matching computed and desired plasma boundary are shown.

Figure 24 :

 24 Figure 24: Computed voltages driving the desired plasma boundary evolution.

  r 3 (∇λ i ∇ψ h )(∇λ j ∇ψ h )drdz boundary conditions. Dirichlet boundary conditions, a vector ψ d of size N b number of vertices on the boundary, are applied using matrix A d of size N × N b with coefficients

where we have noted S f (ψ N ) = f (ψ) which depends on u B since it is obtained from integration of B as follows

B 0 is the vacuum toroidal field at r = r 0 . During assembling of vector M(ψ, u B ) first order derivatives, M ψ (ψ, u B ) and M u B (ψ, u B ) are also computed. We can now formulate the cost function used for equilibrium reconstruction from magnetic and MSE measurement

which has to be minimized under the constraint of the equilibrium model [START_REF] Santos | Plasma boundary reconstruction in ISTTOK using magnetic diagnostic data[END_REF]. The computation of ∆y and S in ( 33) and ( 34) are unchanged. Concerning the gradients of the cost function we have

Very similar to what is done using interfero-polarimetry the following approximated second order derivative quantities are considered to form matrix M and vector m of the SQP algorithm

Pressure

The way pressure measurements can be obtained combining different diagnostics is out of the scope of this paper but when available they can be used for equilibrium reconstruction in NICE. Pressure measurements can be given as function of space, that is to say N p values p i given at points p i = (r i , z i ). They can also be given as a function ψ N that is to say N p values p i given for N p values ψ i N . In either case it leads to the construction of a pressure observation operator P(ψ, u A ) with components

where ψ i N is either a fixed given quantity or computed as ψ N (ψ h (p i )) which leads to differences in the computation of derivatives P ψ (ψ, u A ) and P u A (ψ, u A ) done during the assembling. We can now formulate the cost function used for equilibrium reconstruction from magnetic and pressure measurements

which has to be minimized under the constraint of the equilibrium model [START_REF] Santos | Plasma boundary reconstruction in ISTTOK using magnetic diagnostic data[END_REF]. The computation of ∆y and S in [START_REF] Groebner | Progress in quantifying the edge physics of the H mode regime in DIII-D[END_REF] and [START_REF] Faugeras | Equilibrium reconstruction at JET using Stokes model for polarimetry[END_REF] are unchanged. Concerning the gradients of the cost function we have

Similar to what is done for interfero-polarimetry and MSE the following approximated second order derivative quantities are considered to form matrix M and vector m of the SQP algorithm

scaling factor λ. The state variable is augmented with λ, y := (ψ, λ) and the model equation augmented with an equation imposing a given I p value e(y, u) :=

The direct static I p -free or I p -fixed computation consists in finding y such that e(y, u) = 0 using Newton's method.

The associated inverse problem is to find the currents in the coils, hence u such that the plasma has a desired shape and position. A number N p +1 of points x i = (r i , z i ) defining the desired plasma boundary are provided. They are either given from another simulation or provided from the parametric representation taken from [START_REF] Grandgirard | Modélisation de l'équilibre d'un plasma de tokamak[END_REF]Chap. 4,p84]. A typical cost function is

where

is the misfit term imposing an isoflux line to go through the given desired points and the second term is a penalization term in which the diagonal matrix R holds weights R ii = 1/σ 2 i . As for the equilibrium reconstruction problem the inverse problem min y,u J(y, u) such that e(y, u) = 0 is solved thanks to the QSQP method. The needed ingredients are

for the I p -free case and

for the I p -fixed case.

Quasi-static equilibrium evolution computation

In order to compute the quasi-static equilibrium evolution on time interval [0, T ] an implicit Euler scheme with time step ∆t is used starting from an initial condition y 0 . The computation of the equilibrium at t k from the one at t k-1 is done by solving the following problem: For the the I p -fixed case we have

with y k := (ψ k ; λ k ).


Laetitia Virlouvet

Charlotte Henriet

Marie-Pierre Jacquemot

Sandrine Chaignon

Thierry Balliau

Marie Javelle

Christophe Sallaud

Sandrine Balzergue

Marie-Laure Magniette

Johann Joets

Exploring the impact of ZmMYB31 overexpression in maize under contrasted soil water regimes

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Exploring the impact of ZmMYB31 overexpression in maize under contrasted soil water regimes

Laëtitia Virlouvet a , Charlotte Henriet a , Marie-Pierre Jacquemot a , Sandrine Chaignon a , Thierry Balliau b , Marie Javelle c , Christophe Sallaud c , Sandrine Balzergue d , Marie-Laure Martin-Magniette d,e , Johann Joets b , François Tardieu f , Claude Welcker Because leaf growth is one of the first processes affected by water deficit, many physiological studies concentrated in short-term responses and associated mechanisms. They demonstrated the roles of cellular and metabolic processes such as changes in cell turgor, hydraulic conductance and cell wall plasticity. However, our understanding of how water deficit impacts the cell wall biosynthesis and dynamics is still fragmentary. Here, we report that ZmMYB31, a subgroup 4 R2R3-MYB transcription factor that acts as a repressor of the lignin pathway in Arabidopsis [1][2], is induced by water deficit in the leaf growing zone and colocates with QTLs for lignin content and growth responses to water deficit in maize. To assess the role of ZmMYB31 in maize upon water limitation, we generated maize transgenic lines overexpressing ZmMYB31. We showed, in preliminary experiments, that the leaf growing zone of T1 ZmMYB31 overexpression plants had a higher content of β-O-4-linked lignin units than that of wild-type sister plants grown in the greenhouse under well-watered conditions. A comparative analysis by quantitative RT-PCR revealed significant differences for six lignin biosynthetic genes in the leaf growing zone of T1 ZmMYB31 overexpression plants compared with that of wild-type plants, consistent with the observed changes in -O-4linked lignin unit content. To obtain additional clues with regard to ZmMYB31 function in maize under water-limited conditions, we are currently combining RNA-seq analysis and shotgun proteomics of the leaf growing zone of homozygous ZmMYB31 overexpression plants and their wild-type siblings grown in the PhenoArch platform under well-watered and moderate water deficit conditions. Integrative statistical approaches shall allow us to infer a ZmMYB31-regulated network and identify variables underpinning maize responses to water deficit.

Keywords: maize, water deficit, ZmMYB31, lignin biosynthesis, omics

  Univ Paris Sud, Univ Evry, Univ Paris-Diderot, Sorbonne Paris-Cité, Université Paris-Saclay, Rue de Noetzlin, F-91405 Orsay, France; e UMR MIA-Paris, AgroParisTech, INRA, Université Paris-Saclay, F-75005 Paris, France ; f INRA, UMR 0759, Laboratoire d'Ecophysiologie des Plantes sous Stress Environnementaux, F-34060 Monptellier Cedex 1, France sylvie.coursol@versailles.inra.fr Water deficit directly impacts the ability of plants to intercept and convert light into biomass.

	f , Michel Zivy b , Valérie
	Méchin a , Sylvie Coursol a
	a Institut Jean-Pierre Bourgin, INRA, AgroParisTech, CNRS, Université Paris-Saclay, RD10,
	F-78026 Versailles Cedex, France; b Génétique Quantitative et Evolution -Le Moulon, INRA,
	Univ Paris-Sud, CNRS, AgroPArisTech, Université Paris-Saclay, F-91190 Gif-Sur-Yvette,
	France ; c Biogemma, Centre de Recherche, F-63720 Chappes, France ; d Institute of Plant
	Sciences Paris-Saclay IPS2, CNRS, INRA,


Laetitia Virlouvet

Sandrine Chaignon

Valérie Méchin

Marie-Pierre Jacquemot

Llorenç Cabrera-Bosquet

Thierry Balliau

Charlotte Henriet

Sandrine Balzergue

Etienne Delannoy

Johann Joets

Llorenç

MYB31 triggers water stress memory responses in maize

à 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.

MYB31 TRIGGERS WATER STRESS MEMORY RESPONSES IN MAIZE

Laetitia Virlouvet 1 , Sandrine Chaignon 1 , Valérie Méchin 1 , Marie-Pierre Jacquemot 1 , Llorenç Cabrera-Bosquet 2 , Thierry Balliau 3 , Charlotte Henriet 1 , Sandrine Balzergue 4,5 , Etienne Delannoy 4,5 , Johann Joets 3 , Marie-Laure Martin-Magniette 4,5,6 , Claude Welcker 2 , Michel Zivy 3 , Christophe Sallaud 7 , François Tardieu 2 , Marie Javelle 7 , Sylvie Coursol Water deficit directly impacts the ability of plants to intercept and convert light into biomass. Because leaf growth is one of the first processes affected by water deficit, many physiological studies concentrated in short-term responses and associated mechanisms. They demonstrated the roles of cellular and metabolic processes such as changes in cell turgor, hydraulic conductance and cell wall plasticity. However, our understanding of how water deficit impacts cell wall biosynthesis is still fragmentary. Here, we report that ZmMYB31, a R2R3-MYB transcription factor implicated in the regulation of lignin biosynthesis, is induced by water deficit in the growing zone of maize leaves. Consistently, we showed that ZmMYB31 colocates with quantitative trait loci for growth responses to water deficit and anthesis-silking interval (ASI) under water deficit conditions. Furthermore, we showed that increasing ZmMYB31 expression in maize had an impact on ASI that was maintained under water deficit conditions. On the basis of these data, we generated transcriptomic and proteomic data to detect transcripts and proteins that show substantial changes in abundance in the leaf growing zone of three representative independent transgene-positive homozygous maize lines compared with that of transgene-negative sibling plants grown in the PhenoArch platform under well-watered and water deficit conditions. We identified 352 genes that were transcriptionally regulated by ZmMYB31 under well-watered conditions and displayed opposite transcription levels under water deficit conditions, suggesting that ZmMYB31 may function in a water stress memory response pathway, which likely helps the transgenic maize enduring the water deficit stress. Our studies further revealed co-expression relationships between these genes and others from the proteome and a cell wall data sets, suggesting potential regulatory influence.

Funding acknowledgment: French National Research Agency (ANR-10-BTBR-01-Amaizing), FranceAgrimer and LabEx Saclay Plant Sciences-SPS (ANR-10-LABX-0040-SPS).

  Sud, Université Evry, Université Paris-Saclay, Bâtiment 630, 91405 Orsay, France; 5 Institute of Plant Sciences Paris-Saclay IPS2, Paris-Diderot, Sorbonne Paris-Cité, Bâtiment 630, 91405 Orsay, France; 6 UMR MIA-Paris, AgroParisTech, INRA, Université Paris-Saclay, 75005 Paris, France; 7 Biogemma, Centre de Recherche, 63720 Chappes, France sylvie.coursol@inra.fr

	1
	1 Institut Jean-Pierre Bourgin, INRA, AgroParisTech, CNRS, Université Paris-Saclay, RD10, 78026
	Versailles Cedex, France; 2 LEPSE, INRA, Université Montpellier, 34060 Montpellier Cedex 1,
	France; 3 GQE-Le Moulon, INRA, Université Paris-Sud, CNRS, AgroParisTech, Université Paris-
	Saclay, 91190 Gif-Sur-Yvette, France; 4 Institute of Plant Sciences Paris-Saclay IPS2, CNRS, INRA,
	Université Paris-