Magnetics only real-time equilibrium reconstruction on ASDEX Upgrade

Real-time reconstruction of the magnetic equilibrium provides fundamental control of plasma shape and position in a tokamak. Details of the implementation of the equilibrium reconstruction code developed for the ASDEX Upgrade tokamak (JANET++) are summarized. Cubic Hermite splines are introduced as current density basis functions for solving the Grad–Shafranov equation. The choice of the optimal Tikhonov regularization parameter is discussed. The code is validated by comparing the results of the equilibrium reconstruction with those of further equilibrium reconstructions available on ASDEX Upgrade (CLISTE and IDE). In a high time resolution study of a discharge with edge localized modes (ELM), the poloidal asymmetry of the fits and magnetic probe measurements suggest that the real-time equilibrium reconstruction captures the essential features of the current density redistribution in an ELMing edge plasma. An efficient algorithm to locate multiple X-points and identify the active one in advanced X-divertor and snowflake divertor configurations is presented.


Introduction
Plasma position and shape control on ASDEX Upgrade is performed using function parameterization [1,2].This method relies on the statistical analysis of a large database of simulated experiments to obtain a functional representation of the intrinsic physical parameters of the plasma position and shape in terms of the values of the magnetic probe and flux loop measurements.In contrast to function Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.parameterization, real-time magnetic equilibrium reconstruction solves the Grad-Shafranov equation using magnetic probe and flux loop measurements as constraints to provide flux surfaces with a well defined plasma boundary.Real-time magnetic equilibrium reconstruction is used for the control of plasmas in a number of tokamaks [3][4][5][6].An equilibrium reconstruction code can also use a Kalman filter based eddy current estimator to improve fast transient accuracy [7,8].The main use of real-time equilibria on ASDEX Upgrade was to point the gyrotron at a rational q surface of the magnetic equilibrium for neoclassical tearing mode (NTM) stabilization [9].Real-time equilibria will also be required to support studies of advanced X-divertor (XD) and snowflake divertor (SF) configurations produced by the new upper divertor coils currently being installed in ASDEX Upgrade [10].The database to extend the operation of function parameterization to these configurations is in progress.Real-time equilibrium reconstruction in these configurations requires an efficient algorithm to locate multiple X-points and identify the active one.
In general, external magnetic measurements alone can only provide global current and pressure profile parameters such as plasma current, poloidal beta and internal plasma inductance [11,12].However, it has been shown that in a toroidal plasma bounded by a magnetic separatrix, additional moments of the toroidal current density profile, strongly localized at the plasma edge, can be identified using external magnetic measurements alone [13].It is known that real-time equilibrium with inputs from magnetic probes and flux loop measurements alone gives a current density profile and safety factor in the plasma centre that are poorly estimated [14].Internal constraints from the motional Stark effect [15,16] or polarimetry [17] are required to accurately determine the current density profile.Equilibrium reconstruction can be improved by adding pressure constraints and the current diffusion equation to perform real-time kinetic equilibrium reconstruction [18][19][20].Recent efforts have been made to use artificial intelligence to improve equilibrium reconstruction and it is concluded that the extent to which this can provide real-time control will depend on how much faster it can be made, particularly in the area of the Bayesian experimental analysis [21].
The real-time equilibrium reconstruction code for the ASDEX Upgrade tokamak [22], has been migrated to a C++17 code (JANET++) [23].It closely follows the offline equilibrium reconstruction algorithm implemented in EFIT [12], with the real-time modifications necessary to enable rtE-FIT to control the plasma on DIII-D and KSTAR [3].The fast loop of the real-time equilibrium reconstruction computes flux surfaces for each solver input and subsequently updates the plasma shape and position control variables.In the slow loop, a new current distribution is generated on the updated flux surfaces and the solver is run only once, but the results are not iterated upon to improve the convergence of the fits to the measurements in the real-time equilibrium reconstruction.
Section 2 summarises the Grad-Shafranov equation describing the magnetic equilibrium of a tokamak and the Shafranov integrals for evaluating the global parameters of the equilibrium.Cubic splines and cubic Hermite splines are introduced as basis functions for solving the Grad-Shafranov equation.In section 3, the variation of the value of the regularization parameter for the weighting of the data fit relative to a smooth fit of the current density profile is studied.In section 4, examples of equilibrium reconstructions from JANET++ are reported and the obtained global values of the plasma parameters are compared with those calculated from offline equilibrium reconstruction.An efficient algorithm for locating multiple X-points and identifying the active one in advanced XD and SF configurations is also presented.In section 5, the magnetic probe and flux loop fits from the equilibrium reconstruction are compared with measurements in a plasma with edge localized modes (ELMs).The movement of the position of the magnetic axis, the location of the lower and upper X-points and the radial position of the inner and outer boundary from JANET++ are compared with those obtained from offline equilibrium reconstructions.In section 6 the current density profiles of L-mode and H-mode discharges are presented and the characteristic current density peak at the edge of an H-mode plasma can be demonstrated.

Grad-Shafranov equation and Shafranov integrals
The magnetic equilibrium of an axisymmetric toroidal plasma in a tokamak is described by the Grad-Shafranov equation [24,25]: )) (1 where ψ is the poloidal magnetic flux enclosed by a magnetic surface, p is the plasma pressure and F = RB ϕ is related to the poloidal current, J pol , by: Current density basis functions are assumed to create an expression for the right hand side of equation ( 1) in terms of n coefficients α i and m coefficients β j for the p ′ and the FF ′ terms of the Grad-Shafranov equation, respectively: where B i (ψ N ) are the current density basis functions, ψ N is the normalized radius defined in terms of, ψ axis , the value of ψ on the axis and, ψ bdy , the value of ψ on the plasma boundary: with ψ N = 0 on the axis and ψ N = 1 on the boundary.For the equilibrium reconstructions presented in section 4, an equal number of α i and β j coefficients in the range of 3-6 were used.
In equation (3), a free parameter is required in the current profile model to allow the freedom of the vertical position of the equilibrium to be adjusted during the calculation of the least squares solution [3].The extra term used in JANET++ with the fit coefficient c o is: with the derivative calculated from the equilibrium in the previous step.
The plasma current, I p , and the current centroid radial and vertical positions, R c and z c , are calculated by the following contour integrals on the last closed flux surface: The Shafranov integrals, S 1 , S 2 and S 3 are defined as surface integrals on the last closed flux surface [26], which can be simplified to a contour integral on the last closed flux surface, C: where V is the plasma volume, S pl is the plasma surface, Bpa = µ o I p /s is the mean poloidal field around the contour on the last closed flux surface, s is the plasma circumference of the poloidal cross-section and n is the vector normal to the plasma surface.These integrals are used to calculate the magnetohydrodynamic (MHD) plasma beta, β MHD [14,27]: The normalized plasma inductance, l i , is defined as: where Γ pl is the plasma volume.The plasma energy, W MHD , is defined as: The pressure profile on the normalized radius, p(ψ N ), is found by integrating: This leads to an expression for p(ψ N ) in terms of the fit coefficients and the integral of the basis functions: The sum of the thermal and fast-ion pressure is mapped onto the flux surfaces and the pressure difference in equation ( 12) can be used as a constraint on the coefficients [14].The inclusion of real-time pressure profiles in JANET++ is planned and will improve the determination of the current density profile and the safety factor in the plasma centre.For magnetics only real-time equilibrium reconstruction, the pressure at the separatrix, p(1), is set to zero and W MHD is calculated by substituting p(ψ N ) into equation (10).Cubic Hermite splines with zero derivatives at the knots and cubic splines can be chosen as current density basis functions in JANET++ for equilibrium reconstruction.Cubic Hermite splines have the advantage over cubic splines that they are always positive along the normalized radius and that there is a simple expression for real-time evaluation.For the region between the 3 knots ( K i−1 , K i , K i+1 ) in normalized flux, ψ N , the cubic Hermite spline, H(z), can be expressed as: The position of the knots can be chosen to suit the physical demands of the experiment.In particular, knots can be positioned at the plasma boundary in plasmas where peak current densities are expected at the plasma boundary.The cubic Hermite and cubic spline current density basis functions for knots at a normalized radius of 0., 0.25, 0.5, 0.75, 0.9, 1.0 are shown in figure 1.The cubic Hermite splines with these knot positions were chosen for the equilibrium reconstruction results presented.When comparing the global quantities derived from the equilibrium reconstruction, a satisfactory level of agreement is observed irrespective of the chosen current density basis function.Cubic Hermite splines with gradients at the knot positions on a unit interval with value y j at t = 0 and value y j+1 at t = 1 can be written as : (14) with starting tangent m j at t = 0 and ending tangent m j+1 at t = 1.The choice of a basis function with slopes at the knot positions is more general than that chosen for the equilibrium reconstruction in this study.However, the introduction of non-zero slopes in the basis functions would introduce two additional fitting coefficients for each interval between knots.Consequently, further regularization constraints would be required, as the number of fitting coefficients and regularization constraints should be approximately the same to obtain results from the least squares fitting discussed in the following section.It is acknowledged that only current density functions with zero derivatives at knot positions can be modeled by the basis functions used in the equilibrium reconstruction.Nevertheless, the results presented indicate that the cubic Hermite spline basis functions with zero slope at the knots used in the equilibrium reconstruction adequately reproduce the results of the offline equilibrium codes with sufficient accuracy for real-time control.

Regularization
In equilibrium reconstruction, a linear least squares problem of the form: is considered, where R is the response matrix, b is the set of magnetic probe and flux loop measurements and x are the coefficients of the current density basis functions [3,12].This minimization problem is commonly replaced by a penalized leastsquares problem known as Tikhonov regularization [28]: where L is the regularization matrix.The value of the regularization parameter, λ 2 , is important for the quality of the equilibrium reconstruction.Too large a value will result in an oversmoothed solution that lacks details that the desired solution may have, while too small a value will result in a calculated solution that over-fits the measured data and unphysical oscillations of the current density may occur.This regularization improves the conditioning of the problem and allows a direct numerical solution.To enforce smoothness, the second derivatives are set to zero.A flat current profile is favoured by setting the first derivative to zero.In JANET++, a weak regularization is applied to the first and second derivatives of the current density profile, linking each set of three adjacent basis function coefficients, x i .The constraint on the first and second derivatives for equally spaced basis functions can be expressed as [16]: and this is applied separately to the p ′ and FF ′ terms of the Grad-Shafranov equation.The constraints on the first and second derivatives must take into account that the knots are unequally spaced.In CLISTE, L is the identity matrix with additional individual penalty terms of the second derivative at each of the knots to account for the case that bootstrap currents at the plasma edge may be present [29].
Magnetics only equilibrium reconstruction with 3-6 basis functions for both the p ′ and FF ′ terms of the Grad-Shafranov equation have been performed for different values of the Tikhonov regularization parameter.The regularization parameter for the p ′ and FF ′ terms can be chosen separately, but in the following they are chosen to be the same.The coefficients of the current density basis functions that gave the smoothest and flattest current density profiles and fitted the 40 probe and 18 flux loop measurements were calculated.A sampling rate of 10 kHz was used for the measurements and the equilibria are calculated at 2 kHz, with the input data averaged over 5 samples.
The JANET++ real-time equilibrium results are compared with the offline equilibrium calculated by CLISTE (magnetics only) [30] and IDE (pressure constrained and current diffusion coupled) [31].The regularization parameter was varied for a discharge (#37 642) with poloidal beta between 0.7 and 1.8 in a plasma with 1 MA plasma current, up to 19 MW neutral beam injection heating and 4.5 MW electron cyclotron heating.As shown in figure 2, l i was found to decrease with increasing regularization parameter.The flattening of the current density profile is a direct consequence of increasing the regularization parameter.In figure 3, the IDE and CLISTE values of l i , β p and W MHD are compared with the JANET++ values as a function of the regularization parameter.There is good agreement for The time dependence of the plasma inductance, l i , the poloidal beta, βp, and the plasma energy, W MHD , calculated by JANET++ for different values of the regularization parameter (0.004 (red), 0.008 (green) and 0.012 (blue)) are compared with the evaluations from IDE (magenta) and CLISTE (black).a regularization parameter of 0.008.This value of the regularization parameter was chosen for the equilibrium reconstruction results presented in the following sections.Numerous discharges have been evaluated with this value of the regularization parameter and the agreement of the JANET++, IDE and CLISTE values of l i , β p and W MHD shown in figure 2 has been repeatedly confirmed.

Equilibrium reconstruction
Real-time equilibrium reconstruction is primarily used for NTM stabilization on ASDEX Upgrade.On DIII-D and KSTAR, magnetics only real-time equilibrium reconstruction is performed by rtEFIT for position and shape control [3,32].In this section, a comparison of the JANET++ magnetic equilibria with those calculated by CLISTE and IDE is made to validate that JANET++ can reproduce the offline equilibrium global values such as plasma volume, l i , β p and W MHD .The results presented here are not from the control system version of JANET++, but from the development version running offline with stored data as input for the equilibrium reconstruction.The grid size of the equilibrium reconstruction in JANET++ is 33 in the radial direction and 65 in the vertical direction.The grid size of the equilibrium reconstruction in CLISTE and IDE is 65 in the radial direction and 129 in the vertical direction.
The time evolution of plasma heating, plasma energy and radiated power (#41 027) is shown in figure 4.There is good agreement between the plasma energy calculated from volume integrals of the plasma pressure, W MHD , and that calculated by JANET++ [14,27].The calculated values of β p are also in good agreement.The predicted value of the diamagnetic flux from IDE and JANET++ (diaflfit) is in good agreement with the 8 measurements of the diamagnetic flux.Despite the reduced grid size, JANET++ is able to reproduce the global values of the equilibrium reconstruction.
ASDEX Upgrade is preparing a hardware modification of its upper divertor to investigate alternative divertor configurations, such as the XD and the SF, which have been discussed as possible solutions to the power exhaust problem in a tokamak [10].These configurations can be produced by a pair of invessel coils, each with 4 windings and opposite currents.The identification of the position of several X-points in these equilibria will be required, and an intelligent solution for real-time operation will be needed.An algorithm based on the theory of 'Morse-Smale' complexes [33] has been implemented.If you go 'uphill' from an X-point in the direction of the gradient, you will either reach an extremum (or the grid boundary), hence this gradient path is called a 'stable manifold' in Morse-Smale language.Starting from all identified X-points and moving in all eight spatial directions (horizontal, vertical and diagonal), the next step in the path is found by selecting the direction with the highest value of ψ at its grid point.The extremum reached is recorded (or nowhere, if the result is a grid edge or a closed path), and if a given X-point finds its way to the O-point on the magnetic axis in one of the starting directions, it is a candidate to be an active X-point.In principle, this is a search of dimension (N: number of grid points in one dimension) * number of candidate X-points * 8.It turns out that the paths around In the JANET++ results, the value of the regularization parameter was varied (0.004 (red), 0.008 (green) and 0.012 (blue)).
inactive (secondary) X-points are mostly very short (they do not reach the O-point on the magnetic axis by far) and therefore do not generate much overhead.
An example of the application of this algorithm to the flux surfaces of a discharge with simulated currents flowing in the pair of upper divertor in-vessel coils is shown in figure 5.The active X-point is the one with a ψ value closest to ψ axis .This algorithm has successfully identified the active X-point on 61 test equilibria.This will enable JANET++ to provide realtime equilibrium X-point and strike point information in the XD and SF experiments starting in ASDEX Upgrade next year.

Probe and flux loop measurements and fits
In this section, a comparison of the JANET++ magnetic equilibria with those of CLISTE and IDE is made to validate that JANET++ can reproduce the magnetic axis, the lower and upper X-points, as well as the radial positions of the inner and outer boundaries.A discharge (#36 165) with 800 kA at 2.5 T and increasing gas puffing rate in the quasi-continuous exhaust regime is considered [34].Magnetic probe measurements near the strike line at toroidal positions separated by 180 • (blue and magenta) and the equilibrium reconstruction fit (red) at 2.2, 4, and 6 s are shown in figure 6.Note that the ELM amplitude, as measured by the tile current near the strike line, decreases with increasing gas puffing rate.
The magnetic probe measurements and fits at 2.2 s are shown in figure 7. The maximum perturbation is poloidally located in the regions of X-point formation where the movement of the flux surfaces in response to the ELMs is greatest.The poloidal position of the magnetic probes is shown in figure 8. Interestingly, the measured perturbations exhibit no toroidal amplitude dependence.However, it is challenging to reconcile the excursions in the magnetic probe measurements at all poloidal locations simultaneously.For example, the largest measured excursion is 10 mT while the fit is only 5 mT.The equilibrium reconstruction assumes that the current densities on the flux surfaces are expressed in terms of basis functions for the p ′ and FF ′ terms of the Grad-Shafranov equation.During the ELM, this assumption may not hold, so that the fits cannot be expected to reproduce the measurements exactly.Nevertheless, it is interesting to note that the poloidal asymmetry in the fits and the magnetic probe measurements indicate that the real-time equilibrium reconstruction captures the essential features of the current density redistribution in an ELMing edge plasma.
The position of the lower X-point, the upper X-point, the inner plasma boundary and the outer plasma boundary change with time in response to an ELM.These movements are shown in figures 9-11 respectively.The corresponding changes in the plasma current and loop voltage are shown in figure 12.
For the ELMs at 2.2 s, the separatrix on the inboard side jumps 1 cm outwards and the separatrix on the outboard side jumps 1 cm inwards.The vertical position of the X-point in the lower single null configuration jumps downwards by 2 cm.A detailed study of the positional changes of the plasma strike points before and after ELMs in JET has been carried out [35].It was found that in a sufficiently diamagnetic edge plasma, the inner edge toroidal current density can be negative, since the FF ′ contribution to the toroidal current density is enhanced by 1/R, while p ′ is multiplied by R.

Current density profiles
The current density profiles of a plasma in the L-mode and H-mode are shown in figure 13.The peeling-ballooning modes, which are related to the large edge parallel current and the steep pressure gradient in the pedestal region, are widely accepted as the mechanism of type-I ELMs, since the stability limit of the pedestal predicted with the peeling-ballooning mode is in agreement with the experiment [36][37][38][39].The current density peak at the edge of the H-mode plasma is clearly seen in figure 13.Adding a pressure constraint to the equilibrium reconstruction, shows that the current density peak inside the separatrix can be underestimated if only magnetic measurements are used in the equilibrium reconstruction [38].In this discharge, soft x-ray tomography shows a poloidal asymmetry in the radiated power density on the q = 1 surface, similar to that in figure 13 of [40].Possibly tungsten accumulation in the magnetic island on the q = 1 surface is the cause of this asymmetry.It is known that the contours of constant pressure in a rotating plasma with high Mach number no longer lie on flux surfaces [41].The inclusion of inertia due to toroidal plasma rotation is an extension of the Grad-Shafranov equation and is not relevant for ASDEX Upgrade discharges with toroidal Mach numbers for deuterium less than 0.2 [31].Since the measured W concentration was of the order of 10 −3 [40], an asymmetry in the tungsten concentration on the flux surfaces will not produce a relevant pressure asymmetry.Poloidal temperature perturbations on a flux surface due to a poloidal asymmetry in the radiated power density will not occur because the parallel heat transport ensures that the electron temperature on a flux surface is constant.Therefore, for this discharge on ASDEX Upgrade, the current density can still be expected to be described by equation (3).

Conclusions
Cubic Hermite splines have been introduced as basis functions for the current density in the real-time equilibrium code, JANET++.A study of the dependence of the Tikhonov regularization parameter as a function of the poloidal beta showed that the plasma inductance, l i , decreases with increasing regularization parameter.The best agreement between the values of β p and W MHD calculated by IDE and CLISTE and by JANET++ are found for a regularization parameter value of 0.008.The current density peak at the edge of the Hmode plasma is clearly seen in the equilibrium reconstruction with magnetic probe and flux loop measurements as the only inputs.
The poloidal signature of the magnetic perturbations due to ELMs was studied.Temporal changes in the inner and outer radius of the separatrix and the upper and lower X-point positions result from the ELM-induced changes in the current density profile.A poloidal asymmetry in the fits and magnetic probe measurements suggests that the real-time equilibrium reconstruction captures the essential features of the current density redistribution in an ELMing edge plasma.A comparison of the JANET++ magnetic equilibria with those of CLISTE and IDE showed that JANET++ can reproduce the values of the magnetic axis, the lower and upper X-points, and the radial positions of the inner and outer boundaries with the required accuracy.The Morse-Smale algorithm for locating multiple X-points and identifying the active one was successfully tested on 61 equilibria generated with simulated currents flowing in the pair of upper divertor in-vessel coils currently being installed in the ASDEX Upgrade.Therefore, real-time equilibrium reconstruction with JANET++ will be able to handle multiple X-points in the advanced XD and SF configurations produced by the upper divertor coils in future experiments on ASDEX Upgrade.has reviewed and edited the content as needed and assumes full responsibility for the content of this publication.
This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 -EUROfusion).Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission.Neither the European Union nor the European Commission can be held responsible for them.

Figure 2 .
Figure2.The time dependence of the plasma inductance, l i , the poloidal beta, βp, and the plasma energy, W MHD , calculated by JANET++ for different values of the regularization parameter (0.004 (red), 0.008 (green) and 0.012 (blue)) are compared with the evaluations from IDE (magenta) and CLISTE (black).

Figure 3 .
Figure 3.The calculations of l i , βp and W MHD by IDE (magenta) are compared with the evaluations by CLISTE (black) and JANET++.In the JANET++ results, the value of the regularization parameter was varied (0.004 (red), 0.008 (green) and 0.012 (blue)).

Figure 4 .
Figure 4. Time evolution of (a) plasma current, Ip, (b) injected NBI, ECRH and radiated power, (c) line integrated density at the core (H1) and edge (H5) of the plasma, (d) the plasma energy, W MHD , from JANET++, CLISTE and IDE and (e) the beta poloidal, βp, from JANET++, CLISTE and IDE.In (f), the 8 measurements of the diamagnetic flux are compared with the calculated values from JANET++ and IDE.

Figure 5 .
Figure 5. Flux matrix with ψ contours of a simulated discharge with current in the upper divertor coils, showing the locations of X-points (yellow dots) and O-points (red dots).The active X-point is the one with a ψ value closest to ψ axis that is connected to the O-point on the magnetic axis.The paths connecting the X-points to the O-points are indicated by magenta lines.

Figure 6 .
Figure 6.Magnetic probe measurement near the inner strike point at toroidal positions separated by 180 • (Bthe25 + Bp15i25) and fit from the equilibrium reconstruction (red) for different time windows.The outer divertor tile current (14DUAu) in arbitrary units shows the ELM timing.

Figure 7 .
Figure7.The magnetic field perturbation measurements toroidally separated by 180 degrees (blue and sea blue) and fits (green) of the ELMs at 2.2 s.The maximum perturbation is near the X-point, where the movement of the flux surfaces in response to the ELMs is greatest.

Figure 8 .
Figure 8.The poloidal position of the magnetic probes, with the poloidal array plotted in red.The 5 equilibria are separated in time by 2.0 ms and show the movement of the flux surfaces in response to an ELM.

Figure 9 .
Figure 9. Movement of the lower X-point in the radial and vertical direction in response to ELMs.

Figure 10 .
Figure 10.Movement of the upper X-point in the radial and vertical direction in response to ELMs.

Figure 11 .
Figure 11.Movement of the inner and outer radii for the equilibrium construction in response to ELMs.

Figure 12 .
Figure 12.Plasma current and loop voltage spikes in response to ELMs.

Figure 13 .
Figure 13.Current density profiles of a plasma in (a) the L mode and (b) the H mode.