Kohn–Vogelius formulation for plasma geometry identification problem

We study the problem of plasma geometry control problem in a tokamak. The domain location and shape are determined using an approach based on the Kohn–Vogelius formulation and topological asymptotic method. We present a one-shot numerical procedure based on the developed asymptotic formula and use it on different test configurations.

We remark that since p is unknown, V is also unknown. This problem has been considered in different works, in which control or parametric optimization methods are used [6][7][8]. In our previous work [9], we used the topological sensitivity method to find the unknown plasma boundary p . We present in this work a new formulation of the inverse problem based on the Kohn-Vogelius method. It leads to define, for any given plasma domain D ⊂ , two forward problems.
• A first problem, associated with the Neumann boundary condition : (3) • A second one, associated with the measured velocity ψ m : Note that if ∂D and the plasma boundary p coincides, then the misfit between the solutions vanishes, and then ψ n D = ψ d D . Using this remark, the studied inverse problem can be solved by minimizing the following energy functional [10][11][12]: To solve this optimization problem and to reconstruct the unknown plasma domain p , in this paper, we propose a fast and accurate topological optimization algorithm based on the topological sensitivity analysis concept. An outline of the paper is as follows: -The main ideas of the proposed geometry reconstruction approach are illustrated in the next section. -In Sect. 3, we derive a topological asymptotic expansion describing the variation of the Kohn-Vogelius functional K with respect to the creation of a small geometric perturbation inside the domain . -The implementation of the proposed numerical optimization algorithm is done in Sect. 4. We present some identification results and discuss the influence of some parameters (the location, size, and shape) on the accuracy of the reconstructed plasma domain. -Finally, a conclusion is presented in Sect. 5.

Proposed reconstruction method
Classical shape optimization methods are based on the shape derivative concept [13][14][15]. The optimization process is limited to determine the optimal location of the initial domain boundary. Its main drawback is that is does not allow any topology changes. The initial and optimal domains have the same topology. To overcome this difficulty, the homogenization theory [16][17][18][19] represents important developments in this direction. One of its main features consists in considering the domain as a composite material with density ranging from zero to one, rather than defining domains via their characteristic functions, i.e., by a discrete 0-1 valued function. Rather than optimizing the shape by deforming the boundary, the problem was relaxed by considering composite materials defined by a material density distribution and a microstructure. However, the range of application of this approach is mainly restricted to linear elasticity and particular objective functions. For these reasons, global optimization techniques like genetic algorithms and simulated annealing have been proposed (see, e.g., [20]), but these methods have a high computational cost and cannot be applied to industrial problems. To improve the optimization process, the level set method has been used in the field of shape and topology optimization [21,22]. This approach has been initially introduced by Osher and Sethin [23] for numerically tracking fronts and free boundaries. This method can handle some topology changes. Indeed, the level set can easily remove holes but cannot create new holes in the middle of a shape. In practice, this affect can checked by varying the initialization, which yields different optimal shapes with different topologies.
During the last two decades, the topological sensitivity notion [24][25][26][27][28] has received much attention and has been successfully used for solving various shape and topology optimization problems in different fields. The topological sensitivity-based methods are concerned with the variation of a cost function with respect to a topology modification of a domain. The topology of the computational domain can change during the optimization process by creation of holes. The most simple way of modifying the topology consists in creating a small hole in the domain. In the case of structural shape optimization, creating a hole means simply removing some material. In the case of fluid dynamics, where the domain represents the fluid, creating a hole means inserting a small obstacle. The situation is similar in electromagnetism.
In this paper, we use the topological sensitivity concept, and we built a fast and accurate topological optimization algorithm for solving the minimization problem (5) and identifying the unknown plasma domain p from boundary measurement of the potential field. The idea is to study the topological sensitivity of the function K regarding a small topological domain perturbation ω z,ε of the domain . It leads to an asymptotic expansion of the form where g(ε) is a positive scalar function going to zero with ε. This formula is called the topological asymptotic expansion, and δK is called the topological gradient. Using this expansion, the minimum of K corresponds to the location in where δK is the most negative. In fact, if δK(z) < 0, then K( \ω z,ε ) < K( ) for small ε.
In general, this topology optimization approach leads to build a sequence of geometries ( n ) n≥0 , with 0 = . At the nth iteration the geometry n+1 is obtained by creating a new geometric perturbation ω n in the domain n ; n+1 = n \ω n . The perturbation ω n is constructed as where δK n is the associated topological gradient (computed in n ), and c n is a negative constant chosen so that K( n \ω n ) -K( n ) decreases as most as possible.
• The stopping criteria is given by the natural optimal condition This condition coincides with that obtained by Buttazzo and Dal Maso [29] for the Laplace equation using the homogenization theory.
• This topology optimization algorithm can be viewed as a descent method, where the descent direction is determined by the sensitivity function δK n (the topological gradient), and the step length is given by the volume variation |ω n | = | n \ n+1 |. • The determination of the optimal step length |ω n | (i.e., determination of the optimal value of the constant c n ) is still an open question. In practice, this constant is adjusted using the binary search method. In this paper, we extend this approach for solving the problem of the plasma magnetic equilibrium in a tokamak. Since the real-time identification of the plasma domain is a key point to access high-performance regimes, we propose in this study a fast and accurate reconstruction algorithm (one-shot algorithm) based on the topological sensitivity analysis method.

Topological sensitivity analysis
Using the axisymmetric configuration and an horizontal cut, we can rewrite the used formulation as follows: find the unknown domain p occupied by the plasma as the optimal solution to the optimization problem where ψ n D and ψ d D are solutions to the following systems: To solve this problem, we derive a topological sensitivity analysis for the Khon-Vogeliustype functional K. It consists in studying the variation of K regarding the presence of a small plasma domain ω z,ε around the point z ∈ with Dirichlet boundary condition on ∂ω z,ε .
Using the above notation, we define the function K( \ω z,ε ) by where ψ n ε is the solution to the perturbed Neumann problem The function ψ d ε is the solution to the perturbed Dirichlet problem We can derive the following asymptotic expansion for the shape function K (see [30]): Let ω z,ε be a small topological perturbation in of the form ω z,ε = z + εω ⊂ . Then the function K admits the following asymptotic expansion: where ψ n 0 and ψ d 0 are solutions to the following systems:

Numerical investigations 4.1 Reconstruction algorithm
In this section, we present a fast and simple one-shot numerical reconstruction algorithm. Our numerical procedure is based on the asymptotic formula presented in Theorem 3.1, which describes the variation of the Kohn-Vogelius-type functional K with respect to the creation of a small hole ω z,ε inside the domain . From Theorem 3.1 it follows that the function K satisfies the following topological asymptotic expansion: where δK is the topological gradient defined as with solutions ψ n 0 and ψ d 0 to the Neumann and Dirichlet problems (10)- (11). According to the main idea of the topological sensitivity analysis method, the unknown plasma domain p is likely to be located at the zone where the topological gradient δK is negative. To present our proposed procedure, we introduce some notations. Let δ min be the most negative value of the topological gradient δK in (i.e δ min = min x∈ δK(x)). For all ρ ∈ [0, 1], we define the zone The aim is reconstructing the location and shape of the unknown plasma domain p ⊂ from a given measurement of the potential field on the boundary = ∂ . The main steps of our numerical reconstruction procedure are summarized in the following "one-shot" algorithm.
In this one-iteration algorithm: • The location of the unknown plasma is given by the point z ∈ where the topological gradient δK is most negative (i.e., z = argmin x∈ δK(x)). In other words, the topological gradient tell us that the plasma domain is located around the point z ∈ . • The size of the domain occupied by the plasma is adjusted via the choice of the parameter ρ ∈ (0, 1). Its boundary is approximated by a level-set curve of the topological gradient δK (i.e., the plasma domain is delimited by a well-chosen isovalue of δK) • The computation of the plasma domain size (in Step 3) can be considered as a descent method where the step length is given by the volume variation |V k | = |ω ρ k \ω ρ k+1 |. The determination of ρ can be viewed as a line search step. However, determining the optimal value of the parameter ρ and optimizing the size of the domain ω solution to min ρ∈(0,1) is still an open question. To speed up the convergence of our optimization algorithm, we perform a binary search approach (dichotomy method) for approaching as most as possible the best value of ρ starting from ρ = 1 2 .
Remark 4.1 In the particular case where the exact plasma domain ex p is known, the best value ρ * of the parameter ρ can be determined as the minimum of the following error functional: where meas(B) is the Lebesgue measure of a set B ⊂ R 2 .

Numerical implementation
In our numerical implementation the measurements data ψ m are synthetic, that is, generated by numerical simulations. More precisely, let ex p be the exact plasma domain to be identified. The Dirichlet measurement data ψ m on the boundary are retrieved as ψ m = ψ n ex , where ψ n ex is the solution to the Neumann problem (see (6)) in the domain \ ex Here the exact plasma domain is used only for generating the measured data ψ m . In the numerical simulation the vacuum vessel region is defined by the disc with center x 0 = (2, 0) and radius r = 1, i.e., = B(x 0 , 1). The resolution of the boundary value problems (P n 0 ) and (P d 0 ) is based on a classical P 1 finite element method. The approximated solutions are computed using a uniform mesh on the boundary ∂ (defined by the circle C(x 0 , 1)) with step size h = π/300 ≈ 0.01. The numerical procedure is implemented using the free software FreeFem++ (see http://www.freefem.org/ff++/).
Next, we present some reconstruction results showing the efficiency of the proposed one-shot algorithm.

Numerical simulations
Here we apply our one-shot numerical algorithm for reconstructing the plasma domain in different situations. We start our numerical investigation by the identification of an unknown elliptical plasma domain from boundary measured data (see Sect. 4.3.1). After that, we will discuss the influence of some numerical parameters on the accuracy of the reconstructed results: the effect of the location in Sect. 4.3.2, the effect of the size in Sect. 4.3.3, and the effect of the shape in Sect. 4.3.4.

Reconstruction results
In this section, we apply our one-shot algorithm for reconstructing an elliptical plasma domain from boundary measurement. The measured data ψ m is generated using the ellipse ex p = B(x 0 , r 1 , r 2 ) with center x 0 = (2, 0) and radii r 1 = 0.4 and r 2 = 0.5. The obtained numerical results are described in Figs. 1 and 2: -In Fig. 1, we plot the isovalues (level-set curves) of the topological gradient δK. We can note that the zone where the topological gradient g is negative (the red region) nearly coincides with the exact plasma domain ex p (see Fig. 2). From this figure we can observe that the most negative value of δK is given by δ min = -0.0092, which is reached at the point z = (1.9, 0).
-To compute the optimal value ρ of the parameter ρ and estimate the size of the reconstructed plasma domain, we study the variation of the error function, which allows us to deduce that the optimal value of ρ is given by ρ = 0.12, i.e., 1ρ = 0.88. According to our algorithm, the reconstructed plasma domain is defined as p = x ∈ ; δK(x) ≤ 1ρ δ min , with ρ = 0.12 and δ min = -0.0092.  Fig. 2, we show the exact (delimited by the red line) and reconstructed (delimited by the black line) plasma domains. We conclude here that our numerical algorithm provides an efficient reconstruction result in only one iteration.

Effect of the location
In this section, we discuss the effect of the plasma location on the reconstruction result. The aim is reconstructing a circular plasma ex p = z + 0.2B(0, 1) located at the point z ∈ . We present in Fig. 3 the isovalues of the topological gradient δK for different plasma locations z i , 1 ≤ i ≤ 4. As we can observe here, the topological gradient detects the exact plasma location, i.e., the most negative region coincides with the exact plasma domain.

Effect of the size
In this section, we discuss the effect of the plasma size on the reconstruction result. We illustrate in Fig. 4 the isovalues of the topological gradient δK for different plasma sizes r i , 1 ≤ i ≤ 3.

Effect of the shape
In this section, we discuss the effect of the plasma shape on the reconstruction result. The first test is devoted to an elliptical-shaped plasma. The second test is concerned with a rectangular-shaped plasma. As we can observe in Fig. 5, the topological gradient gives a good approximation of the unknown plasma in both cases.

Conclusion
This paper is concerned with a topology optimization problem related to the plasma magnetic equilibrium in a tokamak. Since the real-time identification of the plasma domain is a key point to access high-performance regimes, we propose a fast and accurate reconstruction algorithm (one-shot algorithm). Our approach is based on the Kohn-Vogelius formulation and the topological sensitivity analysis method. Compared to the classical geometric reconstruction algorithm, our proposed algorithm has several advantages: • Only one iteration is needed to identify the location of the unknown plasma domain, which significantly reduces the running times. • The shape of the unknown plasma domain is approximated by a level-set curve of a scalar function, which gives it more flexibility and a wide range of applications.
• Unlike the classical geometric reconstruction algorithms such as level-set method [21,22] or homogenization method [16], our numerical approach is not sensitive to the initial geometry. The topological sensitivity analysis method consists in creating some geometric perturbations in the initial domain.