Michael S Macfadden
email: macfadden.michael@solute.us

Ning Gu
email: ninggu@fudan.edu.cn

Claudia-Lavinia Ignat

Haifeng Shen
email: haifeng.shen@flinders.edu.au

David Sun
email: davidqsun@gmail.com

Chengzheng Sun
email: czsun@ntu.edu.sg

The Fourteenth International Workshop on Collaborative Editing Systems

Keywords: Collaborative editing systems, groupware. ACM Classification C.2.4 [Computer-Communication Networks]: Distributed Systems -Distributed applications, H.1.2 [Models and Principles]: User/Machine Systems - Information Systems Applications]: Office Automation -Groupware, H.5.3 [Information Interfaces and Presentation]: Group and Organization Interfaces, I.7.1 [Document and text processing]: Document and text editing Algorithms, Design, Human factors, Performance, Security, and Theory

Interest in collaborative editing (CE) has seen a dramatic rise in recent years. The ubiquity of cloud services, crowdsourcing, and mobile devices means that today's Internet citizens are increasingly accustomed to producing and editing data in a shared network environment. While systems such as Google Drive, Microsoft Web Apps, Apache Wave and Codoxware allow users to collaboratively edit shared information, they have just begun to scratch the surface of CE's full potential. In the coming years, users will expect to be able to collaboratively create, share, and edit documents and data in a dynamic, realtime, and intuitive manor. This workshop aims to connect researchers, developers, and users to help explore the future of CE in tomorrow's information landscape. This year's workshop focuses on how researchers and industry practitioners can work together to accelerate delivery of CE capabilities to meet the needs of the typical information-age user.

Introduction

Collaborative Editing (CE) technologies allow users to create and edit a wide array of shared content over computer networks. CE methods have been demonstrated to support a multitude of content types such as rich text documents, spreadsheets, 3D digital media design documents, virtual environments, and multimedia documents. CE techniques focus on the technical underpinnings to maintain consistency and integrity in shared information, and also facilitate a natural and intuitive collaboration environment by providing enhanced awareness of multiple users actions. Proper application of CE systems can lead to increased efficiency in creating shared content, improved team dynamics, and higher quality end products due to increased information exchange between collaborators.

Due to their traditional cost and complexity, implementations of CE Systems have historically been the domain of seasoned computer professionals, academic researchers, and large corporations. However, many of today's average technology consumers have persistent Internet connections, extensive web experience, and mobile Internet enabled devices. Internet users are increasingly creating and sharing content through cloud services and social media sites. As users migrate their content from personal computers and storage devices to hosted services in the cloud, it has become more natural for users to share and work with each other's content. Whether it simply is tagging each other photo's on Facebook, contributing to articles in Wikipedia, or collaboratively editing documents on Google Drive, users are beginning to expect to share and modify their content in a collaborative mode.

Many of these products offer only limited collaborative functionality, including locking, asynchronous, non-realtime collaboration where only one user is able to edit a shared document at any given time. Many services also lack intuitive collaboration awareness functionality to allow users to work together effectively. Meanwhile, over the past 20+ years many research efforts have studied various key technical and social issues surrounding CE. Examples of social and technical issues in CE include the impact of culture or organizational roles on CE, melding of real-time and non-real-time CE in social interaction, privacy and security, consistency maintenance, group undo, group awareness, interaction control, and syntactic/semantic conflict resolution. Significant attention has been given to algorithmic and human factors techniques that allow users to freely and concurrently modify shared documents while minimizing and resolving conflicts. These efforts have demonstrated and documented several effective techniques relevant to user's needs of collaborative systems.

While in recent years commercial industry's participation in the CE domain has steadily increased, it has not kept pace with many users' expectations.

CSCW 2014 • Workshop Summary

February 15-19, 2014, Baltimore, Maryland, USA

Several products such as Google Drive, Novell Vibe, and Codoxware have proven that CE functionality can greatly increase the usability and user impression of mature products. These systems have begun to raise the general public's awareness of CE and demonstrated the marketability of CE functionality. These systems, however, have only begun to scratch the surface of the market and there exists a great opportunity to both increase the capability of CE systems as well increase ability to reach end users.

Workshop Theme

This workshop builds on the success of the previous thirteen annual workshops. The theme of each workshop evolved over time to focus on the current issues surrounding collaborative editing based on new research, technological advancements in computing and networking, and socio-technological changes impacting the way users interact with computer supported collaboration systems.

Last year's workshop, IWCES13 held at CSCW'12, focused on traditional technical CE topics such as theoretical foundation, algorithm design, and system/application development but in the context of complex document formats such as rich text, html/xml, spreadsheet, 2D/3D digital media, CAD, video, etc. Specific focus was given to a discussion of user requirement analysis, evaluation, and usability studies of CE systems for complex real-world problems.

To maintain the momentum generated by IWCES13, this year's workshop will focus on aligning the interests of end users, industry product developers, and academic researchers to accelerate the incorporation of collaborative editing features in future services. As users are pushing their data in to the cloud, there is a growing demand for lively, real-time collaboration features.

Workshop Participation

We encourage participation from academic and industry practitioners interested in researching various aspects of CE. Participation from industry representatives who have experience and interest in building real-world CE systems is also encouraged. Finally we welcome human factors experts, interface designers, social interaction researchers, and end-users that are interested in novel applications of CE.

Each participant is requested to submit a position paper describing the presentation or discussion they wish to contribute to the workshop. While not mandatory, those prospective participants submitting a position paper will be given priority if space becomes limited. We invite submissions pertaining to all topics relevant to CE and CE systems. Each position paper should contain 2-6 pages summarizing the content to be presented and/or discussed. While certainly not required, if your topic pertains to a system and that system is easily demonstrated, a demonstration is encouraged. All papers should be formatted using the standard ACM SIGCHI format and should include an

  The goal of the workshop is to discuss the latest research findings, review CE adoption within industry, and analyze users' expectations. Through a continuing and open dialogue we hope to facilitate transition of research results to industry, to better understand the needs and expectations of users, and to identify new areas of required research that can aid in future adoption. Topics of discussion will be a mix of technical, social, usability, and commercial implications of CE.

abstract of no more than 150 words. Topics of interest include, but are not limited to:

• Concurrency control and consistency maintenance • Visions for the future of CE.

All submissions should be emailed to Michael MacFadden at macfadden.michael@solute.us. Submissions will be reviewed by the workshop organizing committee based on relevancy, novelty, and quality. Accepted papers will be published online on the workshops web site. It is required for at least one author of an accepted paper to register and attend the workshop in order to present the material.

Workshop Agenda

The proposed workshop will be a full-day event. The workshop agenda will contain presentations from participants, focused discussion sessions, and an open forum session. The workshop will likely attract 15-20 participants.

Important Dates

• Paper submission: November 22, 2013.

• Notification of acceptance: December 6, 2013.


H H E U M A N N

J B L U M

C B O U L B E

B F A U G E R

G S E L I G

J. -M A N É2

S B R É M O N D

V G R A N D G I R A R D

E N A R D O N

Quasi-static Free-Boundary Equilibrium of Toroidal Plasma with CEDRES++: Computational Methods and Applications

We present a comprehensive survey of the various computational methods in CEDRES++ for finding equilibria of toroidal plasma. Our focus is on free-boundary plasma equilibria, where either poloidal field coil currents or the temporal evolution of voltages in poloidal field circuit systems are given data. Centered around a piecewise linear finite element representation of the poloidal flux map, our approach allows in large parts the use of established numerical schemes. The coupling of a finite element method and a boundary element method gives consistent numerical solutions for equilibrium problems in unbounded domains. We formulate a new Newton method for the discretized nonlinear problem to tackle the various non-linearities, including the free plasma boundary. The Newton method guarantees fast convergence and is the main building block for the inverse equilibrium problems that we can handle in CEDRES++ as well. The inverse problems aim at finding either poloidal field coil currents that ensure a desired shape and position of the plasma or at finding the evolution of the voltages in the poloidal field circuit systems that ensure a prescribed evolution of the plasma shape and position. We provide equilibrium simulations for the tokamaks ITER and WEST to illustrate the performance of CEDRES++ and its application areas.

Introduction

Computer codes that address the equilibrium of toroidal plasmas are central tools in tokamak fusion science. They are essential, both for detailed simulations with sophisticated magnetohydrodynamic (MHD) models as well as for experimenters that need to control real tokamak reactors. Detailed MHD simulations, which model the plasma on very short timescales, are used to study the various effects of turbulence and instability. They rely on a given plasma equilibrium as initial condition. Experimenters use equilibrium codes to set up discharge scenarios, to study breakdowns and disruptions, or to design the layout of new machines. They also use such codes, in connection with transport codes [START_REF] Hinton | Theory of plasma transport in toroidal confinement systems[END_REF][START_REF] Hirshman | Two-dimensional transport of tokamak plasmas[END_REF][START_REF] Artaud | The CRONOS suite of codes for integrated tokamak modelling[END_REF][START_REF] Coster | The european transport solver[END_REF][START_REF] Parail | Self-consistent simulation of plasma scenarios for iter using a combination of 1.5d transport codes and free-boundary equilibrium codes[END_REF], to design and validate plasma feedback controller for real tokamak machines and to verify the feasibility of scenarios in terms of operational limits (e.g. coil currents or forces).

Hence, equilibrium codes are essential tools for tokamak scientists, and applicants expect a certain degree of maturity and robustness. In the design of discharge scenarios or in the validation of feedback controller, for example, a robust, fast and automated computation of equilibria allows to shift the focus of research towards the difficulties of coupling with complex physics or improved control algorithms. CEDRES++ deals with equilibrium problems that are related to a quasi-static description of plasma evolution, which asserts balance of forces at each instant of time. A code that treats such quasi-static free-boundary equilibrium problems needs to solve non-linear elliptic or parabolic problems with non-linear source terms representing the current density profile, that vanishes outside the unknown free boundary of the plasma. The computational challenges in the design of free-boundary equilibrium codes are a problem setting in an unbounded domain with a non-linearity due the current density profile in the unknown plasma domain and the non-linear magnetic permeability if the reactor has ferromagnetic structures.

The simulation on the unbounded domain can be reduced to computations on a finite domain thanks to analytical Green's functions [START_REF] Lackner | Computation of ideal MHD equilibria[END_REF]). The numerical solution on the finite interior domain is coupled through boundary conditions to the Green's function representation of the solution in the exterior domain. This approach is today fairly standard in many other application areas such as electromagnetics [START_REF] Hiptmair | Coupling of finite elements and boundary elements in electromagnetic scattering[END_REF][START_REF] Zhao | Solving electromagnetic problems using a novel symmetric fem-bem approach[END_REF] or elasticity [START_REF] Costabel | Coupling of finite and boundary element methods for an elastoplastic interface problem[END_REF][START_REF] Bielak | Symmetric finite element and boundary integral coupling methods for fluid-solid interaction[END_REF][START_REF] Stephan | Coupling of finite elements and boundary elements for some nonlinear interface problems[END_REF]). The boundary element method [START_REF] Chen | Boundary element methods[END_REF][START_REF] Nédélec | Acoustic and electromagnetic equations[END_REF] is the name of this general framework. The boundary element method reduces problems on unbounded domains to problems on boundaries, that can then be coupled to any numerical method for the interior of a bounded domain.

The non-linearity due to the current profile in the unknown plasma domain poses the major difficulties according to our experience. It is a peculiarity of plasma equilibrium problems, that the domain of the plasma is an unknown. Speaking differently, the boundary of the plasma is a free boundary, defined either by a contact with a limiter which prevents the plasma from touching the vacuum vessel, or defined as being a separatrix in the case of a poloidal divertor configuration. On top of this fairly unusual kind of nonlinearity, also the current profile in the plasma itself is a non-linear function. Moreover, in the so-called iron transformer tokamaks, a third type of non-linearity appears due to the non-linear magnetic permeability. All these non-linearities will require some iterations towards the numerical solution. Simple fixed-point iterations usually suffer from very slow convergence or even fail to converge, which made researchers move towards Newton-type methods. The latter use the information of gradients, sometimes also referred to as sensitivities, to speed up the convergence, and they can converge in cases where fixed-point iterations don't converge -a very important example is vertically unstable plasmas.

There are basically two different families of solution methods for axisymmetric plasma equilibrium problems. The first family are the so-called flux or Lagrangian coordinate methods, determining the localization of level lines that have equidistant flux-values [START_REF] Lao | Reconstruction of current profile parameters and plasma shapes in tokamaks[END_REF][START_REF] Lao | Variational moment solutions to the Grad-Shafranov equation[END_REF][START_REF] Ling | The Princeton spectral equilibrium code: PSEC[END_REF][START_REF] Turkington | Multiconstrained variational problems in magnetohydrodynamics: Equilibrium and slow evolution[END_REF][START_REF] Gruber | Computation of MHD equilibria by a quasiinverse finite hybrid element approach[END_REF][START_REF] Degtyarev | Adaptive mesh computation of magnetic hydrodynamic equilibrium[END_REF][START_REF] Delucia | An iterative metric method for solving the inverse tokamak equilibrium problem[END_REF][START_REF] Jardin | Dynamic modeling of transport and positional control of tokamaks[END_REF][START_REF] Delucia | An iterative metric method for solving the inverse tokamak equilibrium problem[END_REF][START_REF] Degtyarev | An inverse variable technique in the MHDequilibrium problem[END_REF] (see also [START_REF] Jardin | Computational methods in plasma physics[END_REF], section 5.5)). A second family of methods uses standard finite difference methods on rectangular grids [START_REF] Feneberg | Multipole tokamak equilibria[END_REF][START_REF] Helton | MHD equilibrium in non-circular tokamaks with field-shaping coil systems[END_REF][START_REF] Johnson | Numerical determination of axisymmetric toroidal magnetohydrodynamic equilibria[END_REF][START_REF] Lackner | Computation of ideal MHD equilibria[END_REF] or finite element methods on triangular grids [START_REF] Blum | The self-consistent equilibrium and diffusion code SCED[END_REF]. The main difference between most methods of both of these families is the treatment of the so-called fixed boundary equilibrium problem, i.e. a problem where the plasma domain is known. The computational issues related to the unknown boundary have received less attention.

The CEDRES++ code uses a finite element formulation for the axisymmetric free-boundary equilibrium problem in the interior domain. This allows first, for standard, well established coupling methods to the boundary element formulation on the exterior domain [START_REF] Albanese | On the solution of the magnetic flux equation in an infinite domain[END_REF]). Second, we can derive a perfect Newton method, that uses the information about all non-linearities, e.g. also those related to the free-boundary setting. We consider this to be the most distinctive feature of CEDRES++ among many other equilibrium codes. Up to our knowledge there is no other equilibrium solver that uses this information to speed up the convergence. Furthermore, accurate derivatives are vital for inverse free-boundary equilibrium problems, which aim at finding the values of control parameters that ensure that the plasma attains a certain desired state, i.e. shape or position. Inverse free-boundary equilibrium problems are formulated as constrained optimization problems and only accurately computed derivatives can guarantee that the optimization algorithms find indeed the optimum. For the moment, CEDRES++ uses linear Lagrangian elements, which due to the low regularity of the solution, seem to be the obvious choice. We would like to refer to Section 5 for a general discussion on this topic. CEDRES++ inherits the basic ideas of the free-boundary equilibrium codes SCED [START_REF] Blum | The self-consistent equilibrium and diffusion code SCED[END_REF] and Proteus [START_REF] Albanese | Numerical studies of the Next European Torus via the PROTEUS code[END_REF]) but relies on object oriented and modular programming principles. CEDRES++ uses well established and tested external modules for e.g. mesh generation [START_REF] Shewchuk | Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator[END_REF], linear algebra [START_REF] Renard | GetFEM++, an open-source finite element library[END_REF] and algebraic solver [START_REF] Davis | Umfpack[END_REF]. The very first conception of CEDRES++, that used the same methods as SCED and Proteus, was developed in [START_REF] Grandgirard | Modélisation de l'équilibre d'un plasma de tokamak[END_REF]. Various simulations with this old version of CEDRES++ are reported in [START_REF] Grandgirard | Modélisation de l'équilibre d'un plasma de tokamak[END_REF] and [START_REF] Hertout | The cedres++ equilibrium code and its application to ITER, JT-60SA and Tore Supra[END_REF].

The current version of CEDRES++ contains a new module that, when coupled to a transport code, simulates a quasi-statically evolving equilibrium: the classical Grad-Shafranov equation, a non-linear elliptic partial differential equation, is satisfied at each instant of time. This mode assumes that the evolution of voltages in poloidal field circuits and the non-linearities in plasma current profile are known. The new mode is referred to as the evolution mode as opposed to the static mode that takes poloidal field coil currents and the current density as input. Within the new evolution mode, we solve the full parabolic partial differential equation system. We do not have to estimate the nonlinear mutual inductance of the plasma with the electromagnetic reactor components as the approach in [START_REF] Albanese | The linearized create-l plasma response model for the control of current, position and shape in tokamaks[END_REF] and [START_REF] Ariola | Magnetic Control of Tokamak Plasmas[END_REF], Chapter 2) would require. All the dynamics of the plasma core related to resistive diffusion of magnetic flux and transport of particle density and temperatures, are supposed to be treated by external tools and are not subject of this report. We refer to (Falchetto et al. 2014) for the coupling of CEDRES++ (Couplage Equilibre Diffusion Résistive pour l'Etude des Scénarios i.e. Coupling of Equilibrium and Resistive Diffusion for the Evaluation of Scenarios) with the transport code ETS [START_REF] Coster | The european transport solver[END_REF]. CEDRES++ is also coupled to the transport code CRONOS [START_REF] Artaud | The CRONOS suite of codes for integrated tokamak modelling[END_REF]. The evolution mode used with prescribed evolution of the current profile is also a good practical approach for vertical stability studies, where the timescale of interest is much shorter than the current diffusion timescale of the plasma.

Further, CEDRES++ can solve inverse free-boundary equilibrium problems. The inverse problem in the static mode aims at finding poloidal field coil currents that ensure a desired shape and position of the plasma. The inverse problem in the evolution mode aims at finding the evolution of the voltages in the poloidal field circuits that ensure a prescribed evolution of the plasma shape and position. We use standard algorithms for constrained optimization to solve the inverse problems. Therefore it will be straightforward to add in the near future further constraints, such as constraints on the flux consumption or the currents in the coils. In table 1 we summarize the basic CEDRES++ modes and their areas of application.

Previous implementations of the Newton method in SCED [START_REF] Blum | The self-consistent equilibrium and diffusion code SCED[END_REF] and Proteus [START_REF] Albanese | Numerical studies of the Next European Torus via the PROTEUS code[END_REF] relied on the discretization of a Newton method formulated on a continuous level. It is not clear, whether this formulation remains valid for equilibria with plasma boundaries in the case of a poloidal divertor configuration. The distinctive new feature of CEDRES++ is a Newton method, that solves the discretized non-linear equations. Our new approach has more rigorous mathematical foundations and is supposed to have slightly faster convergence. Moreover, it is only this new approach, which guarantees that the optimization algorithms for solving the inverse problems converge to the correct solution. Section 3 gives more explanations on that.

The users of CEDRES++ do not need to know about details of the algorithms and the parameters. CEDRES++ is a robust, fast and accurate and an easily usable tool. CEDRES++ focuses for the moment on the solution of the so called axisymmetric freeboundary plasma equilibrium with isotropic pressure and without flow. The assumption of perfect axial symmetry is a common model reduction in many equilibrium applications and the treatment of 3D plasma equilibria [START_REF] Park | Plasma simulation studies using multilevel physics models[END_REF][START_REF] Hirshman | Preconditioned descent algorithm for rapid calculations of magnetohydrodynamic equilibria[END_REF] requires still a lot of computational power. We are planning to include in the near future numerical methods for plasma equilibria with flow and plasma equilibria with anisotropic pressure [START_REF] Grad | Toroidal containment of a plasma[END_REF][START_REF] Maschke | An analytic solution of the stationary {MHD} equations for a rotating toroidal plasma[END_REF][START_REF] Goedbloed | Stationary symmetric magnetohydrodynamic flows[END_REF][START_REF] Zwingmann | Equilibrium analysis of tokamak discharges with anisotropic pressure[END_REF][START_REF] Guazzotto | Numerical study of tokamak equilibria with arbitrary flow[END_REF][START_REF] Cooper | Three-dimensional anisotropic pressure free boundary equilibria[END_REF][START_REF] Pustovitov | Anisotropic pressure effects on plasma equilibrium in toroidal systems[END_REF][START_REF] Fitzgerald | EFIT tokamak equilibria with toroidal flow and anisotropic pressure using the two-temperature guiding-centre plasma[END_REF]. Toroidal equilibria with anisotropic pressure and flow are an active area of research that will benefit from our contribution to the computation of free-boundary equilibria. CEDRES++ is not considered to be used as a so-called equilibrium reconstruction code [START_REF] Hofmann | Tokamak equilibrium reconstruction using Faraday rotation measurements[END_REF][START_REF] Lao | Equilibrium analysis of current profiles in Tokamaks[END_REF][START_REF] Mc Carthy | The CLISTE Interpretive Equilibrium Code[END_REF][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 relies on measurements during the discharge to compute the magnetic fields and estimates of current profiles and other characteristics of plasma equilibria.

The outline of the article is the following: In the first section we recall briefly the basic equations that describe the free-boundary plasma equilibrium in a tokamak and state the four main problems that can be solved with CEDRES++. The subsequent section contains detailed descriptions of the various numerical methods that are implemented in CEDRES++. This is followed by a short section containing tests for the numerical validation and various application examples.

Quasi-Static Free-Boundary Equilibrium of Toroidal Plasma

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

grad p = J × B , div B = 0 , curl 1 µ B = J , (2.1)
where p is the plasma kinetic pressure, B is the magnetic field, J is the current density and µ the magnetic permeability. In the quasi-static approximation these static equations are augmented by Faraday's law

-∂ t B = curl E , (2.2)
with E the electric field, and by Ohm's laws in plasma, coils and passive structures.

For the calculations in CEDRES++ we will differentiate between static problems and evolution problems, where the keyword static indicates that the equations do not give a time-varying solution. The static problems and the evolution problems are treated by CEDRES++ static mode and evolution mode, respectively. Force balance, solenoidal condition and Ampère's law in (2.1) yield the static problem we will introduce in detail in Section 2.1, while the evolution problems introduced afterwards in Section 2.3, take also into account the Faraday's and Ohm's law in the poloidal field coils and in the passive structures. All the dynamics due to Faraday's and Ohm's laws in the plasma, as well as the dynamics related to transport of heat and particles are supposed to be treated by external tools. Under the common assumption of perfect axial symmetry, it is convenient to put (2.1) and the quasi-static approximation of Maxwell's equations in a cylindrical coordinate system (r, ϕ, z) and to consider only a meridian section of the tokamak. Then, Faraday's and Ampère's laws decouple into corresponding laws in the toroidal direction and the poloidal plane and give rise to the toroidal and poloidal Poynting theorems:

div B tor × E pol µ = J pol • E pol + 1 2µ ∂ t B 2 tor (2.3) and div B pol × E tor µ = J tor • E tor + 1 2µ ∂ t B 2 pol , (2.4)
where the subscripts tor and pol indicate the toroidal and poloidal components of the fields. Poynting theorems can be used to check the conservation of energy in the simulation of quasi-static plasma equilibria, thus providing a global control of the accuracy. We introduce

Ω ∞ = [0, ∞] × [-∞, ∞],
the positive half plane, to denote the meridian plane that contains the tokamak centered at the origin. The geometry of the tokamak determines the various subdomains (see Figure 1 • the domain Ω Fe ⊂ Ω ∞ corresponds to those parts that are made of iron; for an air-transformer tokamak Ω Fe = ∅;

• the domains Ω ci,j ⊂ Ω ∞ , 1 i L, 1 j N i , correspond to the L i=1 N i = N poloidal field coils. The coils are grouped into L poloidal field circuits and the ith circuit contains N i coils. The intersection of the jth coil in the ith circuit with the poloidal plane is Ω ci,j , and it has n i,j wire turns, total resistance R i,j and cross section area S i,j ;

• the domains Ω ps k ⊂ Ω ∞ , k = 1, ..., N ps corresponding to N ps passive structures with conductivity σ k ;

• the domain Ω L ⊂ Ω ∞ , bounded by the limiter, corresponds to the domain which is accessible by the plasma;

• the domain Ω p ⊂ Ω L , is the domain covered by the plasma. The classical primal unknowns for toroidal plasma equilibria described by (2.1) are the poloidal magnetic flux ψ = ψ(r, z), the pressure p and the diamagnetic function f . The poloidal magnetic flux ψ := rA • e ϕ is the scaled toroidal component of the vector potential A, i.e. B = curl A and e ϕ the unit vector for ϕ. The diamagnetic function f = rB • e ϕ is the scaled toroidal component of the magnetic field. It can be shown that both the pressure p and the diamagnetic function f are constant on ψ -isolines, i.e. p = p(ψ) and f = f (ψ). We refer to standard text books, e.g. [START_REF] Freidberg | Ideal Magnetohydrodynamics[END_REF], [START_REF] Blum | Numerical simulation of the plasma equilibrium in a tokamak[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] Goedbloed | Advanced magnetohydrodynamics: with applications to laboratory and astrophysical plasmas[END_REF]) and [START_REF] Jardin | Computational methods in plasma physics[END_REF] for the details and state in the following paragraphs only the final equations describing the static and evolution problems solved in CEDRES++.

Direct Static Problem

Force balance, solenoidal condition and Ampère's law in (2.1) yield in axisymmetric configuration the following set of equations:

-∇ • 1 µr ∇ψ =      rp (ψ) + 1 µ0r f f (ψ) in Ω p (ψ) ; Ii,j Si,j
in Ω ci,j ; 0 elsewhere , ψ(0, z) = 0 ; lim (r,z) →+∞ ψ(r, z) = 0 ;

(2.5)

where ∇ is the gradient in the two dimensions (r, z), I i,j is the total current (in Ampère turns) in the j-th coil of the ith circuit and µ is a functional of ψ µ = µ Fe (|∇ψ| 2 r -2 ) in Ω Fe µ 0 elsewhere .

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

-∂ ∂r

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

(2.8)

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

ψ N (r, z) = ψ(r, z) -ψ ax (ψ) ψ bnd (ψ) -ψ ax (ψ) .
(2.9)

These two functions, subsequently termed S p and S f f , have, independently of ψ, a fixed domain [0, 1]. Further, in many applications, one assumes that the current profile, i.e. the function rS p + 1 µ0r S f f , is only known up to some scaling constant λ. In those cases the set of equations in (2.5) has to be augmented by an additional equation that matches the scaling with the given total plasma current I P .

Let us state the two problems that we will consider in the following.

Problem 1 (Direct static). Let S p : [0, 1] → R and S f f : [0, 1] → R be two known functions and let the currents I i,j in the coils be given. We want to find the ψ such that (2.5) holds with p (ψ) = S p (ψ N ) and f f (ψ) = S f f (ψ N ).

Problem 2 (Direct static, with given plasma current I P ). Let S p : [0, 1] → R and S f f : [0, 1] → R be two known functions and let the currents I i,j in the coils be given.

Additionally we assume that the total plasma current I P is given. We want to find ψ and λ such that (2.5) holds with p (ψ) = λS p (ψ N ) and f f (ψ) = λS f f (ψ N ) together with

I P = λ Ωp(ψ) rS p (ψ N (r, z)) + 1 µ 0 r S f f (ψ N (r, z)) drdz.
(2.10)

The functions S p and S f f are usually given as piecewise polynomial functions. Another frequent a priori model is

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

Inverse Static Problem

The direct problem in the previous section computes a free-boundary equilibrium for given coil currents. In many applications, in particular in the area of tokamak operation, the inverse problem is equally relevant: What are the currents that give a certain desired shape to the plasma? A popular approach to answer such a question is its formulation as an optimal control problem. The currents I i in the poloidal field coils are the control variables and the magnetic flux map ψ describing the equilibrium is the controlled variable. Then we introduce a cost function for the magnetic flux ψ and the coil currents I i penalizing the deviation from a desired plasma shape and position, and we minimize this cost function under the constraint that the magnetic flux ψ and the currents in the coils solve the equilibrium problem (2.5). A regularization term ensures well-posedness of the inverse problem. Here again, the current profile in the plasma is supposed to be known data.

In CEDRES++ we prescribe a plasma state by a desired plasma boundary Γ desi . Let Γ desi ⊂ Ω L denote a closed line, contained in the domain Ω L that is either smooth and touches the limiter at one point or has at least one corner. The former case prescribes a desired plasma boundary that touches the limiter. The latter case aims at a plasma with X-point lying entirely in the interior of Ω

L . Further let (r desi , z desi ) ∈ Γ desi and (r 1 , z 1 ), . . . (r N desi , z N desi ) ∈ Γ desi be N desi + 1 points on that line. We define a quadratic cost functional K(ψ) that evaluates to zero if Γ desi is a ψ-isoline, i.e. if ψ is constant on Γ desi : K(ψ) := 1 2 N desi i=1 ψ(r i , z i ) -ψ(r desi , z desi ) 2 .
(2.12)

Another quadratic cost functional, that will serve as regularization, is

R(I 1,1 , . . . , I L,N L ) := L i=1 Ni j=1 w i,j 2 I 2 i,j .
(2.13)

The coefficients w i,j 0 are called regularization weights.

Let us state the two inverse problems that we will consider in the following.

Problem 3 (Inverse static). Let S p : [0, 1] → R and S f f : [0, 1] → R be two known functions. We solve the following minimization problem:

min ψ,I1,1,...I L,N L K(ψ) + R(I 1,1 , . . . I L,N L ) subject to (2.5) (2.14) with p (ψ) = S p (ψ N ) and f f (ψ) = S f f (ψ N ).
Problem 4 (Inverse static, with given plasma current I P ). Let S p : [0, 1] → R and S f f : [0, 1] → R be two known functions and assume additionally that the total plasma current I P is given. We solve the following minimization problem:

min λ,ψ,I1,1,...I L,N L K(ψ) + R(I 1,1 , . . . I L,N L ) subject to (2.5) and (2.10) (2.15) with p (ψ) = λS p (ψ N ) and f f (ψ) = λS f f (ψ N ).
Clearly, it is also possible to define other cost functions forcing the plasma to have other characteristics. CEDRES++ can be easily extended in this direction. Furthermore it is possible to add both equality and inequality constraints. We are planning to include for example upper and lower bounds on the currents and the forces in the coils.

Another class of inverse problems related to static equilibrium, appears in real time tokamak control. There, it is important to reconstruct both the plasma boundary as well as the current profile functions p and f f in the plasma from external measurements. Frequent and fast prediction of the current state of the plasma in the tokamak machine are essential information for feedback control system. Hence, the computational challenges in solving these inverse problems are much different, and lead to the development of a separate class of equilibrium codes [START_REF] Hofmann | Tokamak equilibrium reconstruction using Faraday rotation measurements[END_REF][START_REF] Lao | Equilibrium analysis of current profiles in Tokamaks[END_REF][START_REF] Mc Carthy | The CLISTE Interpretive Equilibrium Code[END_REF][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].

Direct Evolution Problem

In contrast to the static problems, the evolution problems in CEDRES++ take also into account the Faraday's and Ohm's laws in the poloidal field coils and in the passive structures. The dynamics due to Faraday's and Ohm's law in the plasma, as well as the dynamics related to transport of heat and particles are supposed to be treated by external tools. Alternatively, one can prescribe the profiles S p and S f f as functions of time. The Poynting theorems (2.4) and (2.3) could provide a global mean to check whether the coupling between CEDRES++ and such external tools is accurate. However, due to discretization, one needs to resort to integrated versions of the Poynting theorems for the accuracy check. Later, in section 3.6, we will present detailed formulas of such integrated Poynting theorems.

The N poloidal field coils are gathered into L poloidal field circuits which contain in total M supplies. Each of the L poloidal field circuits contains a subset of the N coils and a subset of the M supplies. We denote by I i the vector of size M i + N i which contains the currents at the M i supplies and in the N i coils of the circuit with index i, 1 i L. The circuit equations in the ith circuit can be written in the form:

I i = S i V i + R i Ψ i (∂ t ψ) , (2.16)
where the matrices S i ∈ R (Mi+Ni)×Mi and R i ∈ R (Mi+Ni)×Ni depend on the wire turns n i,• , the total resistances R i,• and the cross sections S i,• of the poloidal field coils in the circuit i and on the topology of the circuit. Details on the computation of matrices S i and R i are given in Appendix (A). The vectors V i ∈ R Mi contain the voltages applied to the supplies, and the vectors

Ψ i (ψ) ∈ R Ni are Ψ i (ψ) = ( Ωc i,1 ψdrdz, . . . , Ωc i,N i ψdrdz) T .
In the case of a simple circuit composed of one supply connected to the coil Ω ci,1 the equation (2.16) writes

I i,1 = n i,1 V i,1 (t) R i,1 -2π n 2 i,1 R i,1 S i,1 Ωc i,1 ∂ψ ∂t drdz .
The free-boundary equilibrium problem on the time interval [0, T ] for the time dependent poloidal flux ψ = ψ(t) = ψ(r, z, t) is:

-∇ • 1 µr ∇ψ =            rp (ψ, t) + 1 µ0r f f (ψ, t) in Ω p (ψ) ; S -1 i,j S i V i + R i Ψ i (∂ t ψ) j in Ω ci,j , 1 i L, 1 j N i ; -σ k r ∂ψ ∂t in Ω ps k ; 0 elsewhere , ψ(0, z, t) = 0 ; lim (r,z) →+∞ ψ(r, z, t) = 0 ; ψ(r, z, 0) = ψ 0 (r, z) ,
(2.17) The equation in the passive structures Ω ps k is deduced from Ohm's law and Faraday's law, σ k being the equivalent axi-symmetric conductivity.

Problem 5 (Evolution, direct). Let S p : [0, 1] × [0, T ] → R and S f f : [0, 1] × [0, T ] → R be two known functions. Let the evolution of the voltages V 1 (t), . . . , V L (t)
in the poloidal field circuits and the initial data ψ 0 be given. We want to find the evolution of ψ(t) such that (2.17) holds with p (ψ(t), t) = λS p (ψ N (t), t) and f f (ψ(t), t) = λS f f (ψ N (t), t).

Problem 6 (evolution, with given plasma current I P (t), direct). Let S p : [0, 1]×[0, T ] → R and S f f : [0, 1]×[0, T ] → R be two known functions. Let the evolution of the voltages V 1 (t), . . . , V L (t) in the poloidal field circuits and the initial data ψ 0 be given. Additionally we assume that the evolution of the total plasma current I P (t) is given. We want to find the evolution of ψ(t) and λ(t) such that (2.17) holds with p (ψ(t), t) = λS p (ψ N (t), t) and f f (ψ(t), t) = λS f f (ψ N (t), t) together with

I P (t) = λ(t) Ωp(ψ(t)) rS p (ψ N (r, z, t), t) + 1 µ 0 r S f f (ψ N (r, z, t), t) drdz. (2.18)
To model a consistent quasi-static evolution of plasma equilibrium the equations in (2.17) have to be augmented by the diffusion of density, temperature and magnetic flux. In that case, both functions S p and S f f in the Problems 5 and 6 appear as unknowns of the full system of equations.

Inverse Evolution Problem

The inverse evolution problem is the problem of determining external voltages such that the evolution of the plasma has certain prescribed properties. We will state this problem again as an optimal control problem.

Let Γ desi (t) ⊂ Ω L denote the evolution of a closed line, contained in the domain Ω L that is either smooth and touches the limiter at one point or has at least one corner. The former case prescribes a desired plasma boundary that touches the limiter. The latter case aims at a plasma with X-point that is entirely in the interior of Ω L . Further let (r desi (t), z desi (t)) ∈ Γ desi (t) and (r 1 (t), z 1 (t)), . . . , (r N desi (t), z N desi (t)) ∈ Γ desi (t) be N desi + 1 points on that line. We define a quadratic functional K(ψ) that evaluates to zero if Γ desi (t) is an ψ(t)-isoline, i.e. if ψ(t) is constant on Γ desi (t):

K(ψ(t)) := 1 2 T 0 N desi i=1 ψ(r i (t), z i (t), t) -ψ(r desi (t), z desi (t), t) 2 dt. (2.19)
Another functional, that will serve as regularization, is

R( V 1 (t), . . . , V L (t)) := L i=1 w i 2 T 0 V i (t) • V i (t)dt.
(2.20)

It penalizes the strength of the voltages V i and represents the energetic cost in the coil system. The coefficients w i 0 are called regularization weights.

Problem 7 (Evolution, inverse). Let S p : [0, 1]×[0, T ] → R and S f f : [0, 1]×R → R be two known functions. We solve the following minimization problem:

min ψ(t), V1(t),... V L (t) K(ψ(t)) + R( V 1 (t), . . . V L (t)) subject to (2.17) (2.21) with p (ψ(t), t) = S p (ψ N (t), t) and f f (ψ(t)) = S f f (ψ N (t), t).
Problem 8 (Evolution, with given plasma current

I P (t), inverse). Let S p : [0, 1] × [0, T ] → R and S f f : [0, 1] × R → R
be two known functions. Additionally we assume that the evolution of the total plasma current I P (t) is given. We solve the following minimization problem:

min λ(t),ψ(t), V1(t),... V L (t) K(ψ(t)) + R( V 1 (t), . . . V L (t)) subject to (2.17) and (2.18) (2.22) with p (ψ(t), t) = λ(t)S p (ψ N (t), t) and f f (ψ(t)) = λ(t)S f f (ψ N (t), t).

Computational Methods and Applications

The main challenges for solving the Problems 1-8 numerically are their formulation on an infinite domain, the non-linear right-hand side, the non-linear permeability in iron and the non-linearity due to the free plasma boundary. In the following we will use finite element methods [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] to discretize the problems 1-8, and we will see that this approach is flexible enough to tackle all those challenges at once. First, finite element methods are favored approximation methods due to their flexibility on domains with complex geometry. Second, they allow for a straight forward implementation of Newton methods to handle the strong non-linearities related to the free boundary setting. The convergence speed of such Newton methods is superior to the convergence speed of fixed-point approaches that are otherwise applied for such kind of problems. As a variational formulation is the starting point for any finite element method the section starts with the variational formulations of Problems 1, 2, 5 and 6. The subsequent paragraph on the spatial discretization, a standard finite element method with linear Lagrangian basis functions on triangles, focuses mainly on the special treatment of the free plasma boundary. It gives the important formulas, required to derive the new Newton method afterwards. Having these Newton methods at hand, it is straightforward to tackle the inverse problems. The overview on the computational methods finishes with two paragraphs describing the interfaces of CEDRES++ for the coupling with transport codes and presenting volume integrated Poynting theorems.

Variational Formulation on the Truncated Domain

We chose a semi-circle Γ of radius ρ Γ surrounding the iron domain Ω Fe , the coil domains Ω ci,j and the passive structures domain Ω ps k . The truncated domain, we use for our computations, is the domain Ω ⊂ Ω ∞ having the boundary ∂Ω = Γ ∪ Γ r=0 , where Γ r=0 := {(r, z) , r = 0}. The variational formulations of Problems 1, 2, 5 and 6 use the following Sobolev space:

V := ψ : Ω → R, Ω ψ 2 r drdz < ∞, Ω (∇ψ) 2 r -1 drdz < ∞, ψ |Γr=0 = 0 ∩C 0 ( Ω) (3.1)
and are obtained by multiplying equations in 1, 2, 5 and 6 by testfunctions ξ ∈ V and integrating by parts over Ω. They are called the variational formulations since they are the Euler equations of the mininization of the energy. Then we define

• two mappings A : V × V → R and J p : V × V → R that are linear in the last argument:

A(ψ, ξ) := Ω 1 µ(ψ)r ∇ψ • ∇ξ drdz J p (ψ, ξ) := Ωp(ψ) rp (ψ) + 1 µ 0 r f f (ψ) ξ drdz (3.2) • two bilinear forms j ps , j c : V × V → R j ps (ψ, ξ) := - Nps i=1 Ωps k σ k r ψξ drdz j c (ψ, ξ) := L i=1 Ni j=1 S -1 i,j R i Ψ i (ψ) j Ωc i,j
ξ drdz

(3.3) • N bilinear mappings i,j : R × V → R: i,j (I, ξ) := S -1 i,j I Ωc i,j ξ drdz (3.4)
• a bilinear form c : V × V → R, accounting for the boundary conditions at infinity [START_REF] Albanese | On the solution of the magnetic flux equation in an infinite domain[END_REF]):

c(ψ, ξ) := 1 µ 0 Γ ψ(P 1 )N (P 1 )ξ(P 1 )dS 1 + 1 2µ 0 Γ Γ (ψ(P 1 ) -ψ(P 2 ))M (P 1 , P 2 )(ξ(P 1 ) -ξ(P 2 ))dS 1 dS 2 .
(3.5) with

M (P 1 , P 2 ) = k P1,P2 2π(r 1 r 2 ) 3 2 2 -k 2 P1,P2 2 -2k 2 P1,P2 E(k P1,P2 ) -K(k P1,P2 ) N (P 1 ) = 1 r 1 1 δ + + 1 δ - - 1 ρ Γ and δ ± = r 2 1 + (ρ Γ ± z 1 ) 2 ,
where P i = (r i , z i ) and K and E the complete elliptic integrals of first and second kind, respectively and

k Pj ,P k = 4r j r k (r j + r k ) 2 + (z j -z k ) 2 .
We refer to (Grandgirard 1999, Chapter 2.4) for the details of the derivation. The bilinear form c(•, •) follows basically from the so called uncoupling procedure in [START_REF] Gatica | The uncoupling of boundary integral and finite element methods for nonlinear boundary value problems[END_REF] for the usual coupling of boundary integral and finite element methods. In our case, it can be shown that for all P 1 , P 2 the integral term (ψ(P 1 )-ψ(P 2 ))M (P 1 , P 2 )(ξ(P 1 )ξ(P 2 )) remains bounded. The Green's function that is used in the derivation of the boundary integral method for our problem was used earlier in finite difference methods for the Grad-Shafranov-Schlüter equations [START_REF] Lackner | Computation of ideal MHD equilibria[END_REF]). We derive the following variational formulations of the direct Problems 1 and 2.

Variational Formulation 9 (Static). Let S p : [0, 1] → R and S f f : [0, 1] → R be two known functions and let the currents I i,j in the coils be given. We set

p (ψ) = S p (ψ N ) and f f (ψ) = S f f (ψ N ) in (3.2). We want to find ψ ∈ V such that A(ψ, ξ) -J p (ψ, ξ) + c(ψ, ξ) = L i=1 Ni j=1 i,j (I i,j , ξ) (3.6)
holds for all ξ ∈ V .

Variational Formulation 10 (Static, with given plasma current I P ). Let S p : [0, 1] → R and S f f : [0, 1] → R be two known functions and let the currents I i,j in the coils be given. We set p (ψ) = S p (ψ N ) and f f (ψ) = S f f (ψ N ) in (3.2). Additionally we assume that the total plasma current I P is given. We want to find ψ ∈ V and λ ∈ R such that

A(ψ, ξ) -λ J p (ψ, ξ) + c(ψ, ξ) = L i=1 Ni j=1 i,j (I i,j , ξ), I p -λ J p (ψ, 1) = 0, (3.7)
holds for all ξ ∈ V .

The variational formulation of the evolution problems 5 and 6 is based on an implicit Euler timestepping scheme 0 := t 0 < t 0 + ∆t 1 = t 1 < . . . t n-1 + ∆t n = t n = T . Other choices are possible. Since we will anyway employ only low order spatial discretization, the implicit Euler is the obvious choice.

Variational Formulation 11 (Evolution

). Let S p : [0, 1] × [0, T ] → R and S f f : [0, 1] × [0, T ] →
R be two known functions. Let the evolution of the voltages V 1 (t), . . . , V L (t) in the poloidal field circuits and the initial data ψ 0 be given. We set

p (ψ) = S p (ψ N (t), t) and f f (ψ) = S f f (ψ N (t), t) in (3.2). We want to find ψ k ∈ V approximating ψ(t k ) such that ∆t k A(ψ k , ξ) -∆t k J k p (ψ k , ξ) -j ps (ψ k , ξ) -j c (ψ k , ξ) + ∆t k c(ψ k , ξ) = ∆t k L i=1 Ni j=1 i,j ((S i V i (t k )) j , ξ) -j ps (ψ k-1 , ξ) -j c (ψ k-1 , ξ) , ψ 0 = ψ 0 , (3.8) holds for all ξ ∈ V with J k p (•, •) = J p (•, •) |t=t k .
Variational Formulation 12 (Evolution, with given plasma current I P (t)). Let S p : [0, 1] × [0, T ] → R and S f f : [0, 1] × [0, T ] → R be two known functions. Let the evolution of the voltages V 1 (t), . . . , V L (t) in the poloidal field circuits and the initial data ψ 0 be given. We set p (ψ) = S p (ψ N (t), t) and f f (ψ) = S f f (ψ N (t), t) in (3.2). Additionally we assume that the evolution of the total plasma current I P (t) is given. We want to find ψ k ∈ V and λ k ∈ R approximating ψ(t k ) and λ(t k ) such that

∆t k A(ψ k , ξ) -∆t k λ k J k p (ψ k , ξ) -j ps (ψ k , ξ) -j c (ψ k , ξ) + ∆t k c(ψ k , ξ) = ∆t k L i=1 Ni j=1 i,j ((S i V i (t k )) j , ξ) -j ps (ψ k-1 , ξ) -j c (ψ k-1 , ξ), I p (t k ) -λ k J k p (ψ k , 1) = 0, ψ 0 = ψ 0 .
(3.9)

holds for all ξ ∈ V with J k p (•, •) = J p (•, •) |t=t k .

A Galerkin Discretization

We use a standard linear Lagrangian finite element to discretize the non-linear operators in the previous section. Finite element methods are particularly well suited to treat complex geometries, such as the one of the tokamak (plasma, passive structures, poloidal field coils.) We refer to section 5 for a general discussion on the choice of the order of the finite element method. For this we introduce a triangulation Ω h of the domain Ω that resolves the subdomains Ω L , Ω Fe , Ω ci,j , Ω ps k . The finite element approximation ψ h of ψ in the finite element space V h is an expansion in basis functions λ i :

ψ h (r, z) = i ψ i λ i (r, z) with ψ i ∈ R.
(3.10) Each Lagrangian basis function λ i (r, z) is piecewise linear and vanishes at all vertices except one. The domain of the plasma Ω p (ψ h ) of a finite element function ψ h is bounded by a continuous, polygonal, closed line. The critical points (r bnd (ψ h ), z bnd (ψ h )) and (r ax (ψ h ), z ax (ψ h )) are the coordinates of certain vertices of the mesh. The saddle point of a piecewise linear function ψ h is some vertex (r 0 , z 0 ) with the following property: if (r 1 , z 1 ), (r 2 , z 2 ) . . . (r n , z n ), denote the counterclockwise ordered neighboring vertices the sequence of discrete gradients ψ 0 -ψ 1 , ψ 0 -ψ 2 . . . ψ 0 -ψ n changes at least four times the sign. It remains to specify the quadrature rule that is used to approximate integrals over triangles T and integrals over intersection T ∩Ω p (ψ h ) of triangles with the plasma domain

T f (r, ψ h )λ i drdz and T ∩Ωp(ψ h ) g(r, ψ h )λ i drdz .
(3.11)

The second type of integrals appears in J P due to the fact that the mesh does not resolve the boundary of the plasma domain Ω p . In any case we will use the centers of gravity

b T := (r T , z T ) and b T (ψ h ) := (r T (ψ h ), z T (ψ h )) := (r T ∩Ωp(ψ h ) , z T ∩Ωp(ψ h ) ) (3.12)
of the integration domains T or T ∩ Ω p (ψ h ) as quadrature points. The corresponding quadrature weights are the size of the corresponding domain |T | and |T ∩ Ω p (ψ h )|. The barycenter for the second type of integrals depends itself on ψ h . Our choice of quadrature rule introduces a consistency error of order O(h 2 ), where h is the diameter of the triangle, i.e. the quadrature is exact for linear integrands.

For a triangle T with vertex coordinates a i , a j , a k ∈ R 2 the center of gravity corresponds to the barycenter: If the domain of integration is T ∩ Ω p (ψ h ), we have to distinguish the two cases, where T ∩ Ω p (ψ h ) is either a triangle or a quadrilateral. Without loss of generality we assume that a i , a j , a k is a counterclockwise ordering of the vertex coordinates of T and that ∂Ω p (ψ) intersects ∂T at two points m k and m j at the edges opposite to the vertices a k and a j (See Figure 2). The barycentric coordinates of the intersecting points m k and m j are functions of ψ h :

(r T , z T ) = 1 3 (a i + a j + a k ). (3.13) a i a k a j ∂Ω p m k m j a j a k a i ∂Ω p m k m j
λ j (m k (ψ h )) = ψ bnd (ψ h ) -ψ i ψ j -ψ i , λ k (m j (ψ h )) = ψ bnd (ψ h ) -ψ i ψ k -ψ i , (3.14)
and, clearly, we have 2, left) we find:

λ k (m k (ψ h )) = λ j (m j (ψ h )) = 0. If T ∩ Ω p (ψ h ) is a triangle and a i that vertex of T that is contained in T ∩ Ω p (ψ h ) (See Figure
(r T (ψ h ), z T (ψ h )) = a i + 1 3 λ j (m k (ψ h ))(a j -a i ) + 1 3 λ k (m j (ψ h ))(a k -a i ) (3.15)
and 2, right) we find:

|T ∩ Ω p (ψ h )| = |T |λ j (m k (ψ h ))λ k (m j (ψ h )). (3.16) If T ∩Ω p (ψ h ) is a quadrilateral and a i that vertex of T that is not contained in T ∩Ω p (ψ h ) (See Figure
(r T (ψ h ), z T (ψ h )) = a i + 1 3 1 -λ 2 j (m k (ψ h ))λ k (m j (ψ h )) 1 -λ j (m k (ψ h ))λ k (m j (ψ h )) (a j -a i ) + 1 3 1 -λ j (m k (ψ h ))λ 2 k (m j (ψ h )) 1 -λ j (m k (ψ h ))λ k (m j (ψ h )) (a k -a i ) (3.17)
and

|T ∩ Ω p (ψ h )| = |T | (1 -λ j (m k (ψ h ))λ k (m j (ψ h ))) . (3.18)
In the next paragraph we will present a Newton method for the discretized non-linear problems, and it is important to work out accurately all the non-linear dependencies on the finite element solution ψ h . Only then we can compute the correct derivatives.

For the sake of brevity we do not write down explicitly the discrete versions of the operators from the previous paragraph, but introduce the subcript h to denote the discretized non-linear operators. A h for example is the discretized version of A. We get fully discrete non-linear formulations.

Galerkin Formulation 13 (Static). Let S p : [0, 1] → R and S f f : [0, 1] → R be two known functions and let the currents I i in the coils be given. We set p (ψ) = S p (ψ N ) and f f (ψ) = S f f (ψ N ) in (3.2). We want to find ψ h ∈ V h such that

A h (ψ h , ξ h ) -J p,h (ψ h , ξ h ) + c h (ψ h , ξ h ) = L i=1 Ni j=1 i,j,h (I i,j , ξ h ) (3.19)
holds for all ξ h ∈ V h .

Galerkin Formulation 14 (Static, with fixed plasma current I P ). Let S p : [0, 1] → R and S f f : [0, 1] → R be two known functions and let the currents I i in the coils be given. We set p (ψ) = S p (ψ N ) and f f (ψ) = S f f (ψ N ) in (3.2). Additionally we assume that the total plasma current I P is given. We want to find ψ h ∈ V h and λ ∈ R such that

A h (ψ h , ξ h ) -λ J p,h (ψ h , ξ h ) + c h (ψ h , ξ h ) = L i=1 Ni j=1 i,j,h (I i,j , ξ h ), I p -λ J p,h (ψ h , 1) = 0, (3.20) holds for all ξ h ∈ V h .
Galerkin Formulation 15 (Evolution). Let S p : [0, 1] × [0, T ] → R and S f f : [0, 1] × [0, T ] → R be two known functions. Let the evolution of the voltages V i (t) in the poloidal field circuits and the initial data ψ 0 be given. We set p (ψ) = S p (ψ N (t), t) and f f (ψ) = S f f (ψ N (t), t) in (3.2). We want to find

ψ k h ∈ V h approximating ψ(t k ) such that ∆t k A h (ψ k h , ξ h ) -∆t k J k p,h (ψ k h , ξ h ) -j ps h (ψ k h , ξ h ) -j c h (ψ k h , ξ h ) + ∆t k c(ψ k h , ξ h ) = ∆t k L i=1 Ni j=1 i,j,h ((S i V i (t k )) j , ξ h ) -j ps h (ψ k-1 h , ξ h ) -j c h (ψ k-1 h , ξ h ), ψ 0 h = ψ 0 (3.21) holds for all ξ h ∈ V h with J k p,h (•, •) = J p,h (•, •) |t=t k .
Galerkin Formulation 16 (Evolution, with given plasma current

I P (t)). Let S p : [0, 1] × [0, T ] → R and S f f : [0, 1] × [0, T ] →
R be two known functions. Let the evolution of the voltages V i (t) in the poloidal field circuits and the initial data ψ 0 be given. We set p (ψ) = S p (ψ N (t), t) and f f (ψ) = S f f (ψ N (t), t) in (3.2). Additionally we assume that the evolution of the total plasma current I P (t) is given. We want to find

ψ k h ∈ V h and λ k ∈ R approximating ψ(t k ) and λ(t k ) such that ∆t k A h (ψ k h , ξ h ) -∆t k J k p,h (ψ k h , ξ h ) -j ps h (ψ k h , ξ h ) -j c ,h (ψ k h , ξ h ) + ∆t k c(ψ k h , ξ h ) = ∆t k L i=1 Ni j=1 c i,j,h ((S i V i (t k )) j , ξ h ) -j ps h (ψ k-1 h , ξ h ) -j c ,h (ψ k-1 h , ξ h ), I p (t k ) -λ k J k p,h (ψ k h , 1) = 0, ψ 0 h = ψ 0 , (3.22) holds for all ξ ∈ V h with J k p,h (•, •) = J p,h (•, •) |t=t k .
The Galerkin formulations assume that the function µ Fe is known. In practical applications µ Fe needs to be estimated from experimental data. We refer to [START_REF] Glowinski | Analyse numérique du champ magnétique d'un alternateur par éléments finis et sur-relaxation ponctuelle non linéaire[END_REF]) and, more recently [START_REF] Pechstein | Monotonicity-preserving interproximation of B-H-curves[END_REF], for details.

Newton's Method and the Free Plasma Boundary

Newton's methods for solving a non-linear problem F (x) = 0 for x is the following iterative scheme:

F (x i )(x i+1 -x i ) = -F (x i ) ⇔ F (x i )x i+1 = F (x i )x i -F (x i ).
(3.23)

If F is sufficiently smooth, standard theory for Newton methods asserts that this iteration converges quadratically fast to the solution x. In our case the magnetic flux ψ or its finite element approximation ψ h plays the role of the unknown x. If we want to apply this method to either our continuous non-linear variational formulations 9, 10, 11 and 12 or the discretized versions, namely the Galerkin Formulations 13, 14, 15 and 16, we need to compute derivatives of the non-linear operators.

For the continuous formulations we need to calculate all the directional derivatives D ψ A(ψ, ξ)( ψ), D ψ J p (ψ, ξ)( ψ), D ψ j ps (ψ, ξ)( ψ), D ψ j c (ψ, ξ)( ψ) and D ψ c(ψ, ξ)( ψ). This calculation is simple for the bilinear mappings j c , j ps , e.g.,

D ψ j ps (ψ, ξ)( ψ) = j ps ( ψ, ξ), D ψ j c (ψ, ξ)( ψ) = j c ( ψ, ξ), (3.24)
and the non-linear mapping A (see (2.6)):

D ψ A(ψ, ξ)( ψ) = Ω 1 µ(ψ)r ∇ ψ•∇ξ drdz-2 ΩFe µ Fe (| grad ψ| 2 r -2 ) µ 2 Fe (| grad ψ| 2 r -2 )r 3 ∇ ψ•∇ψ∇ψ•∇ξ drdz.
The remaining derivative of J p was given in (Blum 1987, Lemma I.4):

D ψ J p (ψ, ξ)( ψ) = Ωp(ψ) ∂j p (r, ψ N (ψ)) ∂ψ N ∂ψ N (ψ) ∂ψ ψ ξ drdz, - Γp(ψ) j p (r, 1)|∇ψ| -1 ( ψ -ψ(r bnd (ψ), z bnd (ψ)))ξ dΓ + Ωp(ψ) ∂j p (r, ψ N (ψ)) ∂ψ N ∂ψ N (ψ) ∂ψ ax ψ(r ax (ψ), z ax (ψ))ξ drdz, + Ωp(ψ) ∂j p (r, ψ N (ψ)) ∂ψ N ∂ψ N (ψ) ∂ψ bnd ψ(r bnd (ψ), z bnd (ψ))ξ drdz , (3.25)
where Γ p is the plasma boundary ∂Ω p and

j p (r, ψ N (ψ)) = rS p (ψ N (ψ)) + 1 µ 0 r S f f (ψ N (ψ)) . (3.26)
The derivation involves shape calculus [START_REF] Murat | Sur le contrôle par un domaine géométrique[END_REF][START_REF] Delfour | Shapes and geometries[END_REF] and the non-trival derivatives:

D ψ ψ ax (ψ)( ψ) = ψ(r ax (ψ), z ax (ψ)) and D ψ ψ bnd (ψ)( ψ) = ψ(r bnd (ψ), z bnd (ψ)) .
The formula of the derivative relies on certain smoothness assumptions on ψ. Up to our knowledge, there is no theoretical evidence that this formula holds also for plasma equilibria with boundaries that contain X-points. In particular the second term on the right hand side seems to blow up if ψ reaches a critical point. Also in [START_REF] Blum | Numerical simulation of the plasma equilibrium in a tokamak[END_REF], it is shown that the derivative of J p (ψ, ξ) in the direction ψ vanishes: D ψ J p (ψ, ξ)(ψ) = 0. Then the Newton scheme for solving Problem 10 is the following iteration: Let (ψ n , λ n ) be the solution at the n-th iteration. For given (ψ n , λ n ) we introduce the linear form:

F n (ξ) := -A(ψ n , ξ) + D ψ A(ψ n , ξ)(ψ n ) + L i=1 Ni j=1 c i,j (I i,j , ξ) = -2 ΩFe µ Fe (|∇ψ| 2 r -2 ) µ 2 Fe (|∇ψ| 2 r -2 )r 3 |∇ψ n | 2 ∇ψ n • ∇ξ drdz + L i=1 Ni j=1 c i,j (I i,j , ξ).
(3.27)

and the Newton update (ψ n+1 , λ n+1 ) is the solution of the infinite dimensional linear system

D ψ A(ψ n , ξ)(ψ n+1 ) -λ n D ψ J p (ψ n , ξ)(ψ n+1 ) + c(ψ n+1 , ξ) -J p (ψ n , ξ)λ n+1 = F n (ξ), ∀ξ ∈ V λ n D ψ J p (ψ n , 1)(ψ n+1 ) + J p (ψ n , 1)λ n+1 = I p
After each iteration we need to recompute

ψ ax (ψ n ) = ψ n (r ax (ψ n ), z ax (ψ n )), ψ bnd (ψ n ) = ψ(r bnd (ψ n ), z bnd (ψ n ))
and Ω p (ψ n ). For the computation of the initial flux function ψ 0 we choose a constant permeability in iron and replace the non-linear form J p (ψ, ξ) with some linear form Ωp j init ξ drdz, where Ω p is a given ellipse and j init a given constant current density. Hence ψ 0 is the solution to a linear problem and determines the plasma axis and the plasma boundary in the first Newton iteration. The Newton iterations for the Problems 9, 11 and 12 follow likewise. The equilibrium codes SCED [START_REF] Blum | The self-consistent equilibrium and diffusion code SCED[END_REF] and Proteus [START_REF] Albanese | Numerical studies of the Next European Torus via the PROTEUS code[END_REF]) are based on discretizations of such Newton iterations. The flux functions ψ n and ψ n+1 are approximated by finite weighted sums of finite element basis functions and the test functions ξ cycle over all test functions. In each Newton iteration, one has to invert an algebraic system whose size is equal to the number of finite element basis functions. But since it is not clear, whether the formula for the derivative of J p (ψ, ξ)( ψ) remains valid for plasma boundaries with X-points, these approaches are not very trustworthy.

In CEDRES++ we prefer to use Newton methods for the Galerkin Formulations 13, 14, 15 and 16. Such Newton methods need the directional derivatives

D ψ h A h (ψ h , ξ h )( ψ h ), D ψ h J p,h (ψ h , ξ h )( ψ h ), D ψ h j ps h (ψ h , ξ h )( ψh ), D ψ h j c h (ψ h , ξ h )( ψh ) and D ψ h c h (ψ h , ξ h )( ψh ).
Here again, this is a straightforward and simple calculation for all mappings except one: the mapping J P,h that is related to the non-linear current profile in the plasma domain. The mapping J P,h is given by

J P,h (ψ h , λ m ) = T J T P,h (ψ h , λ m ) := T |T ∩ Ω p (ψ h )| j p (b T (ψ h )) λ m (b T (ψ h )),
where

j p (b T (ψ h )) = j p (r T (ψ h ), ψ N (ψ h (b T (ψ h )), ψ ax (ψ h ), ψ bnd (ψ h ))
). The directional derivative of J P,h (ψ h , λ m ) in direction λ n is the partial derivative with respect to the expansion coefficient ψ n :

D ψ J P,h (ψ h , λ m )(λ n ) = ∂ ∂ψ n J P,h (ψ h , λ m ) = ∂ ∂ψ n J P,h ( i ψ i λ i , λ m )
Computing the derivative of each terms of J P,h (ψ h , λ m ) is a tedious application of chain and product rules. We distinguish three different cases: 2). With a slight abuse of notation we identify ψ bnd and ψ ax with the corresponding finite element expansion coefficient and use the Kronecker deltas δ n,bnd and δ n,ax .

T ∩ Ω p (ψ h ) = 0, T ∩ Ω p (ψ h ) = T and T ∩ Ω p (ψ h ) ⊂ T (see Figure
(a) T ∩ Ω p (ψ h ) = 0:

∂ ∂ψ n J T P,h (ψ h , λ m ) = 0 ; (b) T ∩ Ω p (ψ h ) = T : ∂ ∂ψ n J T P,h (ψ h , λ m ) =|T | ∂j p (r T , ψ N (b T )) ∂ψ N ∂ψ N (b T ) ∂ψ h λ n (b T ) + ∂ψ N (b T ) ∂ψ bnd δ n,bnd + ∂ψ N (b T ) ∂ψ ax δ n,ax λ m (b T ) (c) T ∩ Ω p (ψ h ) ⊂ T :
Without loss of generality we adopt the notation from section 3.2, introduce λ k j = λ j (m k ) and

λ j k = λ k (m j ) use b T to denote b T (ψ h ). We define AR = |T |λ k j λ j k if T ∩ Ω p (ψ h ) is a triangle and AR = |T |(1 -λ k j λ j k ) if T ∩ Ω p (ψ h ) is a quadrilateral. We find ∂ ∂ψ n J T P,h (ψ h , λ m ) =A T n (ψ h , λ m ) + C T n (ψ h , λ m ) + T T n (ψ h , λ m ) with • the derivative related to the area |T ∩ Ω p (ψ h )|: A T n (ψ h , λ m ) =s|T | ∂λ k j ∂ψ n λ j k + λ j k ∂λ k j ∂ψ n + ∂λ k j ∂ψ bnd λ j k + λ j k ∂λ k j ∂ψ bnd δ n,bnd j p (b T ) λ m (b T ),
where s = 1 if |T ∩ Ω p (ψ h )| is a triangle and s = -1 else.

• the derivative related to the current j p (r T , ψ N (b T )):

C T n (ψ h , λ m ) =AR ∂j p (r T , ψ N (b T )) ∂r ∂r T ∂ψ n λ n (b T ) + ∂r T ∂ψ bnd δ n,bnd λ m (b T ) + AR ∂j p (r T , ψ N (b T )) ∂ψ N ∂ψ N (b T ) ∂ψ bnd δ n,bnd + ∂ψ N (b T ) ∂ψ ax δ n,ax λ m (b T ) + AR ∂j p (r T , ψ N (b T )) ∂ψ N ∂ψ N (b T ) ∂ψ h λ n (b T ) + ∂ψ N (b T ) ∂ψ h ∇ψ h (b T ) • ∂b T ∂ψ n + ∂ψ N (b T ) ∂ψ h ∇ψ h (b T ) • ∂b T ∂ψ bnd δ n,bnd λ m (b T )
• the derivative related to the test function λ m (b T ):

T T n (ψ h , λ m ) =AR j p (b T ) ∇λ m (b T ) • ∂b T ∂ψ n + ∇λ m (b T ) • ∂b T ∂ψ bnd δ n,bnd
The derivatives of ψ N follow easily from the definition (2.9). We would like to stress that the Galerkin matrix D ψ J T P,h (ψ h , λ m )(λ n ) can be assembled in a fairly standard, i.e. element wise, fashion, provided we compute in a preprocessing step the following information for each element element T belonging to the last case: We need to know the barycententric coordinates of the intersection points λ k (m j ) and λ j (m k ), the barycenter b T (ψ h ) and the derivatives

∂λ k (mj ) ∂ψi , ∂λ k (mj ) ∂ψj , ∂λ k (mj ) ∂ψ k , ∂λ k (mj ) ∂ψ bnd , ∂λj (m k ) ∂ψi , ∂λj (m k ) ∂ψj , ∂λj (m k ) ∂ψ k , ∂λj (m k )
∂ψ bnd and ∂b T ∂ψi , ∂b T ∂ψj , ∂b T ∂ψ k , ∂b T ∂ψ bnd . All this information can be easily computed for given ψ h , ψ bnd and ψ ax using the formulas (3.14), (3.15) and (3.17). All the terms that contain the Kronecker deltas δ n,bnd or δ n,ax lead to non-local entries in the stiffness matrix. They connect the coefficients ψ i1 = ψ bnd and ψ i2 = ψ ax with all coefficients ψ j that are associated to vertices of elements that are intersected by the plasma domain Ω P (ψ h ).

The size of the algebraic systems that we need to solve in each iteration corresponds to the number of vertices of the triangulation. Even for very fine discretizations it is today possible to use direct linear solvers such as UMFPACK [START_REF] Davis | Umfpack[END_REF]. As long as the storage amount for the algebraic system does not exceed the memory, modern direct solvers will outperform in most cases an iterative solver.

Sequential Quadratic Programming for the Inverse Problems

In CEDRES++ we use the following fully discrete reformulation of the inverse Problems 3 and 4, to find optimal currents in the poloidal field coils.

Inverse Problem 17 (Static). Let S p : [0, 1] → R and S f f : [0, 1] → R be two known functions. We set p (ψ) = S p (ψ N ) and f f (ψ) = S f f (ψ N ) in (3.2).
We solve the following minimization problem:

min ψ h ,I1,1,...I L,N L K(ψ h ) + R(I 1,1 , . . . I L,N L ) subject to (3.19) .
(3.28)

Inverse Problem 18 (Static, with given plasma current I P ). Let S p : [0, 1] → R and S f f : [0, 1] → R be two known functions and assume additionally that the total plasma current I P is given. We set p (ψ) = S p (ψ N ) and f f (ψ) = S f f (ψ N ) in (3.2). We solve the following minimization problem:

min λ,ψ h ,I1,1,...I L,N L K(ψ h ) + R(I 1,1 , . . . I L,N L ) subject to (3.20).
(3.29)

The inverse Problems 17 and 18 are finite dimensional constrained optimization problems. The sequential quadratic programming (SQP) method is the fastest method for finite dimensional constrained optimization problems. We refer to the text book (Nocedal & Wright 2006, Chapter 18) for the details and explain here only the basic idea.

Both inverse Problems 17 and 18 are optimization problems of the following type min u,y (3.30) where the quadratic matrices H and K are the discretization of the cost functions K and R, the state variable y is the vector of the finite element coefficients ψ i and the scaling factor λ, the control variable u is the vector of the N currents I i in the poloidal field coils and B and F the Galerkin discretizations of (3.19) or (3.20). The Lagrange function formalism in combination with Newton-type iterations is one approach to derive the SQP-methods: the Lagrangian for (3.30) is

1 2 y T Ky + 1 2 u T Hu s.t B(y) = F (u),
L(y, u, p) = 1 2 y T Ky + 1 2 u T Hu + p T (B(y) -F (u)) (3.31)
and the solution of (3.30) is a stationary point of this Lagrangian:

Ky + D y B T (y)p = 0, Hu -D u F T (u)p = 0, B(y) -F (u) = 0 (3.32)
The superscript T indicates the adjoint operator, which corresponds to matrix transposition in the finite dimensional case. The second line in (3.32) corresponds to the optimality condition for the gradient of the reduced cost functional 1 2 u T Hu+ 1 2 y T (u)Ky(u), where y(u) is implicitly defined by B(y(u)) = F (u). This is the main reason for which gradient type methods for a corresponding unconstrained optimization problem for the reduced cost function are too expensive: one evaluation of the gradient requires the very expensive solution of the non-linear problem in the third line of (3.32). For the SQP-methods on the other hand, the overall computing time in practical examples has about the same magnitude as the computing time for solving the constraint for given control parameters.

A quasi-Newton method for solving (3.32) are iterations of the type

  K 0 D y B T (y i ) 0 H -D u F T (u i ) D y B(y i ) -D u F (u i ) 0     y i+1 -y i u i+1 -u i p i+1 -p i   = -   Ky i + D y B T (y i )p i Hu i -D u F T (u i )p i B(y i ) -F (u i )  
(3.33) We call the iterative scheme (3.33) a quasi-Newton method since we omit the second order derivatives of B and F . The quadratic convergence of Newton's method deteriorates to super-linear convergence. The number of control parameters is much smaller than the number of state coefficients. Therefore the algebraic system in (3.33) is roughly twice as large as the algebraic system of a Newton iteration of the direct problem. Hence, as in the direct case, there is today no need to use iterative linear solver.

This will be different for the inverse problems of evolving free-boundary equilibria. There the size of the algebraic system increases by a factor that corresponds to the number of timesteps. We refer for [START_REF] Blum | Optimal control for quasi-static evolution of plasma equilibrium in tokamaks[END_REF] for details and state here only the finite dimensional inverse problems that are addressed in CEDRES++. They are based on a discrete cost function K h ({ψ k h } n k=1 ):

K h ({ψ k h } n k=1 ) = n k=1 ∆t k 2 N desi i=1 ψ k h (r i , z i ) -ψ k h (r desi (t k ), z desi (t k )) 2 (3.34)
for the finite element approximation {ψ k h } n k=1 at t k and a discrete regularization function:

R h ({ V 1 (t k )} n k=1 , . . . , { V L (t k )} n k=1 ) = L i=1 w i 2 n k=1 ∆t k V i (t k ) • V i (t k ) . (3.35) for the coil voltages { V i (t k )} n k=1 . Inverse Problem 19 (Evolution). Let S p : [0, 1]×R → R and S f f : [0, 1]×R → R be two known functions. We set p (ψ) = S p (ψ N (t), t) and f f (ψ) = S f f (ψ N (t), t) in (3.2).
We solve the following minimization problem:

min {ψ k h , Vi(t k )} n k=1 K h ({ψ k h } n k=1 ) + R h ({ V 1 (t k )} n k=1 , . . . , { V L (t k )} n k=1 ) subject to (3.21) .
Inverse Problem 20 (Evolution, with given plasma current I P ). Let S p : [0, 1] × R → R and S f f : [0, 1] × R → R be two known functions and assume additionally that the total plasma current I P is given. We set p (ψ) = S p (ψ N (t), t) and f f (ψ) = S f f (ψ N (t), t) in (3.2). We solve the following minimization problem:

min {λ k ,ψ k h , Vi(t k )} n k=1 K h ({ψ k h } n k=1 ) + R h ({ V 1 (t k )} n k=1 , . . . , { V L (t k )} n k=1 ) subject to (3.22) .
We would like to highlight that the SQP-method relies on proper derivatives of the nonlinear operators B and F . In our case F is affine, hence the derivative of B remains the most difficult part. On the other hand these derivatives are exactly the same derivatives that we used for the new Newton methods. Hence the implementation of a SQP-method for the inverse problem uses the same main building blocks.

Flux Surface Averages and Geometric Coefficients

As for any equilibrium code numerous outputs can be extracted from the 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 conventions of the European Integrated Tokamak Modelling Project (Falchetto et al. 2014;ITM 2013). 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 S f f S f (ψ N ) = [(r 0 B 0 ) 2 -2(ψ bnd -ψ ax ) 1 ψN S f f (x)dx] 1/2 , (3.36)
where B 0 is the vacuum toroidal field at r = r 0 . 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 (ψ(r,z)) r
drdz and Ω ψ N is the domain bounded by the line of flux Γ ψN . The quantities φ s and ρ s are computed from the discrete ψ h for all ψ s N using a barycentric quadrature rule (cf. Section 3.2):

φ s = T S f (ψ N (b T (ψ h ))) r T (ψ h ) |T ∩ Ω ψ s N |. (3.37)
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 | (3.38)
and ( ∂V ol ∂ρ ) s = V ol s using finite differences. Following [START_REF] Blum | Numerical simulation of the plasma equilibrium in a tokamak[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). (3.39)
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 |. (3.40)
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 .

Volume integrated Poynting theorems

The subset of equations (2.1) and ( 2.2) we used in section 2 to derive the evolution problems 5 and 6 that are solved in CEDRES++ involve the poloidal Faraday and the toroidal Ampère law. Hence, the poloidal Poynting theorem (2.4) can be used to check the accuracy of the solution independently of an additional treatment of the transport equations.

We integrate the poloidal Poynting theorem (2.4) over a volume V ol that contains all the plasma of a given scenario and get, by toroidal symmetry:

∂S ∇ψ • n µ 0 r ∂ t ψds = - S∩ΩP(ψ) rp (ψ) + 1 µ 0 r f f (ψ) ∂ t ψdrdz + S ∇ψ • ∇∂ t ψ µ 0 r drdz,
(3.41) where S denotes the intersection of V ol with the poloidal plane.

If we choose S to be the domain of the plasma Ω P (ψ) then the left hand side of equation (3.41) is V loop I P , where V loop is the plasma loop voltage. The first integral of the right hand side is related to the time rate of variation of the toroidal magnetic energy and to the work done against the plasma pressure gradient. The last term of the righthand side is the time rate of variation of the poloidal magnetic field energy.

We would like to stress that the integrated poloidal Poynting theorem corresponds to a variational formulation of the Grad-Shafranov-Schlüter equations on S. Using the notation (3.2) of our variational formulation from Section 3 we remark that the two integrals on the righthand side correspond to J p (ψ, χ S ∂ t ψ) and A(ψ, χ S ∂ t ψ), where χ S is the characteristic function of S. Hence, it can be shown that the solutions ψ k of the evolution problems 15 and 16 fulfill the volume integrated Poynting theorem up to first order accuracy in the mesh size.

The volume integrated version of the toroidal Poynting theorem (2.3) together with the static inverse mode was used in [START_REF] Ané | Design of next step tokamak: Consistent analysis of plasma performance flux composition and poloidal field system[END_REF] for the optimization of ITER scenarios.

Tests and Examples

Validation and Performance

From the best of our knowledge, there does not exist analytical solutions for the free boundary equilibrium problem considered in this paper. To provide nevertheless some evidence for convergence of the method, we follow a common approach in engineering and study the convergence towards a numerical solution that is computed on a very fine mesh.

We consider a static equilibrium with a given plasma current (Problem 2) in ITER geometry. The plasma current is I p = 15.10 × 10 6 A and the current density profile is prescribed using the model (2.11) with r 0 = 6.2m α = 0.5978, β = 0.5978, γ = 1.395. With these data we solve the Galerkin Formulation 16 on a sequence of five meshes with increasing number of elements. The solution obtained on the mesh with the largest number of triangular elements, is used as a reference solution and is noted ψ ref .

In our case the reference solution ψ ref has 577415 number of unknowns and the mesh consists of 1153174 triangular elements. The reference solution is depicted in Figure 3. For each of the other four meshes, the numerical solution ψ N ukwn is evaluated at N points = 812 different points of the computational domain which are located independently of the mesh (see Figure 4). Then the relative error

E points (N ukwn ) = ( Npoints i=1 |ψ N ukwn (M i ) -ψ ref (M i )| 2 ) 1/2 ( Npoints i=1 |ψ ref (M i )| 2 ) 1/2 (4.1)
is used to quantify the convergence. The values in Figure 4 demonstrate the expected linear convergence. Similarly, in Figure 5, we monitor convergence of the plasma quantities total plasma volume V ol S , the numerical derivative V ol s (see (3.38)) and the geometric coefficient G s := |∇ρ| 2 r -2 s (see (3.39)) that are computed in a post-processing step from the numerical solution ψ h with the methods from Section 3.5. Here again (see Figure 5) we observe approximately linear convergence in the number of unknowns.

In table 2 we give the overall computing time for the previous five computations on an Intel Sandy Bridge 2.6GHz. The computing time scales linearly with the number of unknowns. Given the fact that we solve non-linear problems, the computation time is reasonably small. Application engineers can easily solve a huge amount of different scenarios in short time to do parameter studies for example. CEDRES++ is perfectly suitable in larger workflow environments, such as the European Integrated Tokamak Modelling Project (Falchetto et al. 2014;ITM 2013). One reason for such short running times is the Newton method. The convergence history of the residuum in Table 2 shows per- number of unknowns N ukwn fect quadratic convergence: we need only very few iterations to find a numerical solution solving the discrete non-linear problem within the limits of machine precision.

|V olS -V ol S,ref |/|V ol S,ref | V ol • -V ol •,ref / V ol •,ref G• -G •,ref / G •,ref slope=-1
4.2. Quasi-static plasma equilibrium simulations for WEST

The Tungsten (W) Environment in Steady-state Tokamak (WEST) project [START_REF] Bucalossi | Feasibility study of an actively cooled tungsten divertor in tore supra for ITER technology testing[END_REF] represents a major change in the magnetic configuration of Tore Supra, moving from a circular limited configuration to a diverted (or X-point) configuration. CEDRES++ is one of the main modeling tools used for the preparation of WEST. It has been employed in particular for the definition of reference equilibria, the dimensioning of the plasma vertical position feedback system, the design of the plasma shape controller, breakdown studies, disruption simulations, etc. We give below a few examples of CEDRES++ simulations for WEST. Note that Tore Supra is an iron core tokamak and that the iron is taken into account in all of these simulations. The six return arms of the iron core are represented in CEDRES++ by an axisymmetric equivalent model, which gives the 1/R shape of the return arms visible in Figure 6. We are using the experimental data for the poloidal magnetic field B pol and the poloidal magnetizing field H pol from table 3, do piecewise linear interpolation of these data and reconstruct the permeability for arbitrary magnetic field values via µ FE (B 2 pol ) = |B pol ||H pol | -1 . We present in the following sections three different examples from research for WEST that use the static direct, static evolution and inverse static modes of CEDRES++. First simulations with the inverse evolution mode are presented in [START_REF] Blum | Optimal control for quasi-static evolution of plasma equilibrium in tokamaks[END_REF]. The inverse static mode of CEDRES++ is also extremely useful in order to define and optimize reference equilibria. We will give details for WEST in a forthcoming publication.

The current-focused case: direct static mode

Figure 6 shows a typical WEST poloidal flux map calculated by CEDRES++ in current-focused mode. The X-point is visible at the bottom of the plasma. Here, CE-DRES++ solves the direct static Problem 2 with prescribed total plasma current I P = 700kA and with the parametrized current profiles S p and S f f in (2.11), using α = 1, β = 1.5, γ = 0.9, r 0 = 2.6m. The vacuum toroidal field is B 0 = 3.524T at r = 2.6m. A few output parameters are: β p = 1.70, li = 0.93, q 95 = 3.33, q 0 = 1.17.

The voltage-evolution-focused case, direct evolution mode

Starting from the equilibrium in the previous section, we run CEDRES++ in direct evolution mode to solve problem 6. We keep all the input parameters fixed and we apply a constant voltage to the coils, equal to the resistive voltage being the product of coil resistance and current, except for the divertor coils, where we perturb the resistive voltage with ∆V = +0.1 V olt in the lower coil and ∆V = -0.1 V olt in the upper coil. This is in order to trigger a vertical instability (otherwise the plasma would stay in place). The simulation is run with a time step of 20ms. Figure 7 shows the plasma boundary at intervals of 100ms. The plasma moves down and the diverted configuration is lost after a few 100ms when the plasma comes in contact with the baffle of the pumping system. Figure 7 shows the vertical position of the magnetic axis z ax as a function of time. The early evolution is exponential (as one expects) with a time constant τ z = 95ms, while the later evolution is rather linear. The error of the integrated Poynting theorem from section 3.6 is approximately 5% and decreases if the time step size and the triangle size are refined.

0 2 4 -3 -2 -1 0 1 2 3 R (m) Z (m) 1.5 2 2.5 3 -1 -0.5 0 0.5 1 R (m) Z (m)

The current-focused case, inverse static mode

We present here an example that requires to solve an inverse free-boundary equilibrium problem: the post-processing of CRONOS [START_REF] Artaud | The CRONOS suite of codes for integrated tokamak modelling[END_REF]) simulations of WEST scenarios. Indeed, CRONOS simulations are typically run in fixed-boundary equilibrium mode (using the HELENA equilibrium code [START_REF] Huysmans | Isoparametric bicubic Hermite elements for solution of the Grad-Shafranov equation[END_REF]) with a prescribed boundary geometry). One may however need to know the magnetic field outside the plasma boundary (for example in order to prepare JOREK [START_REF] Huysmans | MHD stability in X-point geometry: simulation of ELMs[END_REF] non-linear MHD simulations), or to assess whether the scenario is feasible in terms of current limits in the coils for example. These questions may be addressed by solving the inverse Problems 17 and 18, using as desired boundary the boundary used in the CRONOS simulation, and as current profile those profiles S p and S f f that are calculated by CRONOS. These profiles are shown in Figure 8, where one can notice peaks at the edge of the plasma which are characteristic for H-mode profiles, while Figure 8 the result of the CEDRES++ calculation. It can be seen that the CRONOS boundary (black dots) is well matched.

Conclusions and Perspectives

We have presented in detail the computational methods of CEDRES++. It enables to compute quasi-static equilibrium configurations, the currents in the poloidal field coils or the voltages applied in the circuits of the poloidal field system being prescribed. In its inverse mode the code computes these currents and voltages that ensure a certain prescribed plasma shape that might evolve in time.

Due to its stability and robustness, CEDRES++ is a perfect tool to be coupled with transport codes [START_REF] Hinton | Theory of plasma transport in toroidal confinement systems[END_REF][START_REF] Hirshman | Two-dimensional transport of tokamak plasmas[END_REF], so that the evo-lution of plasma equilibrium is simulated at the resistive timescale consistently with transport processes. Reciprocally the transport codes take into account the precise geometry of the magnetic flux lines. The numerical stability of such a coupling is challenging and subject of ongoing research. This is particularly important for the simulation and optimization of scenarios in new generation tokamaks. CEDRES++, when coupled to the transport codes CRONOS [START_REF] Artaud | The CRONOS suite of codes for integrated tokamak modelling[END_REF] and ETS [START_REF] Coster | The european transport solver[END_REF], is in use for simulating such self consistent plasma evolution. The evolution mode itself, when plasma current profiles are given, is a good practical approach for vertical stability studies where the timescale of interest is much shorter than the current diffusion timescale of the plasma.

Furthermore, the modular and clear structure of CEDRES++ and the emphasis on accurate Newton methods, will make CEDRES++ very useful to implement fully automatic approaches to the optimization of scenarios. It will be easy to study and predict operational limits, and to devise control strategies that circumvent such limits.

It is possible to extend the methods presented in this work to higher order finite elements. Nevertheless, there are a couple of obstacles in order to obtain entirely higher order accurate methods:

• We are solving here a non-linear elliptic problem with discontinuous coefficients (in the case of iron-transformer tokamaks) and discontinuous righthand side. The standard convergence theory for finite elements and elliptic regularity theory does not yield improved approximation results for polynomials of degree higher than 1. Nevertheless practical experience and the theoretical results in [START_REF] Feistauer | Finite element approximation of nonlinear elliptic problems with discontinuous coefficients[END_REF] for low-order approximation of non-linear problems and in [START_REF] Li | Optimal a priori estimates for higher order finite elements for elliptic interface problems[END_REF] for high-order approximation of linear problems suggest improved accuracy. From the regularity theory for magneto-statics, we know in the case of iron-transformer tokamaks (such as JET or Tore Supra) that the solution lacks of regularity in the vicinity of the iron parts, in particular close to the interfaces air/iron and the corners of subdomains with iron. In such nonregular settings, it is required to switch to the so-called hp-version of the finite element method [START_REF] Schwab | p-and hp-finite element methods : theory and applications in solid and fluid mechanics[END_REF], that uses small triangles and low order polynomials in regions with non-regular solution and large triangles and high order polynomials elsewhere.

• Moreover, the general setting suffers from a fairly large modeling error due to the experimental permeability curve and to the axisymmetric representation of the ferromagnetic circuit that in reality consists of an iron core and a certain number of nonaxisymmetric return limbs. It is not clear whether such modeling errors might surpass the discretization error.

• Higher order accuracy requires also sufficiently accurate quadrature rules in the definition of the Galerkin methods. While such higher order quadrature rules are standard for (iso-parametric) finite elements, we foresee technical difficulties alongside with a considerable increase of computational complexity for the integrals over the intersection of the plasma with triangles. We need to implement sufficiently accurate quadrature rules for polygonal domains with non-straight boundaries. On top of this, we need to implement for the Newton method the derivatives of such quadrature rules.

We are planning to investigate such topics in the near future and to compare with alternative approaches. One promising alternative might be to switch after a couple of Newton iterations to a different discretization scheme that uses separate meshes and separate polynomial degrees for the representation of the flux in the plasma domain and its exterior. The non-linear coupling of the plasma and its exterior will lead to iteration schemes that induce small variations of the two meshes from one iteration to the next. With this, both the location of the magnetic axis and X-point are a priori not limited to a finite set of points, as it is the case for linear Lagrangian elements with fixed global mesh.

Any higher order method will involve more sophisticated algorithms and it is hard to predict if the accuracy at a fixed computing time will drastically improve. If high accuracy is required at the moment, this can still be achieved with reasonable effort by simple mesh refinement. The running time of the current version of the code is not yet optimized and still it is possible to do calculations with over half a million unknowns in less than 7 minutes on a workstation.

Acknowledgments

This work, supported by the European Communities under the contract of Association between EURATOM-CEA was carried out within the framework of the Task Force on Integrated Tokamak Modelling of the European Fusion Development Agreement. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

The authors would like to thank the reviewers for their comments that help improve the manuscript. In a tokamak, the poloidal field system is made of a number of circuits each comprising a number of coils and power supplies (see Figure 9). In the following we will present the derivation of the circuit equations for one single circuit.

Let N N be the number of nodes, N C the number of coils, and N S the number of supplies of one circuit. To get a model for this circuit taking into account all those connections, the idea is to write the potential difference for each supply and each coil of the circuit. For a supply S, we have

U k -U l = V, (6.1)
where V is the applied voltage and U k and U l the potential at the nodes that enclose the supply. For a coil C, the potential difference writes

U i -U j = R I n + 2πn S C ∂ψ ∂t drdz. (6.2)
where U i and U j are the potentials at the nodes that enclose the coil, R is the total resistance of the coil C, S its cross section, n its number of wire turns and I the total current (in Ampère-turns). The average of ∂ψ ∂t over the coil C is an approximation of the discrete sum of inductive terms seen by the various turns of the coil. This approximation is perfectly valid, if one has a homogeneous distribution of turns in the coil.

We also consider Kirchhoff's current law at each node of the circuit adding N N equations to the system. To fix the potential, we suppose that U 1 = 0. Thus, we get a set of N eq = N S + N C + N N + 1 equations which can be written in the form

A U = B V + C I + D Ψ(∂ t ψ). (6.3)
where the N S first equations represent (6.1) and the following N C equations are (6.2). The last equation of the system fixes the potential U 1 . The matrices A ∈ R Neq×N N , B ∈ R Neq×N S , C ∈ R Neq×(N S +N C ) and D ∈ R Neq×N C are called potential matrix, voltage matrix, current matrix and induction matrix, respectively. The vectors U ∈ R N N , V ∈ R N S and I ∈ R N S +N C contain the electric potential at the circuit nodes, the voltages applied at the supplies and currents at the coils and the supplies. The components of the vector Ψ(∂ t ψ) ∈ R N C are the integrals C ∂ t ψdrdz over the domain that is occupied by a coil C of the circuit. For given V , I and Ψ(∂ t ψ), there is a unique U which satisfies (6.3), hence A T A is regular. We find (6.4) plug this into (6.3) and get

U = (A T A) -1 A T B V + (A T A) -1 A T C I + (A T A) -1 A T D Ψ(∂ t ψ),
E I + F V + G Ψ(∂ t ψ) = 0 , (6.5) with E = A(A T A) -1 A T C -C F = A(A T A) -1 A T B -B G = A(A T A) -1 A T D -D.
The system (6.5) of N eq equations is over determined and I can be computed using the normal equation

I = S V + R Ψ(∂ t ψ) (6.6) with S = -(E T E) -1 E T F and R = -(E T E) -1 E T G.

Figure 1 .

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

Figure 2 .

 2 Figure 2. Integration over T ∩ Ωp(ψ h ). The green dots indicate the location of the quadrature point. The integration domain T ∩ Ωp(ψ h ) is either a) empty, b) the whole element T , c) a triangular domain or quadrilateral domain.

Figure 3 .

 3 Figure 3. The reference solution for the ITER test problem.

Figure 4 .

 4 Figure 4. Left: the location of the points in the definition of Epoints. Right: convergence of Epoint.

Figure 5 .

 5 Figure 5. Convergence of plasma volume V olS, the numerical derivative V ol s = ( ∂V ol ∂ρ )s, s = 0, ...S and the geometric coefficient Gs = |∇ρ| 2 r -2 s, s = 0, ...S. The reference quantities V ol S,ref , V ol s,ref and G s,ref correspond to the quantities computed for ψ ref .

Figure 6 .

 6 Figure6. Poloidal cross section showing ψ-isolines for a WEST typical equilibrium. Left: global view; right: zoom on the inner vacuum vessel The iron is displayed in blue, the poloidal field coils in orange and the passive structures (vacuum vessel and vertical stabilization components) in grey. The black curve is the limiter curve, i.e. the domain accessible to the plasma. The red curve is the plasma boundary and the red star the magnetic axis.

Figure 7 .Figure 8 .

 78 Figure 7. Left: Plasma boundary at intervals of 100ms in a vertical instability simulation for WEST. Right: Time evolution of the vertical position of the magnetic axis zax in a vertical instability simulation for WEST.

Figure 9 .

 9 Figure 9. Example of PF-circuit

Table 1 .

 1

	CEDRES++ mode	functionality	application areas
	static, direct	simulates plasma equilibria for	reference equilibria, initial condi-
		given poloidal field coil currents	tions for short timescale plasma
			models
	static, inverse	finds the poloidal field coil cur-	preparation of scenarios
		rents which allow the best match
		with a given plasma shape at a
		fixed time
	evolution, direct	simulates the quasi-static evolu-	coupling to transport codes, de-
		tion of plasma equilibrium for	sign and test environment for
		given poloidal field circuit volt-	feedback controller, study of
		ages	breakdowns and disruptions.
	evolution, inverse	finds optimal poloidal field circuit	feedforward control for discharge
		voltages and desired evolution of	scenario optimization.
		the plasma

  S,ref , V ol s,ref and G s,ref correspond to the quantities computed for ψ ref .

	number of triangles number of unknowns	computing time (in s)	iteration relative residual
	12099 23723 58744 328693 1153174	6134 11985 29556 164887 577415	2 5 11 88 368	1 2 3 4 5	2.667473 × 10 +00 9.157459 × 10 -02 1.781645 × 10 -03 0.525234 × 10 -06 3.935226 × 10 -12

Table 2 .

 2 Left: calculation time. Right: convergence history of the Newton iteration in the calculation of ψ ref .

  aims at equipping Tore Supra with an actively cooled tungsten divertor. This

		|B pol |	|H pol |		|B pol |	|H pol |
	1	0.00	0	9	1.76	7.968 × 10 3
	2	0.50	3.833 × 10 2	10	2.06	4.821 × 10 4
	3	0.70	3.982 × 10 2	11	2.25	1.628 × 10 5
	4	0.80	4.102 × 10 2	12	3.05	8.090 × 10 5
	5	0.88	4.270 × 10 2	13	4.05	1.588 × 10 6
	6	1.00	4.703 × 10 2	14	6.05	3.178 × 10 6
	7	1.20	6.274 × 10 2	15	98.20	7.651 × 10 7
	8	1.52	2.474 × 10 3	16	10 5	7.957 × 10 10

Table 3 .

 3 The data for the poloidal magnetic field |B pol | and the poloidal magnetizing field |H pol | that is used to reconstruct the magnetic permeability.

† Email address for correspondence: holger.heumann@inria.fr