TOKAMAK PLASMA BOUNDARY RECONSTRUCTION USING TOROIDAL HARMONICS AND AN OPTIMAL CONTROL METHOD
CASTOR Team-Project, INRIA and Laboratoire J.A. Dieudonné CNRS UMR 7351
Université Nice Sophia-Antipolis, Parc Valrose, 06108 Nice Cedex 02, France
This paper proposes a new fast and stable algorithm for the reconstruction of the plasma boundary
from discrete magnetic measurements taken at several locations
surrounding the vacuum vessel. The resolution of this inverse problem takes two steps.
In the first one we transform the set of measurements
into Cauchy conditions on a fixed contour close to the measurement points.
This is done by least square fitting a truncated series of toroidal harmonic functions to the
measurements. The second step consists in solving a Cauchy problem for the elliptic equation
satisfied by the flux in the vacuum and for the overdetermined boundary conditions on
previously obtained with the help of toroidal harmonics.
It is reformulated as an optimal control problem on a fixed annular domain of external boundary
and fictitious inner boundary . A regularized Kohn-Vogelius
cost function depending on the value of the flux on and measuring the discrepency
between the solution to the equation satisfied by the flux obtained using Dirichlet conditions on
and the one obtained using Neumann conditions is minimized.
The method presented here has led to the development of a software, called VacTH-KV,
which enables plasma boundary reconstruction in any Tokamak.
Keywords: inverse problem, boundary identification, magnetic measurements
In order to control the plasma during a fusion experiment in a Tokamak it is important to be able to compute its boundary in the vacuum vessel. This boundary is deduced from the knowledge of the poloidal flux which itself is computed from a set of discrete magnetic measurements of the poloidal magnetic field and flux scattered around the vacuum vessel.
This paper presents a new fast and stable algorithm for the reconstruction of the poloidal flux in the vacuum surrounding the plasma and of the plasma boundary.
Let us first briefly recall the equations modelizing the equilibrium of a plasma in a Tokamak [Wesson:2004].
Assuming an axisymmetric configuration one considers a 2D poloidal cross section of the vacuum vessel in the system of coordinates. Throughout the manuscript bold face letters represent vectors or matrices. In this setting the poloidal flux is linked to the magnetic field through the relation and, as there is no toroidal current density in the vacuum outside the plasma, satisfies the following linear elliptic partial differential equation
where denotes the elliptic operator
The unknown plasma boundary is determined from the equation
, being the value of the flux at an X-point or the value of the flux for the outermost closed flux line inside the limiter (see Figure 1). The reconstruction of this unknown plasma boundary is a Cauchy problem for the operator . A supplementary difficulty is that one does not have directly access to a complete set of Cauchy data on an external measurement countour but only to discrete magnetic measurements from sensors not necessarily positioned on a single contour.
The proposed resolution method for this inverse problem takes two steps. In the first one we transform the set of discrete measurements into Cauchy conditions, and its normal derivative, on a fixed contour chosen close to the measurement points. This is done by least square fitting a truncated series of toroidal harmonic functions to the magnetic measurements. In principle this first step can be sufficient to determine the plasma boundary [Lee:1981, Deshko:1983, Bondarenko:1984, Alladio:1986, Kurihara:1992, VanMilligen:1994, Fischer:2011, ACL.B.Faugeras.14.2] (also see the review article [Braams:1991]). Indeed the toroidal harmonics expansion is valid anywhere in the vacuum surrounding the plasma. The method of toroidal harmonics presented in [ACL.B.Faugeras.14.2] proved to be very efficient for the WEST Tokamak configuration [Bucalossi:2011]. The equivalent of magnetic measurements were generated with the equilibrium code CEDRES++ [ACL.B.Faugeras.15.1]. Even when noise is added to the synthetic measurements the boundary reconstructions are accurate with a judicious choice of the order of the toroidal harmonics used.
However since then new numerical experiments conducted for other Tokamaks with more elongated plasma shapes than WEST and with real experimental measurements showed that due to the ill-posedness of the inverse problem, the plasma boundary reconstructed after this first step involving toroidal harmonics only can in some cases be inaccurate. Nevertheless the interpolation of the discrete magnetic measurements on a contour chosen to be close to the measurement points is allways very accurate even though the extrapolation of the flux towards the plasma boundary can suffer from some perturbations due to the singularity of internal toroidal harmonics. This is the reason why we propose the following second step which can be seen as a regularization procedure for the toroidal harmonics expansion.
The second step consists in solving a Cauchy problem for the elliptic equation satisfied by the flux and for the overdetermined boundary conditions on obtained with the help of toroidal hamonics as explained above. The proposed method consists in a reformulation of the problem as an optimal control problem on a fixed annular domain of external boundary and an fictitious inner boundary (see Figure 1) located inside the plasma domain. We aim at minimizing a regularized Kohn-Vogelius [Kohn:1985] cost function depending on the value of the flux on and measuring the discrepency between the solution to the equation satisfied by obtained using Dirichlet conditions on and the one obtained using Neumann conditions.
The analysis of this second step alone was presented by the author in [ACL.B.Faugeras.12.2] but only tested on a very academic case with synthetic Cauchy data. The algorithm proposed in this paper is a combination of the two previous works [ACL.B.Faugeras.14.2] and [ACL.B.Faugeras.12.2]. It combines the precise data interpolation provided by the toroidal harmonics method to generate Cauchy conditions on and the robustness of the optimal control approach to extrapolate these Cauchy conditions towards the plasma and compute the plasma boundary. An implementation for fast numerical computations is also proposed.
The introduction of an optimal control problem on fixed domain appears in [Blum:1987, Blum:1989] where the cost function however is not of the same type as the one proposed in this paper. The advantage of the cost function proposed in this work is that it uses the given Dirichlet and Neumann boundary conditions on the external contour in a symmetric way. For the analysis of section II.B to be valid the control variable has to be the Dirichlet conditions on the internal contour .
The use of a fixed inner contour also appears in [Feneberg:1984] where a surface current sheet modelizes the plasma and in [Hakkarainen:1984, Kurihara:1993, Kurihara:2000] in the framework of boundary integral equations which might not be easily applicable in the case of iron core Tokamaks.
The next section presents the two main steps of the proposed method implemented in the code VacTH-KV. Section III then gives details on the numerical methods used and shows some numerical results.
Ii. Overview of the Method
II.A. Step1: interpolation of discrete magnetic measurements with toroidal harmonics
The poloidal flux at any point of the vacuum surrounding the plasma is represented by
The first term represents the flux generated by the poloidal field coils. It is computed using Green functions. Each coil is modelized by a number of filaments of current of given intensity.
The second term as a solution to the equation can be uniquely decomposed in a series of toroidal harmonics. These toroidal harmonics are found from computations [Braams:1986, Fischer:2011] involving a quasi separation of variable technique which lead to define external and internal harmonics. They form a complete set of solutions to this equation in any annular domain [Fischer:2011]. In particular this implies that this expansion can be used for an iron core Tokamak such as WEST.
The term is approximated by the following truncated expansion
The toroidal harmonic functions are computed in a toroidal coordinate system which depends on the choice of a pole . The internal harmonics involve the Legendre functions of the first kind, of degree one and half-integer order, which are singular at the pole . The external harmonics involve the Legendre functions of the second kind, of degree one and half-integer order, which are singular on the axis . If we denote by the vector of the coefficients of the expansion we have an analytic expression for . Its gradient can be computed analytically as well leading to an expression for the magnetic field in terms of , . It remains to identify from the set of magnetic measurements available. These discrete magnetic measurements are of three types:
B probes provide measurements of the poloidal field at points and directions given by unit vectors such that where the dot represents the scalar product
Flux loops provide flux measurements at points such that
Saddle loops provide flux variations measurements between two points such that
The unknown vector of coefficients is computed minimizing the least square cost function
Where and the known contributions of the poloidal field coils are substracted from the measurements:
The weights and , correspond to the assumed measurement errors. is quadratic and is minimized by solving the normal equation giving the optimal set of coefficients .
Once is computed an approximation of the flux can be obtained at any point of the vacuum surrounding the plasma by
In particular one can evaluate and its normal derivative in order to provide Cauchy boundary conditions on any fixed closed countour . In this work is a contour chosen to be close to the measurements. If all sensors were located a single contour this sensor or measurement contour would be a natural choice for .
II.B. Step 2: an optimal control method
From the method described in the previous paragraph one can compute a complete set of Cauchy conditions, on and on . It is then possible to employ the optimal control method discussed in [ACL.B.Faugeras.12.2] which is summed up below.
The poloidal flux satisfies
In this formulation the annular domain contained between and (see Fig. 1) is unknown since the free plasma boundary is unknown. Moreover the problem is ill-posed as there are two Cauchy conditions on the boundary .
In order to compute the flux in the vacuum and to find the plasma boundary we define a new problem approximating the original one. Let us define a fictitious boundary fixed inside the plasma. We seek an approximation of the poloidal flux satisfying in the domain contained between the fixed boundaries and . The problem becomes one formulated on the fixed annular domain . Freezing the domain to by introducing the fictitious boundary enables to remove the nonlinearity of the problem. The plasma boundary can still be computed as an iso-flux line and thus is an output of our computations.
Let us insist here on the fact that this problem is an approximation to the original one since in the domain between and , should satisfy the Grad-Shafranov equation. The relevance of this approximating model is consolidated by the Cauchy-Kowalewska theorem [Courant:1962]. For smooth enough the function can be extended in the sense of in a neighborhood of inside the plasma. Hence the problem formulated on a fixed domain with a fictitious boundary not "too far" from is an approximation of the free boundary problem. If were identical with then by the virtual shell principle [Shafranov:1972] the quantity would represent the surface current density (up to a factor ) on for which the magnetic field created outside the plasma by the current sheet is identical to the field created by the real current density spread throughout the plasma.
The shape of the internal contour is chosen to be a circle or an ellipse. Although no thorough sensitivity analysis has been conducted concerning the shape and size of the numerical results presented in the next section did not prove to be very dependent on these factors.
No boundary condition is known on and boundary conditions are over determined on . One way to deal with this and to solve such a problem is to formulate it as an optimal control one. We are going to compute a function such that the Dirichlet boundary condition on is such that the Cauchy conditions on are satisfied as nearly as possible in the sense of the Kohn-Vogelius error functional defined below which enables to use Dirichlet and Neumann conditions in a symmetric way.
We consider two separate well-posed sub-problems. In the first one we retain the Dirichlet boundary condition on only, assume is given on and seek the solution of the boundary value problem:
This solution can be decomposed in a part linearly depending on and a part depending on only. We have the following decomposition:
In the second sub-problem we retain the Neumann boundary condition only and look for satisfying the boundary value problem:
in which can be decomposed in a part linearly depending on and a part depending on only. We have the following decomposition:
For given Cauchy data and we look for satisfying where is the regularized error functional defined by
In this expression the first term measures a misfit between the Dirichlet solution and the Neumann solution whereas the second one is a regularization term in which is a regularization parameter.
Let and define the bilinear symmetric forms
and the linear form
as well as the quantity depending on the Cauchy data only
Following the analysis provided in [ACL.B.Faugeras.12.2] the functional can be rewritten as
and its minimum is given by the first order optimality condition or Euler equation which is linear and reads
Solving this equation gives the optimal Dirichlet conditions on from which one can compute or and then look for the iso-contour defining the plasma boundary.
Iii. Numerical Implementation and Results
III.A. The toroidal harmonics expansion
The Legendre functions of degree 1 and half integer order involved in the toroidal harmonic expansion are computed thanks to the Fortran code provided together with reference [Segura:2000].
In the first step of the method one also needs to provide a pole for the toroidal coordinate system and the number of harmonics to be used.
Because the internal toroidal harmonics functions are singular at the pole of the coordinate system, the expansion procedure provides meaningful results if the pole lies inside the unknown plasma region and not too close from the boundary. An excellent candidate for the choice of this pole is the current center defined from the moments of the plasma current density [Zakharov:1973, Braams:1991].
where is any closed contour containing the plasma and the domain it defines. These quantities can be directly approximated as weighted sums of magnetic measurements and then precisely recomputed as integrals on the contour at every point of which the flux and the field can be evaluated.
The number of harmonics to be used need not be large to obtain precise fit to the measurements [ACL.B.Faugeras.14.2]. Typically the order of the harmonics is chosen to be , or . Figure 2 shows a WEST example case of boundary reconstruction with the code using toroidal harmonics only (called VacTH). Two reconstructions are performed using harmonics of order 4 (case 1) and 6 (case 2). In both cases the root mean square errors for the fit to the measurements is of for B sensors and for flux loops. The Cauchy conditions computed in both cases differ from only small values: and . However the plasma boundary reconstructed with VacTH differs quite a lot from one case to the other. In case 1 it is accurate whereas in case 2 it presents unrealistic oscillations. This shows two things. On the one hand the interpolation of the magnetic measurements using toroidal harmonics and the evaluation of Cauchy data on the contour is accurate regardless of the order of the harmonics used. On the other hand the reconstructed plasma boundary is quite dependent on the order of the harmonics used and particularly on the order of the internal harmonics. For all WEST cases that have been treated up today the choice of order 4 harmonics allways gave accurate results. However it can be understood from this example that the choice of this order might not be easy or even possible for other Tokamak configurations. Figure 3 shows such an issue. We used VacTH to reconstruct the plasma boundary using ASDEX UpGrade experimental data for shot 25374 from the ITM-WPCD database [ITM]. The same choice for the order of the toroidal harmonics gives a good reconstruction at time s but this is not the case at time s. This illustrates the need for the optimal control method proposed in this paper and implemented in a code called VacTH-KV.
III.B. The optimal control problem
The resolution of the optimal control problem in VacTH-KV is based on a classical finite element method [Ciarlet:1980]. Although over choices would have been possible the advantage of this numerical method over finite differences or boundary integral methods is its simplicity when it comes to adapt the mesh to any internal and external contours, and to extract the plasma boundary from the finite dimensional representation of the flux.
Given the fixed domain let us consider the family of triangulation of , and the finite dimensional subspace of defined by
Let us also introduce the finite element space on
Consider a basis of and assume that the first mesh nodes (and basis functions) correspond to the ones situated on . A function is decomposed as and its trace on as .
Given boundary conditions on and , on one can compute the approximations , , and with the finite element method. These computations involve two different linear systems. One for the Dirichlet-Dirichlet problem and one for the Dirichlet-Neumann problem with two symmetric positive definite matrices and .
In order to minimize the discrete regularized error functional, we have to solve the discrete optimality condition which reads
which is equivalent to look for the vector solution to the linear system
where the matrix represents the bilinear form and is defined by
and is the vector .
The matrices which constitues are evaluated by
where is the trivial extension which coincides with on and vanishes elsewhere.
In the same way the right hand side is evaluated by
It should be noticed here that matrices , and depend on the geometry of the problem only. Therefore an internal contour being given they can be computed once as well as their Cholesky decomposition and be used for the resolution of successive problems with varying input data. Only the righthand side has to be recomputed.
An important point however is that in the same way that the pole of the toroidal coordinate system moves following the plasma current center, the inner contour might also need to be displaced. We have adopted the following strategy consisting in defining a number of points in the vacuum vessel. Each of them is the center of a circle defining a contour . A different mesh is precomputed for each of these contours as well as the associated finite elements matrices (see Figure 4 for an example). Each time a new set of magnetic measurements is provided the plasma current center is computed in the first step of the algorithm involving the toroidal harmonics expansion. Then one looks for the internal contour whose center is the closest to this point and use the associated and precomputed mesh and matrices in the second step in order to minimize the Kohn-Vogelius functional. This enables fast computations. The computing time is of the order of for one boundary reconstruction in the WEST configuration on a Laptop with two quadcore processors of frequency 2.40 Ghz.
Figure 5 shows the reconstructed plasma boundaries for the same WEST example case as in figure 2 but this time with VacTH-KV. Using toroidal harmonics of order either 4 or 6 for the evaluation of Cauchy conditions on leads to the same reconstruction almost superimposed with the reference boundary from the equilibrium code CEDRES++. In addition to the unrealistic results given by VacTH for the two AUG cases with real experimental measurements, figure 3 also shows the excellent results provided by VacTH-KV on these same cases. As a last example illustrating the non tokamak dependance of the method figure 6 shows the plasma boundary reconstructed by VacTH-KV on a JET case (shot 74221 at time 46s). In this case VacTH alone was not able to reconstruct a proper plasma boundary.
In all our numerical experiments the regularization parameter is given the value of as suggested by the L-curve method analysis [Hansen:1998] conducted in [ACL.B.Faugeras.12.2].
Finally we would like to point out the fact that the optimal control method presented in section II.B. is completely unchanged if poloidal field coils or passive structures with measured currents are included in the domain . In fact this is the case for the WEST, AUG and JET numerical results presented here.
Iv. Summary and Conclusion
In this paper a method for plasma boundary reconstruction is proposed. It is decomposed in two steps. A toroidal harmonics expansion is first used to interpolate discrete magnetic measurements on an external contour. This interpolation using functions which are exact analytic solutions of the equation satisfied by the flux is accurate. On the other hand we have shown that the extrapolation of the flux towards the plasma boundary using toroidal harmonics only can in some cases be inaccurate. This is the reason why we introduce a second step in which a Cauchy problem on a fixed domain is solved by an optimal control method. These two steps together provide an accurate and robust plasma boundary reconstruction method. It is generic and can be used for any Tokamak having ferromagnetic structures or not. An implementation using a finite elements discretization is proposed. The precomputation of several meshes and all associated linear algebra quantities allows fast computation. A code called VacTH-KV has been developped and is available on the ITM-WPCD platform. This paper does not deal with the computation of error bars on the plasma boundary reconstruction results. This is of course of interest as for any ill-posed inverse problem and might be conducted in the future for example using the epsilon-nets technique [Zaitsev:2014]
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 (ITM-WPCD). The views and opinions expressed herein do not necessarily reflect those of the European Commission.
The author would like to thank his colleagues at the university of Nice, Jacques Blum, Cedric Boulbe and Holger Heumann for many useful discussions. The input magnetic measurements for WEST were provided by colleagues from the IRFM-CEA at Cadarache, Eric Nardon and Remy Nouailletas, and for AUG and JET by the ITM-WPCD project.
Finally the author would like to thank the two anonymous reviewers for their constructive criticism from which the paper has benefitted significantly.