ATEQ: Adaptive Toroidal Equilibrium code

A radially adaptive numerical scheme is developed to solve the Grad-Shafranov equation for axisymmetric magnetohydrodynamic equilibrium. A decomposition with independent solutions is employed in the radial direction and Fourier decomposition is used in the poloidal direction. The independent solutions are then obtained using an adaptive shooting scheme together with the multi-region matching technique in the radial direction. Accordingly, the Adaptive Toroidal Equilibrium (ATEQ) code is constructed for axisymmetric equilibrium studies. The adaptive numerical scheme in the radial direction improves considerably the accuracy of the equilibrium solution. The decomposition with independent solutions effectively reduces the matrix size in solving the magnetohydrodynamic equilibrium problem. The reduction of the matrix size is about an order of magnitude as compared with the conventional radially grid-based numerical schemes. Also, in this ATEQ numerical scheme, no matter how accuracy in the radial direction is imposed, the size of the matrices basically does not change. The small matrix size scheme gives ATEQ more flexibility to address the requirement of the number of Fourier components in the poloidal direction in the tough equilibrium problems. These two unique features, the adaptive shooting and small matrix size, make ATEQ useful to improve tokamak equilibrium solutions.


I. INTRODUCTION
Solving the magnetohydrodynamic (MHD) equilibrium problem is fundamental in plasma physics and plays an essential role, in particular, in the magnetic confinement approach to fusion. In axisymmetric geometries, the equilibrium problem is reduced to solving the Grad-Shafranov equation [1,2]. Since this equation is nonlinear, a numerical solution is necessary in general. In advanced tokamaks [3], two circumstances conspire to make the solution of the Grad-Shafranov equation particularly challenging. First, advanced tokamaks rely on broad current distributions to increase β, the ratio of kinetic to magnetic pressure. This leads to a current profile peak near the edge and to the sensitivity of the stability limit to details in the geometry of the plasma edge. Second, they rely on the H-mode for confinement. The pressure gradients in H-mode drive localized, peaked bootstrap currents near the edge that add to the difficulty in two ways, first by increasing the stiffness of the Grad-Shafranov problem and second by increasing the accuracy needed to calculate the stability of Edge Localized Modes (ELM) [4] as well as Resistive Wall Modes (RWM) [5].
Great efforts have been made previously to develop numerical solvers for the Grad-Shafranov equation. The applications for these numerical solvers are diverse, ranging from the interpretation of experimental observations [6,7] to the design of operation scenarios [8], real-time control of experiments [7,9,10], the analysis of the stability and transport properties of various configurations [11], and the optimization of machine designs [3,12].
The diversity of the applications leads to different requirements regarding properties of the algorithm such as speed, accuracy, stability, and flexibility. These different requirements are partly responsible for the multiplicity of solution strategies. The 1991 review article by Takeda and Tokuda [13] describes early codes including J-Solver [14], VMEC [15], TOQ [16], and others [17]- [20]. Subsequent efforts led to the development of the codes CHEASE [21], CORSICA [22], and EFIT [6]. Refs. [23]- [34] describe further works. As reviewed in [13], the methods for solving the Grad-Shafranov equation are categorized into two types: the Eulerian or "direct," and the Lagrangian or "inverse" numerical schemes. The finite difference, finite element, and Fourier decomposition methods are employed to discretize the equation. In all cases, iteration is used to handle the nonlinearity.
Despite the great successes achieved with the existing codes in various scenarios, challenges remain for solving the equilibrium problem, especially for the cases with high beta, strong shaping, and diverter geometries that give rise to separatrices. The need for adaptive solvers was realized a long time ago. It has, for example, led to the development of the VMEC code for 3D equilibria [15]. Later, the edge equilibrium code (EEC) was developed in order to address the numerical challenges pertaining to the tokamak edge equilibrium problem [35].
In this work, we introduce a new adaptive numerical scheme to solve the Grad-Schafranov equation and describe its implementation in the ATEQ (Adaptive Toroidal EQuilibrium) code for tokamaks. The code uses a decomposition with independent solutions in the radial direction and Fourier decomposition in the poloidal direction. It then obtains the independent solutions with adaptive shooting together with the multi-region matching technique in the radial direction. The adaptive numerical scheme in the radial direction improves considerably the accuracy of the equilibrium solution. The decomposition with independent solutions effectively reduces the matrix size in solving the magnetohydrodynamic equilibrium problem. The adaptive numerical scheme has been successfully used in the linear MHD and kinetic stability codes, AEGIS [36] and AEGIS-K [37].
In addition to its adaptive nature, the reduction of the matrix size by ATEQ is about an order of magnitude, as compared to the conventional radially grid-based numerical schemes.
Also, in this ATEQ numerical scheme, no matter how accuracy in the radial direction is imposed, the size of the matrices basically does not change. Note that all numerical schemes for solving the Grad-Shafranov equation ultimately reduce to solving matrix equations.
The size of matrices then matters. To achieve high accuracy, especially for tough problems related to the axis, X-point, or pedestal, etc. one has to increase the grid density in the radial and poloidal directions in the grid-based codes, or the radial grid density and the number of poloidal Fourier components in the Fourier-decomposition based codes. The dramatic reduction of matrix size by ATEQ is important for this research. The Fourier-decomposition based codes remain to be important tools in this field, for example CORSICA is used for ITER, VMEC is still popular. The small matrix size scheme gives ATEQ more flexibility to address the requirement of the number of Fourier components in the poloidal direction for tough equilibrium problems.
The remainder of this paper is organized as follows: Sec. II introduces the MHD equilibrium equations; Sec. III describes the formulation of numerical equations; Sec. IV gives the numerical procedure and results; Sec. V presents the benchmark studies and comparison with the existing equilibrium codes; Lastly, Sec. VI presents the conclusions and discussion.

II. MHD EQUILIBRIUM EQUATIONS
In this section, we describe the MHD equilibrium equations and the goal of this work.
Force balance, Ampére's law, and the absence of magnetic charge form the basic set of equations describing the MHD equilibrium for a static plasma (V = 0) [39]: where B is the magnetic field, J represents the current density, p denotes the pressure, µ 0 is the magnetic constant, and boldface denotes the vectors.
The paper addresses axisymmetric toroidal equilibria. For such equilibria it is convenient to use a cylindrical coordinate system (X, Z, φ), where φ is the toroidal angle, Z denotes vertical coordinate, and X is radial coordinate from the toroidal axisymmetric axis on the φ = 0 plane. In this coordinate system, the magnetic field in the axisymmetric case can be represented as [39]: where χ is the poloidal magnetic flux. Both pressure p(χ) and g(χ) are flux functions.
where prime denotes the derivative with respect to the poloidal flux χ. The MHD equilibrium is fully determined by χ.
The two free functions p(χ) and g(χ) need to be specified to determine χ from Eq. (5).
In practice, one usually specify p and g as the functions of normalized fluxχ = χ/χ a , where χ a is the poloidal flux at the edge or on the last closed flux surface and the poloidal flux is assumed to be zero at the magnetic axis.
The goal of the present paper is to lay out a new numerical scheme to solve Eq. (5) and describe the ATEQ code that implements this scheme. The paper restricts attention to the fixed boundary problem, i.e., solving Eq. (5) with the plasma boundary specified. We defer consideration of the free boundary problem to future work.

III. FORMULATION OF NUMERICAL EQUATIONS
In this section, we describe the numerical scheme to solve the Grad-Shafranov equation A. Decomposition of the Grad-Shafranov equation In this subsection, we introduce the radial, poloidal, and toroidal coordinates and project the Grad-Shafranov equation onto this coordinate system. We then use Fourier decomposition to decompose the equations.
To solve the Grad-Shafranov equation, we introduce the coordinate system (ψ, θ, φ), with ψ labelling the radial grids and θ being the poloidal angle. The coordinates ψ and θ is general, only requiring that the Jacobian J = 1 ∇ψ × ∇θ · ∇φ remains finite. In this coordinate system one can obtain where Therefore, one has Here, we have denoted ∂ ∂θ = iM with M being the matrix specifying the poloidal Fourier numbers, since the Fourier decomposition with θ will be introduced later on. Using this decomposition the Grad-Schafranov equation (5) can be reduced to the set of first order differential equations: where To solve the set of equilibrium equations, Eqs. (7) and (8), the following Fourier decompositions are introduced, Here, M max represents the maximum Fourier component to be used. Introducing the Fourier decomposition in Eq. (9) the set of equations (7) and (8) becomes the set of matrix equations with the coefficients becoming the maxtrices as defined as follows Note that for the non-up-down symmetric system the Fourier components are complex. The set of matrix equations in complex can be written as Here, χ and A 1 are the vectors in the Fourier space with the total components M = 2m max +1 for each and F ij are the matrices with dimension M × M . Therefore, Eq. (10) represents a set of 2M differential equations. The matrix equation, Eq. (10), can be rewritten concisely as follows where the source term s is usually a nonlinear function of u.

B. Computation of the metric parameters
In this subsection, we describe how the matrix F is computed in the ATEQ code. This is related to the determination of the metric parameters, such as |∇ψ| 2 , ∇ψ · ∇θ, etc.
As in the PEST code [40], we introduce the polar coordinates to compute the metric parameters: where x = X − X 0 and z = Z with X 0 being the major radius at the magnetic axis locating at Z = 0. Noting that, since X(ψ, θ) and Z(ψ, θ) are given when introducing the (ψ, θ) grids, one can also determine S r (ψ, θ) and Θ(ψ, θ). Consequently, one can derive both Using these results one can compute the metric parameters in the (S r , Θ, φ) coordinate system.
We first work on the Jacobian J . Note that Note further that One obtains the Jacobian expression in the polar coodinates Next, we work on other metric parameters. By straightforward reduction one can obtain where it has been noted that Noting further that one obtains Using the toroidal symmetry property, we can also find that |∇φ| 2 = 1 X 2 , ∇ψ · ∇φ = 0, and ∇θ · ∇φ = 0.
The expressions of metric parameters given above can be used to compute the matrix F and the vector s in Eq. (11).

C. Iteration scheme and boundary conditions
With the computation of metric parameters given in the last subsection, we describe the iteration scheme to solve the Grad-Shafranov equation with proper boundary conditions.
Since the equation, Eq. (11), are nonlinear, an iteration process is necessary. One can follow the usual iteration scheme to get the converged solution: Here, n denotes the iteration step.
Equation (12) is a set of inhomogeneous differential equations of first order. Its general solutions at step n + 1 can be expressed as where c k are the complex constants to be determined by the boundary conditions, u k are the independent solutions to the homogeneous equations and u s is the specific solution to take into account the source term s on the right hand side of Eq. (12). For brevity, the step index n has been dropped.
Since the number of equations is 2M , the solutions are completely determined by the M boundary conditions in complex at the magnetic axis and M boundary conditions in complex at plasma edge χ a . The boundary conditions at plasma edge χ a are specified by the given shape of the last closed flux surface in the fixed boundary value problem. The boundary conditions at the magnetic axis are just the requirement that the independent solutions are "small" in terms of the terminology of differential equation theory. The "large" solution causes the system energy to diverge, while the "small" solution is square-integrable with respect to the energy integral. Near the magnetic axis, the homogeneous part of the Grad-Shafranov equation can be approximated by the cylinder model. In this limit the solutions are as follows [39]: where r is the minor radius and a m and b m are constants. Therefore, the boundary conditions for small solutions are simply b m = 0. This yields The boundary conditions for A 1m can be obtained using the definition of A 1 in Eq. (6).
Note that the general solution to the set of differential equations is the summation of homogeneous solutions and specific solution and the boundary conditions are satisfied by the constants c k tied to the homogeneous solutions. Therefore, the boundary conditions for specific solution are arbitrary.

D. Solution of equilibrium equations
The principle to solve Eq. (12) is laid out in subsection III C. The actual implementation is more complicated. One needs to divide multiple regions in the radial direction and then match the solutions in the individual regions to get the global solution. In this subsection, we will outline the actual numerical process in the ATEQ code to solve the Grad-Shafranov equation.
The M boundary conditions at the magnetic axis in Eqs. (14) and (15) where Matrix Y is a band matrix. By inverting it one can obtain the solution of Eq. (16) With the constants obtained from Eq. (17), the solutions in each region are then simply These give the numerical scheme being implemented in the ATEQ code to solve the Grad-Shafranov equation.

IV. NUMERICAL PROCEDURE AND RESULTS
In this section, we describe how to implement the numerical scheme in section III. This leads to the development of ATEQ code. The computational flow chart of ATEQ is given in Fig. 1. To be more specific to describe the computational flow, we use an ITER-like equilibrium as an example. The major radius 6.2 m, minor radius 2 m, elongation 1.78, triangularity 0.4, the vacuum magnetic field at the geometric center of plasma column is 6 T, the total current 15.9 MA, and the volume average beta value is 3.371%. Figure 2 shows the cross section with the "a" part showing the initial grid setup and the "b" part showing the magnetic surfaces computed by the ATEQ code. The case will also be used for the benchmark studies with the TOQ code. Further details about the equilibrium will be described there.
First, one needs to set up radial and poloidal grids (ψ, θ) as shown in Fig. 2a. The grids are constructed to surround the magnetic axis (x axis , z axis ). Because the magnetic axis is unknown beforehand, iteration is needed. The value of the previous step (n) is used to construct the grids to advance to the next step n + 1. Following the iteration scheme in Eq. (12), the source term on the right hand side of Eq. (12) is evaluated by using the solution for poloidal flux u (n) (ψ (n) , θ (n) ) in the previous step. Note that the pressure and current profiles are prescribed by the normalized poloidal flux. The total poloidal flux χ a needs also to be determined iteratively. At the first step, the quantities at the previous step are prescribed by initial guessing. The matrices F and s can then be computed with previous step grids as described in subsection III B. Using the splines the matrices F and s are made to be radially continuous functions.
Here, it is noted that the proper choice of initial (ψ, θ) grids can affect how many Fourier components are required. For the usual equilibria without X points included the choice is rather arbitrary, i.e., a wide range of grid choices can work well. For the equilibria with X points included proper choice of initial grids is important. In the ATEQ code, the initial grids are specified as follows. First, the grids with ellipticity k and triangularity δ are set up inside the specified plasma-vacuum boundary according to the formula with r = [(X − x axis ) 2 + (Z − z axis ) 2 ] 1/2 . Here, k and δ can be polynomial functions of ψ.
This means that one can adjust the ellipticity and triangularity from the axis to the outmost surface. In most cases, the linear dependence is sufficient. Next, the difference between the specified plasma boundary and the outmost surface given by Eqs. (19) and (20) are distributed radially. The distribution can be adjusted through an exponential multiplier of ψ. Also, the ψ grids can be packed near the axis and boundary. In our experience, with these flexibilities, roughly 100 Fourier sidebands are sufficient to get a good equilibrium solution with X points included. It is using this type of initial grid setting that the Solovév solution with X points to be described later is reproduced numerically. There is always a possibility to use the (ψ, θ) solution at step n for the grids at step n + 1. Nevertheless, it can only be used if the solution at step n is sufficiently smooth and well-behaved.
Next, the whole radial domain is split into L regions. As described in subsection III D, adaptive shooting is implemented to get the independent solution matrices in each region.
By solving for l c M l using Eq. (17), one can construct the global solution through Eq. (18). At this step, we first check if the magnetic axis (x axis , z axis ) and total poloidal flux χ a converge.
Usually, total poloidal flux converges in one or two steps, using the following formula for prediction Instead, to find the magnetic axis (x axis , z axis ) one needs a few iterations. The code shoots outwardly from the assumed magnetic axis (x axis , z axis ). After achieving the solution, the minimum of poloidal flux χ is determined. The location of this minimum is used as (x axis , z axis ) for the shooting in the next step. This process is repeated until the starting (x axis , z axis ) matches the location of the χ minimum computed to a required accuracy.   Fig. 2 uses 50 sidebands.
With the magnetic axis (x axis , z axis ) and total poloidal flux χ a being converged, one can further iterate to get the converged solution χ(ψ, θ). With this solution, one can obtain the numerical solution for the poloidal magnetic flux χ(X, Z). The magnetic surfaces with χ(X, Z) = const are plotted in Fig. 2b.

V. BENCHMARK STUDIES AND COMPARISON WITH EXISTING CODES
This section describes the benchmark studies. We begin with the analytical Solovév equilibrium with the X points included [41]. The benchmark with the Solovév solution is not a trivial task. This is because the X-points are present in the Solovév equilibrium. The equilibrium computation with the X-points included is challenging because much more Fourier components are needed. The analytical Solovév solution is given as follows where the parameters are given as follows in the benchmark studies: X 0 = 10, a = 1, As pointed out in Ref. [41], the second-order solution in Eq.   Fig. 5. Two curves completely overlap, although the initial guessing as shown in Fig. 4a deviates dramatically from the actual solution in Fig. 4b in the ATEQ computation.
Comparisons with the existing equilibrium codes are also performed. Here, we describe a benchmark example with the TOQ equilibrium code [16]. A typical case is described as follows. A TOQ sample initiation file with equiltype = f f prime is taken. To be more specific to compare with TOQ, here the same numerical parameter notations as in the TOQ manual are used to describe the equilibrium parameters. The shape of boundary type is specified by ishape = 2, which is described as follows: where the basic parameters are specified as follows: The major radius rzero = 6.2 m, the minor radius rmax = 2 m, the elongagtion eshape = 1.78, and the triangularity xshape = 0.4. This leads the equilibrium cross section to be given in Fig. 2.
The pressure gradient (p ) and poloidal current flux parameter (gg ) profiles are specified, respectively, by setting modelp = 3 and modelf = 1, which are described as follows: Note here that the profiles are specified with the normalized poloidal magnetic fluxχ, varying from 0 to 1 from the magnetic axis to plasma boundary. The tolerance is set to be toleq = 10 −5 in the TOQ iteration with successively increasing grid densities. Here, we have used the nonlinear pressure profile in Eq. (24), which is different from the TOQ sample initiation file, in order to avoid the linear profile case considered in the Solovév case. The pressure and poloidal current flux profiles are given in Fig. 6 with respect to the minor radius on the outer vertical mid-plane. Although the p and gg are the same as specified in Eqs. (24) and (25) for TOQ and ATEQ codes, the pressure (p) and poloidal current flux (g) profiles can be slightly different since they are given in the minor radius, instead of the normalized poloidal flux. The slight difference of poloidal magnetic flux solution as discussed later can cause the difference. The volume average beta is 3.371%, the normalized beta is 2.54, and l i = 0.730 in this equilibrium. The decomposition with independent solutions effectively reduces the matrix size for solving the magnetohydrodynamic equilibrium problem, as compared with numerical schemes based on a fixed radial grid. The adaptive numerical scheme is expected to be especially helpful to deal with stiff equilibrium problems. Our results also indicate that the backward substitution method can be necessary to obtain a reliable equilibrium solution.
Let us here further discuss the unique features of the ATEQ numerical scheme. The numerical methods for solving the Grad-Shafranov equation ultimately reduce the problem to solve the matrix equations. The matrix size then matters and reducing the matrix size in discretizing the Grad-Shafranov equation is important. In the grid-based numerical schemes both in the radial and poloidal directions (finite difference or finite element), the size of the matrix is N r × N t . Here, N r and N t are respectively the numbers of grids in the radial and poloidal directions. In the numerical scheme based on the poloidal Fourier decomposition, the matrix size is N r × N f . Here, N f is the number of the poloidal Fourier components. To achieve high accuracy, especially for tough problems related to the axis, X-point, or pedestal, etc. one has to increase the N r and N t (or N f ). Consequently, the size of the matrices becomes large and the matriices become hard to deal with numerically. In the ATEQ numerical scheme, the radial direction is split into L regions with each region addressed by the adaptive shooting of independent solutions. It reduces the radial N r grid problem into a L region matching problem. This cuts down the N r × N t (or N r × N f ) matrix problem in the conventional numerical schemes into a L × N indep problem in the ATEQ numerical scheme.
Here, the number of regions L is about a few 10s and N indep is the number of independent solutions, which is of the same order as N f . The reduction of the matrix size is by the factor L/N r , which is about an order of magnitude, as compared with the conventional radially grid-based numerical schemes. Also, in this ATEQ numerical scheme, no matter how accuracy in the radial direction is imposed, the size of the matrices basically does not change. Such an improvement in the order of magnitude rarely happens. It therefore represents a significant development in this research.
The equilibrium problem is a little bit different from the stability one. If one uses the exact flux solution as the radial grids, only a single Fourier component for the magnetic flux χ is required because it is constant on the surfaces. Therefore, the required number of the Fourier components, N f , in principle can be somewhat optimized by setting proper radial grids. Since the matrix size is reduced in the radial direction in ATEQ, one has more flexibility to increase the umber of poloidal Fourier components if it is required. This is a distinct feature of ATEQ as compared to the conventional Fourier decomposition based codes. This improvement is useful.
It is realized in this field that a good numerical equilibrium solution near the axis and plasma boundary in the presence of the X points is critical. It is a challenging issue for decades. As cited in the introduction, several efforts have been made. Our work provides another possible solution. To directly compare with other codes to treat the X point equilib-rium problem will be our next task. This may require close collaboration with other teams.
Equilibrium codes often need certain specific procedure to execute them. Using the backward substitution method we found that the equilibrium accuracy varies a lot even with the same code. That's why we're wary of doing code-to-code comparisons directly without the other party involved. Each code may have their own particular features. We have compared with TOQ since the example file is in the public domain. Even in this case, we have provided  b) The converged magnetic flux surfaces.