{'fraction_non_alphanumeric': 0.04098929890463714, 'fraction_numerical': 0.013400044193314182, 'mean_word_length': 4.187407892582282, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
ECOM: a fast and accurate solver for toroidal axisymmetric MHD equilibria

We present ECOM (Equilibrium solver via COnformal Mapping), a fast and accurate fixed boundary solver for toroidally axisymmetric magnetohydrodynamic equilibria with or without a toroidal flow. ECOM combines conformal mapping and Fourier and integral equation methods on the unit disk to achieve exponential convergence for the poloidal flux function as well as its first and second partial derivatives. As a consequence of its high order accuracy, for dense grids and tokamak-like elongations ECOM computes key quantities such as the safety factor and the magnetic shear with higher accuracy than the finite element based code CHEASE [H. L\"utjens \textit{et al.}, Computer physics communications 97, 219 (1996)] at equal run time. ECOM has been developed to provide equilibrium quantities and details of the flux contour geometry as inputs to stability, wave propagation and transport codes.


I. INTRODUCTION
Numerically computed magnetohydrodynamic (MHD) equilibria are the starting point of a wide class of numerical solvers that are used to study MHD stability, transport, and heating and current drive in magnetic fusion devices [1][2][3][4]. Static MHD equilibria of toroidally axisymmetric configurations are described by the Grad-Shafranov (G-S) equation [5,6], a nonlinear, second-order elliptic partial differential equation. Stationary equilibria with purely toroidal flows are determined by solving a close variant of the G-S equation [7], the only difference being that for the latter the pressure term does not only depend on the poloidal flux function Ψ, but also on the radial variable R. Numerical codes to solve the G-S equation have been developed since the early days of the magnetic fusion program [7][8][9]. Nevertheless, the development of optimized G-S codes remains a topic of active research, for three main reasons. First, G-S solvers must be able to properly resolve complex two-dimensional geometries [10,11], with boundaries that may have a corner, corresponding to a magnetic field X-point [11]. Second, G-S solvers must be fast. This criterion is particularly relevant in the context of multiphysics integrated simulations [12][13][14][15][16]. Several of these multiscale, multiphysics solvers already include, or will eventually include, in their iterative procedure a step in which the equilibrium configuration is self-consistently updated. A reasonable requirement is that the calculation of the updated equilibrium takes a negligible amount of time and computing resources as compared to the computationally intensive transport, MHD stability and plasma heating solvers. Third, G-S solvers must be accurate. The solution of the G-S equation is the poloidal flux Ψ, but the physical quantities of interest, such as the magnetic field, the safety factor, the magnetic shear, the magnetic curvature, and the current density are all functions of partial derivatives of Ψ. Since there always is some loss of accuracy when computing derivatives, a high level of accuracy for Ψ is desired.
In this article, we present the new Grad-Shafranov code ECOM (Equilibrium solver via COnformal Mapping). ECOM is a fixed boundary, direct solver written in Fortran 77/90 that is based on three key elements: 1) the formulation of the G-S equation as a nonlinear Poisson problem; 2) a spectrally accurate numerical method to compute the conformal map from the smooth plasma cross section of interest to the unit disk; 3) a fast, high order Poisson solver on the unit disk [17]. Its main novelty lies in the last two aspects discussed in the paragraph above, namely accuracy and speed. Regarding the first point, we demonstrate in this article that ECOM has better convergence properties than popular G-S solvers based on finite elements [10,18,19]. In the finite element approach the numerical error of the solution decays as a power of grid size, i.e. N −α where α is an integer and N is the number of grid points in one direction. Often, for magnetic fusion applications α ≤ 4 [9,18], although α ≤ 7 was recently demonstrated [10]. In contrast, convergence in ECOM is exponential: the error decays as β −N for some real number β > 1. Just as importantly, in ECOM the rate of convergence for the derivatives of Ψ is the same as that of Ψ, whereas in the finite element approach the derivatives of Ψ converge more slowly than Ψ [10]. Remarkably, in ECOM numerical accuracy is not obtained at the expense of computational complexity and speed. For a given grid size, our solver is faster than finite element solvers and less demanding in terms of memory. We will show that a drawback of relying on conformal mapping is that ECOM often requires a denser grid than FEM based solvers to achieve a desired accuracy. Yet despite this, we find that for tokamak geometries and medium to high number of grid points, ECOM is more accurate than FEM based equilibrium codes at equal run time.
This article follows an earlier article [17], in which we gave a detailed description of our new numerical algorithm for solving the G-S equation. The focus here is different. One of the main motivations is to present extensions recently added to our G-S solver that make it a practical tool readily usable in fusion applications. The new capabilities of our equilibrium solver include the possibility of computing equilibria with arbitrary toroidal flow profiles, the possibility of specifying current and pressure profiles in various ways, as well as the evaluation of the key physical quantities that are required as inputs in stability, transport and heating codes. Since ECOM is a direct solver that calculates Ψ on a prescribed grid for the poloidal cross section [20], we put a particular emphasis on the accurate computation of the contours of constant flux, and of flux surface quantities such as the safety factor and the magnetic shear. A second motivation for this article is to perform detailed comparisons between the popular G-S code CHEASE [18] and ECOM, and to assess the merits of each solver.
The structure of the article is as follows. In Section II we briefly review the numerical algorithm we use to solve the G-S equation [17]. In Section III, we give a detailed presentation of the equilibrium quantities ECOM computes during the postprocessing phase, and of the numerical methods we implemented to calculate these quantitites with high accuracy. In Section IV we evaluate the speed, accuracy, and convergence properties of our solver, and compare them to those of CHEASE [18]. In Section V we explain how ECOM computes stationary equilibria with toroidal flows, and in Section VI we summarize our main findings, discuss the current limitations of ECOM and future plans. Appendix A presents our method to calculate the Miller parametrization [21] of a numerically computed flux contour, and Appendix B contains a table with all the important variables in ECOM, along with a short description for each of them.

II. NUMERICAL ALGORITHM
In this section, we briefly review the numerical algorithm used in ECOM to solve the G-S equation. A more detailed presentation of each of the steps described below can be found in [17].

A. The Grad-Shafranov equation as a nonlinear Poisson problem
The Grad-Shafranov equation is given by where (R, φ, Z) is the usual cylindrical coordinate system associated with the toroidal geometry, 2πΨ is the poloidal magnetic flux, µ 0 is the permeability of free space, p(Ψ) is the plasma pressure, and F (Ψ) = RB φ , with B φ the toroidal component of the magnetic field. Once the free functions p(Ψ) and F (ψ) are given and Eq. (1) is solved with appropriate boundary conditions, the magnetic field B and the current density J can be computed according to the following formulae (1) is a second-order elliptic nonlinear partial differential equation for Ψ. ECOM solves the fixed boundary problem associated with this equation. Specifically, the boundary curve ∂Ω enclosing the plasma domain Ω of interest is an input to the solver, and ECOM solves Eq. (1) with the Dirichlet data Ψ = Ψ b on ∂Ω, where Ψ b is a constant. This formulation is particularly convenient for multiphysics theoretical studies of the influence of shaping on plasma performance [22][23][24]. Two types of inputs can be used in ECOM to determine the geometry of ∂Ω. One option is to give an exact representation of the plasma boundary, for example in the form of parametric equations [21]. When such a representation is not available, one can also give the coordinates (R n , Z n ) of discrete points on the boundary. At the moment, ECOM can only treat smooth plasma boundaries, and can therefore not compute equilibria whose plasma boundary has a separatrix.
The functional dependence on Ψ of the pressure and toroidal magnetic field profiles is either prescribed or determined from transport equations. In both cases, it is an input to ECOM. In general, these profiles are such that Eq. (1) is nonlinear, and for a wide class of profiles Eq. (1) has to be solved as an eigenvalue problem [8,17,25,26]. This means that Eq. (1) has to be solved by iterating on Ψ [27]. In ECOM, this is done as follows. A normalized flux ψ is defined by ψ = (Ψ − Ψ b )/(Ψ 0 − Ψ b ), where Ψ 0 and Ψ B are the poloidal flux at the magnetic axis and the last closed flux surface respectively, so that ψ = 1 at the magnetic axis and ψ = 0 at the last closed flux surface. The pressure and toroidal magnetic field profiles are also normalized and expressed in terms of ψ according to: dp(Ψ) dΨ = dp(ψ) dψ and Defining where λ plays the role of an eigenvalue. In ECOM, there are several options to specify the profiles dp/dψ and dF 2 /dψ, with corresponding namelist parameter IPTYPE for the pressure and IFTYPE for the poloidal current. If IPTYPE=1 or IFTYPE=1, the profiles are given by an explicit formula in terms of ψ. In ECOM, we often use dp/dψ = p 0ψ (1 − (1 − ψ) pin ) pout , as is also done in CHEASE [18], where the constants p 0ψ , p in , and p out are specified in the namelist. Likewise, we often use dF 2 /dψ = 2F 0ψ (1 − (1 − ψ) Fin ) Fout . Different expressions can be easily implemented, such as formulae describing a steep pressure pedestal [17]. If IPTYPE=2 or IFTYPE=2, the profiles are given by a set of data points and the corresponding values of ψ or of the minor radius. The value of the namelist variable IPTABLE determines whether the tabulated values of the profiles are in terms of ψ or of the minor radius. If IPTABLE=0, the numerical tables of dp/dψ and dF 2 /dψ are specified in terms of discrete values of ψ. If IPTABLE=1,p(ψ) andF (ψ) instead of their derivatives are specified by tables in terms of ψ. If IPTABLE=2, the numerical tables ofp andF are specified in terms of the normalized minor radius ρ. In Section III, we explain in detail how ECOM accurately computes dp/dψ and dF 2 /dψ starting from such tables. Furthermore, ECOM offers the possibility to choose between three different definitions for the minor radius ρ, corresponding to three different values of the namelist variable IRHO. If IRHO=0, the minor radius is defined by is the outer location of the flux surface ψ at Z = 0 and R 0 is the location of the magnetic axis. If IRHO=1, ρ is defined by As already mentioned, Eq. (3) must be solved iteratively. G-S solvers usually iterate on the operator ∆ * [10,18]. One of the key ideas in ECOM is to iterate on the Laplacian operator ∆ instead. This can be done without any loss of generality, and does not incur any additional computational cost. The advantage is that at a given iteration, one now has to solve Poisson's equation, and one can therefore rely on the larger body of numerical methods developed for fast high order Poisson solvers in two dimensions. To solve for ψ and the smallest eigenvalue λ in Eq.(3), ECOM uses a modified version of the inverse iteration method [27]. Specifically, if ψ (i) and λ (i) are known at the iteration step i, then ψ (i+1) and λ (i+1) are computed according to where ||ψ i+1 || ∞ is the extremum ofψ i+1 in the domain Ω. In ECOM, the iterative process terminates when ||ψ (i+1) − ψ (i) || ∞ < δ for some predetermined small δ. In the numerical results presented in this article, we typically had δ = 10 −14 . In order to keep consistency across different grid sizes, ||ψ i+1 || ∞ is not only calculated for all values of ψ i+1 on the mesh, but instead over the entire domain Ω. In ECOM, this is done by finding the location where ∇ψ = 0 with the Newton-Raphson method. It requires knowledge of the Hessian matrix, whose values away from grid points are evaluated by interpolation, based on the high order Fourier and Chebyshev representations ECOM uses for the Poisson solver on the unit disk. The extremum of ψ on the grid is used to provide a very good initial guess, so that in practice very few Newton steps are subsequently required to find ||ψ (i+1) || ∞ in Ω.
It is known empirically that iterative schemes such as the one above converge faster when the right-hand side of the partial differential equation is slowly varying. For better convergence we thus scale the unknown function ψ as ψ = u √ R. Replacing ψ with u, Eq.(4) becomes Unlike Eq. (4), the right-hand side of Eq. (7) does not have derivatives of ψ, and it is therefore smoother. Note that while one solves for u in the Poisson step, the normalization steps (5) and (6) are still computed in terms of ψ.
Thus far, we have explained how ECOM treats the Grad-Shafranov equation as a nonlinear Poisson problem. We now briefly describe how ECOM computes the solution of the Poisson equation (7) at each iteration, i.e. with fixed right-hand side, on a domain Ω of fusion interest and with the Dirichlet boundary conditionũ = 0 on ∂Ω. The Poisson solver in ECOM is based on two elements: 1) a spectrally accurate numerical method to compute the conformal map from the plasma domain Ω to the unit disk; 2) a fast, high order Poisson solver on the unit disk.

B. Conformal maping from the plasma domain to the unit disk
Conformal mapping is an effective method for solving Poisson's equation because a conformal map transforms a Laplacian operator into another Laplacian operator, with a scale factor [9]. Consider the generic Poisson equation in Ω u(R, Z) = 0 on ∂Ω the conformal map W : z = R + iZ → w = α + iβ from Ω to the unit disc D 1 and its inverse map B : w = α + iβ ∈ D 1 → z = R + iZ ∈ Ω. Solving Eq. (8) is equivalent to solving the following Poisson problem in D 1 : where u(R, Z) = v(α(R, Z), β(R, Z)), and the functions R(α, β) and Z(α, β) should be seen as the real and imaginary parts of the inverse map B. Clearly, solving Eq. (9) is easier than solving Eq. (8), provided one has a way to calculate the inverse map B at both boundary and interior points of D 1 . In ECOM, this is done as follows. ECOM first computes the forward map W for points on ∂Ω that are equispaced in arc length through the Kerzman-Stein integral equation based on the Szegö kernel [17,28,29]. Using oversampling and interpolation, ECOM then uses the boundary values of the forward map to calculate the inverse map R(α, β) and Z(α, β) for points on the boundary of D 1 that are equispaced in the polar angle ϑ. Finally, ECOM computes B for points in the interior of D 1 using the Cauchy integral formula and the Fast Fourier Transform [17]. ECOM relies on a somewhat naive implementation of the Kerzmann-Stein integral equation for the computation of W on ∂Ω, that requires O(n 3 1 ) work, where n 1 is the number of discretization points on ∂Ω. There exist methods resulting in an asymptotic O(n 1 ) run time [30], but they are not currently implemented in ECOM because the computation only needs to be done once, and because other steps in ECOM are more expensive. The computation of the inverse map at interior points is based on the Cauchy integral formula and the Fast Fourier Transform, and results in a run time complexity of O(n r n ϑ logn ϑ ), where n r is the number of radial grid points and n ϑ the number of angular grid points for the mesh that ECOM uses to solve Poisson's equation on D 1 The boundary ∂Ω can be defined in several ways, that are specified by the namelist variable IBTYPE in ECOM. If IBTYPE=0, the plasma boundary corresponds to the contour Ψ = 0 of a Solov'ev equilibrium we discuss in Section IV A, parametrized by Eqs. (44)- (45). If IBTYPE=1, the boundary is specified by the Miller parametrization [21] given by Eqs. (52)-(53) of Section IV B, and the elongation κ and triangularity δ m must then be specified. If IBTYPE=2, the boundary is specified by a set of discrete points (R, Z). This allows the computation of equilibria specified by experimental data and is also the method of choice to compute up-down asymmetric equilibria. The conformal mapping routine requires that the points on ∂Ω be equispaced in arc length. When IBTYPE=0 or IBTYPE=1, these points are easily calculated from the parametric equations for the boundary. When IBTYPE=2, ECOM uses Lagrange interpolation to compute these points.
In principle, the conformal map has to be computed only once, at the beginning of the iterative procedure corresponding to Eqs. (5)- (7). However, as one might physically expect, computing flux surface quantities is much more convenient if the point that is mapped to the center of D 1 coincides with the magnetic axis. To facilitate the calculation of these quantities, which takes place after the G-S equation is solved, ECOM recomputes the conformal map several times within the Poisson iterations to adjust the center of D 1 to the magnetic axis. We have empirically observed that the conformal map only needs to be recomputed a few times.
The Riemann mapping theorem guarantees the existence and uniqueness of an analytic map between any simply connected plasma cross section Ω and the unit disk. However, this does not mean that conformal mapping is a practical numerical method for any arbitrary plasma shape. The issue is that points on the boundary of D 1 that are equispaced in the angle ϑ are not necessarily mapped, under the inverse map, to points that resolve the boundary of Ω in the desired fashion. It is well known, for example, that if Ω is an elongated ellipse, uniformly spaced points in ϑ on the boundary of D 1 correspond to a distribution of points on the boundary of Ω which is sparse on the curved parts and crowded on the flat parts of the ellipse [9,17]. In the remainder of this article, we will call this phenomenon the "crowding effect". It has two direct implications for ECOM. First, ECOM can only treat in a robust manner domains that have a smooth boundary, and can therefore not be used to compute equilibria with a magnetic X-point on the plasma boundary. Second, ECOM is particularly efficient for plasma shapes that are not too elongated, as is the case for conventional tokamaks. As elongation is increased, the high order convergence properties are maintained for Ψ and its derivatives, but a higher number of grid points is required to reach a certain level of accuracy [17]. We will go back to this point in Section VI.

C. Fast, high order Poisson solver on the unit disk
We finish this section by describing how ECOM solves Poisson's equation on the unit disk, as required in Eq. (9). The solver uses separation of variables in the usual polar coordinates (r, ϑ) and expands v and g as Fourier series Substituting these expressions into Poisson's equation, we get the following ordinary differential equation for each n: where the boundary condition at r = 0 is obtained by requiring the regularity of the solution at this point. For each n, a solution of Eq. (11) that does not satisfy the boundary condition at r = 1 can be written in terms of convolutions with the Green's function associated with Eq. (11) that has the proper behavior at r = 0 and r → ∞: For each n, the general solution to the homogeneous equation satisfying the regularity condition at r = 0 can also be written explicitly: v H n (r) = c n r |n| (14) where c n is a constant to be determined from the boundary condition at r = 1. Setting c n = −v P n (1), we can then write the solution of Eq. (11) satisfying the proper boundary conditions as follows: One of the major advantages of the Green's function formulation and the formula (15) used in ECOM is that partial derivatives of v can be calculated explicitly from the formulae in Eqs. (12) and (14) [17]. Numerical differentiation is never required, which is one of the main reasons why our numerical method leads to partial derivatives of ψ that have the same order of convergence as ψ.
In ECOM, the computation and sum of the Fourier series are done with the Fast Fourier Transform. The angular grid on the unit disk is uniformly spaced in the polar angle ϑ, to guarantee the spectral accuracy of the representation for smooth data. The number of grid points in the ϑ direction is n ϑ . In the radial direction, ECOM uses a piecewise Chebyshev grid. Specifically, the interval [0, 1] is divided into n L subintervals, and on each of these subintervals a Chebyshev grid of order n ch is constructed. The number of grid points in the radial direction is then n r = n ch n L . n ch = 16 is the default setting in ECOM. The convolutions with the Green's function in Eq. (12) are computed with a 16th order Gaussian quadrature rule. Reference [17] describes ways to avoid the computational issues associated with the rapid growth and decay of the monomials s |n| and s −|n| for large n, as well as the recursive algorithm used to compute these integrals in O(n r ) work. The run time complexity of the Poisson solver on the disk is O(n r n ϑ logn ϑ ): ECOM computes O(n r ) FFTs of size n ϑ at a cost of O(n r n ϑ logn ϑ ) and solves n ϑ radial ordinary differential equations at a cost of O(n r n ϑ ).
The flow diagram in Figure 1 presents a condensed view of the iterative scheme used in ECOM, as a summary of Section II. The initialization step corresponds to the specification of the grid resolution and of the values of the namelist variables presented in this article. The parameters n f , n θD and n θE in Figure 1 refer to discretizations used during postprocessing, after the G-S equation is solved, and are defined in Section III.

III. POSTPROCESSING
The main purpose of postprocessing is to compute equilibrium quantities that play a key role in heating and current drive, stability and transport calculations, and to scale the normalized solution ψ to the physical flux Ψ. In addition, ECOM can compute the parallel current density using reduced models for the Ohmic current and the bootstrap  [31] and the global Troyon limit [32], and calculate the Miller parametrization that best fits a flux surface chosen by the user.

A. Flux functions
The evaluation of flux quantities requires integration along contours of constant poloidal magnetic flux ψ. A flux coordinate system (ψ, θ, φ) is more convenient for such integrals than the (R, φ, Z) coordinate system used by ECOM to solve the G-S equation. For postprocessing, ECOM thus constructs the following flux coordinate system. φ is chosen to be the usual toroidal angle φ, with ∇φ = e φ /R, while the poloidal angle θ is defined by where atan2 is the four-quadrant inverse tangent and (R 0 , Z 0 ) is the position of the magnetic axis. In terms of the (ψ, θ, φ) coordinates, the flux surface average of a function X [7], written X , is defined by where J = (∇ψ × ∇θ · ∇φ) −1 is the Jacobian of the transformation between Cartesian coordinates and the (ψ, θ, φ) coordinate system: As will soon be apparent, several physical quantities are best expressed in terms of the three functions I A (ψ), I B (ψ), and I C (ψ) defined by Rewriting Eq. (3) using flux coordinates, a simple relationship can be obtained that links I A , I B , and I C . Indeed, Eq. (3) takes the following form in flux coordinates [7]: Taking the flux surface average of 1/R 2 × Eq. (19), we obtain the desired relationship between I A , I B , and I C : The plasma volume inside the flux surface ψ is and the total plasma volume is V 0 = V (ψ = 0). The toroidal current within the flux surface ψ is where we have used Eq. (20) and I A (ψ = 1) = 0. The total toroidal current in Ω is I p = I φ (ψ = 0). The volume averaged pressure p V is given by where we have used integration by part, combined with dV /dψ = 2πI B , V (ψ = 1) = 0, and p(ψ = 0) = 0. The total beta is where Our definitions of the poloidal beta and of the internal inductance are the same as Jardin's [7]: The total poloidal magnetic field energy is Finally, for each flux surface, the safety factor is defined as The normalized radius can be used to find the differential volume dV (ψ)/dρ and the differential flux dψ/dρ, which are often used in transport or MHD analysis. Also, the magnetic shear can be defined in terms of ρ bŷ B. Numerical method for contour integrals To evaluate the integrals in Eq. (18) numerically, one needs to find the location of the desired flux contours, which in general do not coincide with the (R, Z) grid of the Poisson solver, and then integrate the integrands along these contours. In ECOM, this is done with the following three steps: (i) ECOM first determines the radial location of n f flux surfaces for each angle ϑ in the unit disk D 1 ; (ii) ECOM then evaluates the integrands at the corresponding (R, Z) points in Ω, and interpolates the integrands to a grid that is equispaced in the angle θ; (iii) ECOM finally computes the integrals along contours of constant ψ in the domain Ω using the trapezoidal rule. Since the numerical methods that are used for interpolation and integration in steps (i)-(iii) are spectrally accurate, ECOM computes the location of the flux contours and the integrals in Eq. (18) without significant loss of accuracy, as we will demonstrate in Section IV. We now describe steps (i), (ii), and (iii) in more detail.
For a description of step (i), we put ourselves in the situation in which the G-S solver has computed the values ψ(r i , ϑ j ) of the flux ψ on the grid of the unit disk D 1 , and we imagine that we want to determine the location of the contour ψ = ψ s on D 1 . ECOM does this as follows. For each angle ϑ j , ECOM first finds the Chebyshev subinterval of the radial grid for which ψ(r t , ϕ j ) > ψ s > ψ(r t+ch−1 , ϕ j ), where t is the index of the first Chebyshev point in that radial subinterval. Once the subinterval is found, ECOM constructs a local continuous approximation ψ c (r, ϑ j ) of ψ in the radial direction from the known values ψ(r i , ϑ j ) on the subinterval of interest and the Chebyshev grid for that subinterval. Specifically, ψ c is written as the following sum: where the functions T k are the Chebyshev polymials associated with the Chebyshev grid of the subinterval, and where the coefficients a k are given by the expression where δ ij is the kronecker delta. ECOM then uses the expansion in Eq. (30) to find the radial position satisfying |ψ c (r) − ψ s | < δ with a Newton-Raphson iterative method: This root finding process usually converges in a few iterations, typically m ≤ 5, and the total cost to find the location of n f contours at n ϑ angles is O(n f n ϑ n ch ). Finally, after the radial position of a given contour is found, ECOM uses (∂ψ/∂r)(r i , ϑ j ) and (∂ψ/∂ϑ)(r i , ϑ j ) for i = t, t + 1, ..., t + n ch − 1 and Chebyshev representations analogous to Eq. (30) to accurately evaluate (∂ψ/∂r)(r, ϑ j ) and (∂ψ/∂ϑ)(r, ϑ j ) at the location of the flux contours. At the end of step (i) the radial position of the specified flux contours are known for each angle ϑ j in D 1 , and so are the values of the integrands in Eq. (18) at these points. Through the backward map, all these quantities are also known in Ω. In Eq. (18), one integrates quantities that are 2π-periodic in θ over the period [0, 2π]. Numerically, this can be done very accurately with a trapezoidal-rule quadrature, provided that the integrands are known on an equispaced θ grid. Since the equispaced ϑ grid in D 1 is not mapped to an equispaced θ grid, the goal of step (ii) is to interpolate the quantities computed in step (i) on a grid in Ω that is equispaced in θ. Because of the crowding effect that is inherent to the conformal mapping technique, the angular grid resulting from the inverse map underresolves certain regions of Ω. For accurate interpolation, it is thus desirable to first oversample the integrands to be interpolated. The oversampling

Namelist
Constraint to rescale λ Constraint on q profile ISCALE=0 None is done by refining the equispaced ϑ grid of D 1 using the FFT, assuming band-limited integrands. As a result of this, the integrands are known at n ϑ k samp angular grid points, where k samp is the oversampling factor. In principle, these values could then be mapped back to Ω, and interpolated. In practice, however, the backward mapping of the n ϑ k samp grid points for n f contours from D 1 to Ω requires O(n f (n ϑ k samp ) 2 ) operations, and results in significant computational time when k samp is as large as desired for accurate interpolation. Note that the run-time complexity is not of the form n f n ϑ k samp log(n ϑ k samp ) as it was for the backward mapping of the (r, ϑ) grid points in D 1 because the calculation of the backward map for the flux contours cannot be trivially accelerated by the FFT. The reason for this is that the radial location of the flux contours varies as a function of ϑ. In order to reduce the computational cost, we only compute the backward map for n θD points among the n ϑ k samp oversampled points, with n θD n ϑ k samp , and chosen such that their mapped positions are the closest to the target grid of n θE points equispaced in θ. These n θD points are found by computing the backward map of a small number of points lying on the contour of interest in D 1 with increasing angle ϑ across the interval [0, 2π], and using k Lag -th order Lagrange interpolation to construct an approximation of the function ϑ(θ) on that contour. We usually take n θD = n θE and k Lag = 8 k samp . Once the n θD non-equispaced points are found, we use trigonometric interpolation for a periodic function [33] to interpolate the integrands at the n θD points to the n θE equispaced points. This requires O(n f n 2 θD ) work to find the barycentric factors, and O(n f n θD n θE ) work to interpolate at n θE points.
Once step (ii) is completed, step (iii) is straightforward. The contour integrals in Eq. (18) are computed from the integrands on the equispaced θ grid using the trapezoidal rule. Since the integrands are smooth and periodic in θ, and since the θ grid is uniform, the trapezoidal rule is spectrally accurate [34]. The required work for the trapezoidal-rule quadrature is very small, O(n f n θE ).

C. Scaling the equilibrium
A single solution of the normalized form of the G-S equation as given in Eq. (3) can describe an infinite sequence of axisymmetric equilibria that have a different total toroidal current I p and a different safety factor q 0 at the magnetic axis. To understand these degrees of freedom, consider that ECOM has just computed the eigenvector-eigenvalue solution (ψ, λ) of Eq. (3). The normalized total toroidal current I N p can be calculated from this solution according to Eq. (22). All is then needed to obtain an equilibrium with the desired total toroidal current I D p is the simple rescaling λ → (I N p /I D p )λ, which is equivalent to rescaling Ψ. From the definitions in Eq. (2), it is clear that dp/dψ and dF 2 /dψ are also scaled by λ, and must be rescaled as well: dp/dψ → (I D p /I N p )dp/dψ and dF 2 /dψ → (I D p /I N p )dF 2 /dψ. Once λ is fixed, there still is a degree of freedom for the determination of q, because q depends on F instead of dF 2 /dψ, as can be seen in Eq. (28). In ECOM, this degree of freedom can be removed by specifying the value of either the poloidal current or the safety factor at a certain radial location.
There are several options for the scaling of the normalized equilibrium in ECOM, with corresponding namelist variable ISCALE. They are summarized in Table I. If ISCALE=0 or ISCALE=1, λ is not rescaled, so that the total toroidal current is the normalized total toroidal current I N p . When ISCALE=0, the degree of freedom associated with q is removed by specifying either F (Ψ = Ψ 0 ) or F (Ψ = Ψ b ). For F (Ψ = Ψ 0 ), the additional namelist variable IFPOL needs to be set to 0, while for F (Ψ = Ψ b ) IFPOL needs to be set to 1. When ISCALE=1, q(Ψ = Ψ 0 ) is specified. If ISCALE=2 or ISCALE=3, λ is rescaled so that the total toroidal current is adjusted to the desired toroidal current I D p . The choices for constraining the q profile are the same as before: when ISCALE=2, either F (Ψ = Ψ 0 ) or F (Ψ = Ψ b ) is given, depending on the value of IFPOL, and when ISCALE=3, q(Ψ = Ψ 0 ) is given. Finally, a last option to fix λ in ECOM is to specify q(Ψ = Ψ 0 ) and either F (Ψ = Ψ 0 ) or F (Ψ = Ψ b ). This option corresponds to ISCALE=4. D. Evaluation of the parallel current density ECOM includes the option to evaluate the neoclassical parallel current density using a reduced description for the bootstrap current, either based on the Hirshman model [35] or on the Sauter model [36], and the Sauter formula for the Ohmic current [36]. The namelist variables associated with these capabilities are IBSCUR and IJBSMODEL. IBSCUR must be set to 1 for ECOM to calculate the parallel current, and the model ECOM uses for the calculation of the bootstrap current depends on the value of the variable IJBSMODEL. If IJBSMODEL=1, the Hirshman model is used, if IJBSMODEL=2, the Sauter model is used.
Consider the parallel current density in the formula The contribution of the ohmic current to the parallel current is determined by the loop voltage and the neoclassical resistivity according to where σ neo is the neoclassical resistivity, and V loop is the loop voltage. In ECOM, the value of the loop voltage is specified by the namelist variable VLOOP0, in unit of volts, and the Sauter model is implemented for the evaluation of the neoclassical resistivity σ neo : where σ spitz is the Spitzer resistivity as defined in [36], Z i is the ion charge, and f tef f (33) is given by Eq. (13b) in [36]. To evaluate f tef f (33) , the effective passing particle fraction f p (ψ) is calculated on each flux surface using the following formula where B max is the maximum value of the magnetic field on the flux surface. The integral is computed numerically with a Chebyshev-Gauss quadrature. If IJBSMODEL=1, ECOM evaluates the contribution of the bootstrap current J R B (ψ) to the parallel current with the Hirshman model [35]. The quantity J · B B is calculated using the formulae (23) to (25) in Reference [18] and the passing particle fraction f p in Eq. (36) of the present article. If IJBSMODEL=2, ECOM relies on the Sauter model [36] to compute J · B B . Specifically, where η i = (d ln n/dψ) −1 (d ln T i /dψ), the electron and ion densities are the same, n e = n i = n, and the electron and ion temperatures are also assumed to be equal for simplicity, T e = T i . L 31 ,L 32 , L 34 , and α are all defined in [36], and are evaluated using f p in Eq. (36).

E. MHD stability
If the namelist variable ISTABILITY is set to 1, ECOM verifies whether the computed equilibrium crosses or not the Troyon limit [32]. This limit is thought to ensure the no-wall stability of the equilibrium to the n = 1 internal kink mode, as well as ballooning, and external ballooning-kink modes, and can be expressed in two equivalent ways [37], where q = 2B 0 A 0 /(µ 0 R 0 I p ), A 0 = 0 1 dψI D (ψ) is the total poloidal cross section area, I D (ψ) = 2π 0 dθ(J/R), and β has been defined in Eq. (24).
In addition, ECOM also checks the Mercier criterion for stability against interchange modes [31] on each flux surface, given by where

F. Miller parametrization of the flux surfaces
If the namelist variable IFITMIL is set to the value 1, ECOM uses a nonlinear least square method to compute the Miller parametrization that best fits a given flux surface of interest. The details of the fitting method can be found in Appendix A. The outputs of the calculation are the Miller parameters κ, δ m , a/R m0 , dR m0 /dρ, dκ/dρ, dδ m /dρ, q, s, and α m , which can for example be used in ballooning stability studies and in gyrokinetic codes. α m is given by the expression [21]

IV. ACCURACY AND SPEED
In this section, we consider two examples to compare the performance of ECOM with that of the popular G-S solver CHEASE [18]. The first example corresponds to a family of equilibria originally studied by Solov'ev [38], for which simple analytic expressions can be written for the solution Ψ. These equilibria are particularly advantageous for detailed error analysis, but lack generality in the sense that the G-S equation is linear and does not have to be solved as an eigenvalue problem. In Section IV B we thus consider a more general equilibrium, with p and F profiles chosen in such a way that the G-S equation is nonlinear and has to be solved as an eigenvalue problem. Every computational test in this article is conducted using a single core 2.6GHz AMD Opteron processor with 8GB of memory.

A. Example 1: Solov'ev profiles
For the first example, we consider the Solov'ev profiles µ 0 p(Ψ) = −C s Ψ and F (Ψ) = F B , where C s and F B are constants. The G-S equation then reduces to ∆ * Ψ = C s R 2 , and an up-down symmetric solution is given by the following expression [18] where C s = F B (κ + 1/κ)/(R 3 0 q 0 ), R 0 and q 0 are the major radius and the safety factor at the magnetic axis, and a and κ are the effective minor radius and elongation of the last closed flux surface, given by Ψ = 0. A particularly convenient parametrization for the surface, which ECOM uses to compute the conformal map from the plasma boundary to the unit disk, is given by: The poloidal flux takes its minimum value Ψ 0 at the magnetic axis R = R 0 , Z = 0, with Ψ 0 given by The safety factor at a given flux surface Ψ = Ψ s can also be calculated exactly: where R max and R min are the solutions of Eq. (43) for Ψ = Ψ s and Z = 0 satisfying R max > R min > 0 and R 2 max + R 2 min = 2R 2 0 . Here, E(k) is the complete elliptic integral of the second kind with modulus k = 1 − (R min /R max ) 2 . E(0) = π/2 and the integral formula 6 in Section 3.156 of Reference [39] were used to derive Eq. (49). For simplicity, we define the magnetic shearŝ(Ψ) in terms of Ψ, giving the exact formulâ where K(k) is the complete elliptic integral of the first kind. In ECOM, this Solov'ev case is computed when the namelist variables are chosen such that IPTYPE=0, IFTYPE=0, and IBTYPE=0. Figure 2 shows the error in the L ∞ norm between the numerical values of Ψ, ∂Ψ/∂R, ∂ 2 Ψ/∂R 2 , and q calculated with ECOM and the exact values computed from Eq. (43) and Eq. (49), for the parameters R 0 = 1.0, a/R 0 = 0.32, κ = 1.7, F B = 1.0 and q 0 = 1.0. The expression "on grid" means that the error is evaluated at the (R, Z) points of the grid on which ECOM solves the G-S equation. The curves labeled "at contours", on the other hand, also include the error induced by the postprocessing steps described in Section III. Specifically, for all the flux contours constructed in the postprocessing phase, we compute the error between the value of Ψ at the contour and the actual value of Ψ at this location according to Eq. (43). The largest of these errors over the whole domain Ω is used to plot the curve we call "Ψ at contours". The curve labeled "q at contours" represents the maximum error between q at the contours as computed by ECOM according to Eq. (28) and the exact value of q at these locations as given by Eq. (49). Figure 2 demonstrates the exponential convergence of the maximum error as the number of grid points is increased, as pointed out in the introduction. The green dashed line in Figure 2 indicates that the convergence rate is approximately 1.05 −N . Figure 2 is also a proof that the derivatives of Ψ have a convergence rate that is similar to that of Ψ. Note finally that the numerical methods used in the postprocessing steps to compute contours of constant flux lead to similar convergence rates and accuracy for Ψ and q on the contours.
At equal grid size, ECOM is much faster than CHEASE, as shown in Figure 3. The run time complexity of the solver is O(n r n ϑ log n ϑ ) instead of O(n 2 r n 2 ϑ ) for typical finite element based codes solving the G-S equation. Note that the "solver" part of ECOM represented by the solid line in Figure 3 includes the run times of both the conformal mapping and the Poisson solver. The run time of the Poisson solver is the major contributor to the total run time in ECOM for typical grid sizes, because the Poisson solver is typically called 20 to 30 times during an equilibrium calculation while the conformal mapping routine is called at most a few times. The operation count of postprocessing depends on the number of contours n f , which is typically smaller than the number of radial grid points: n f = 10 − 30 < n r . The run time of postprocessing is relatively short because most computations are one dimensional, as described in Section III B.
Before comparing the accuracy of ECOM and CHEASE, it is instructive to look at the consequences of the grid crowding effect due to the conformal map on the accuracy of ECOM for shaped plasma equilibria. For fusion applications, a key question is how fast the numerical error evolves as the egg-shaped plasma cross section gets more and more elongated. Figure 4 provides an answer to that question for the Solov'ev equilibrium studied in this section, showing a significant deterioration of the accuracy of ECOM as the elongation of the last closed flux surface is increased. The consequence of this result is that for elongations corresponding to modern tokamaks and spherical tokamaks, κ 1.5 − 2, ECOM will often require a denser grid than FEM based G-S solvers to achieve the same accuracy. Since ECOM is much faster than these solvers in terms of work per grid point, and computes derivatives with high accuracy, ECOM remains very often more desirable that FEM based solvers, as we will show next. However, ECOM is not an attractive option to compute highly elongated equilibria, such as those in Field Reversed Configurations (FRCs) [40], with κ ∼ 10.
We have just shown that ECOM is much faster than CHEASE for a given grid size, but that for elongated plasma shapes ECOM may need a denser mesh to achieve a desired accuracy, due to crowding effects. In this context, a fair  comparison of the performance of the two codes is done by evaluating the accuracy of each solver for a given run time on the same machine. This is precisely the purpose of Figure 5, which shows the numerical error in the safety factor q and the magnetic shearŝ as a function of the normalized radius ρ = (Ψ − Ψ 0 )/(Ψ b − Ψ 0 ) for the Solov'ev equilibrium considered in this section and κ = 1.7. We have chosen to focus on q andŝ because MHD stability and turbulent transport are known to depend sensitively on these quantities. We compare ECOM (red markers) and CHEASE (blue markers) for three different run times. In Figure 5-(a), the run time is 1 second, corresponding to a grid resolution of n r = 32 in CHEASE and n r = 96 in ECOM; in Figure 5-(b), the run time is 3 seconds, corresponding to a grid size of n r = 52 in CHEASE and n r = 192 in ECOM; in Figure 5-(c), the run time is 5 seconds for a grid resolution of n r = 64 in CHEASE and n r = 256 in ECOM. One can see that for small grids, CHEASE computes the safety factor with a better accuracy than ECOM, a direct consequence of the crowding effect. However, even if in that case CHEASE calculates q more accurately, the accuracy for the magnetic shear are comparable in ECOM and CHEASE. One reason for this is that in ECOM, we constructed the n f flux contours so that they would coincide with a global Chebyshev grid of size n f on the interval ψ = 0 and ψ = 1. ECOM can thus use spectral differentiation to computeŝ, leading to a more limited loss of accuracy between q andŝ. The construction of such a Chebyshev grid for the flux variable is particularly convenient in ECOM because the piecewise Chebyshev grid used by the Poisson solver to discretize the radial direction in D 1 is well refined near the end points ψ = 0 and ψ = 1 of the interval. Radial derivatives of flux functions (e.g.ŝ, dV /dρ and dψ/dρ) are therefore calculated without significant loss of accuracy in ECOM.
As the grid size and computation time are increased, ECOM outperforms CHEASE, which is a direct result of the geometric convergence demonstrated in Figure 2. For a run time of 3 seconds, q is computed with similar accuracy in ECOM and CHEASE, but the error onŝ is more than 100 times smaller in ECOM. For a run time of 5 seconds, the error on both q andŝ with ECOM is orders of magnitude smaller than the error obtained with CHEASE.
Although memory aspects rarely lead to severe constraints on G-S solvers in fusion applications, it is interesting to note that memory requirements are much smaller in ECOM than they are in CHEASE. As an empirical illustration of this, we observed that because of the memory limitations of the computer we used for our comparison studies, we were limited to grids smaller than n r n ϑ 2 × 10 4 in CHEASE, whereas we computed equilibria in ECOM with grids n r n ϑ ≥ 10 6 without any difficulty. For the second example, we choose pressure and current profiles in such a way that the G-S equation is an eigenvalue partial differential equation given by Eq. (3). Specifically, we set the namelist variables IPTYPE and IFTYPE to 1, and specify the profile constants according to p in = 2, p out = 1, F 0ψ = −1, F in = 2, and F out = 1. In the three equilibria we study in this section, we will vary p 0ψ : for the first equilibrium, we set p 0ψ = −0.01/µ 0 , leading to a very small Shafranov shift, and we set p 0ψ = −1/µ 0 for the next two equilibria, which consequently have a much larger Shafranov shift. To describe the last closed flux surface of the equilibrium, we set the namelist varibale IBTYPE to 1, corresponding to the Miller parametrization: where the parameter t goes from 0 to 2π, κ is the elongation and δ m is the triangularity. For the purpose of comparison, we specify the same profiles and parametrization of the plasma boundary in the input file of CHEASE, and for both codes we look at the convergence of the poloidal magnetic field energy W p , given by Eq. (27). This global, 0-D quantity has often been used to measure the convergence properties of CHEASE [41]. For the equilibrium under consideration, the analytic expression for ψ is not known, and there does not exist a formula for W p . We thus use the value of W p (n ref ) calculated with a large number of grid points as the reference point for the convergence studies. Figures 6 -8 show the convergence of W p in CHEASE and in ECOM as a function of the number of grid points and as a function of run time, for three different plasma shapes. Figure 6 corresponds to a circular tokamak equilibrium with a small Shafranov shift. Since the domain Ω is a disk, ECOM computes such equilibria without conformal mapping, and this case allows us to focus on the error that is not due to the conformal mapping part of the solver. We observe that for a very small number of grid points, N ≤ 16, CHEASE computes W p with higher accuracy than ECOM. This is somewhat artificial in that CHEASE uses a grid that is refined near the magnetic axis [18], whereas in the absence of conformal mapping the center of the (r, θ) grid in ECOM does not coincide with the magnetic axis. The convergence rate of the poloidal magnetic energy in CHEASE is N −6 as found in [41], while the convergence rate in ECOM is found to be 2.6 −N . As a result, for grids with N > 16 ECOM quickly becomes much more accurate than CHEASE. Finally, if we fix the run time instead of the grid size, we find that ECOM is always more accurate than CHEASE, as shown in Figure 6-(b).
In Figure 7, we consider an equilibrium with a larger Shafranov shift, a significant triangularity but with no elongation. Comparing Figure 6 and Figure 7, we can see that the crowding effect inherent to conformal mapping results in a strong loss of accuracy, with a relative error in W p which is up to 10 4 times larger at low grid resolutions. It also results in a reduction of the convergence rate from 2.6 −N to 1.23 −N . Even if so, ECOM remains more accurate than CHEASE at fixed run time, as shown in Figure 7-(b). Figure 8 corresponds to an ITER-like equilibrium [42], with the same triangularity and pressure profile as Figure 7, but with elongation κ = 1.7. As we would expect from Figure 4, elongation amplifies the crowding effect, leading to further degradation of the accuracy. The convergence rate is reduced from 1.23 −N to 1.05 −N . The strong dependency of the convergence rate on the plasma geometry and on crowding is an undesirable aspect of ECOM. In contrast, Figures 6-8 show that convergence in CHEASE is fairly insensitive of the shape of the plasma boundary. Despite this weakness, ECOM computes W p more accurately than CHEASE for run times longer than 8 seconds. Furthermore, the results in Section IV A suggest that the run time threshold is lower for local quantities, in particular if these quantities depend on derivatives of flux functions, such aŝ s or high order derivatives of the ψ, such as the local magnetic shear.

V. EQUILIBRIA WITH TOROIDAL FLOWS
Large equilibrium flows are observed in tokamak experiments [43][44][45], and flows and flow shear are thought to have a strong influence on the stability and transport properties of the plasma [44][45][46][47][48][49][50][51][52][53]. When the flow speed is of the same order as the sound speed, the inertial term in the pressure balance relation can no longer be ignored, and MHD equilibrium force balance is given by where ρ = m i n, m i is the ion mass, n the ion density, and u the plasma flow. The plasma flow u in Eq. (54) must also satisfy the steady-state version of Faraday's law in the ideal MHD model [9]: Poloidal flows are damped by neoclassical viscosity and expected to be much smaller than the ion sound speed in toroidally axisymmetric equilibria, except perhaps near the edge [54][55][56]. It is therefore a good approximation to only retain the effect of toroidal flows in Eq. (54). The most general toroidal flow u = ue φ satisfying Eq. (55) can be written as u φ = RΩ φ (Ψ). For a purely toroidal flow, MHD force balance thus takes the form where e R = ∇R. In a general axisymmetric geometry, p only depends on two variables. Since we know that for static equilibria p is a function of Ψ only, we choose R and Ψ as the two independent variables for the pressure profile in axisymmetric equilibria. We then have ∇p = ∂p/∂R∇R + ∂p/∂Ψ∇Ψ, and dotting Eq. (56) with B yields an equation for the R dependence of the pressure: where we have used the fact that B · ∇Ψ = 0. When Eq. (57) is satisfied, Eq. (56) can be written as By dotting this equation with J, it is easy to show that RB φ = F (Ψ) as in the static case, and Equation (58) becomes the following modified G-S equation for the flux function Ψ in the presence of a toroidal flow: There are two well known situations for which Eq. (57) can be integrated analytically. The first situation corresponds to the assumption that the entropy S ≡ pρ −γ is only a function of the poloidal flux [9], the second situation corresponds to the assumption that the temperature is a flux function because of the high thermal conductivity along the magnetic field lines in fusion grade plasmas [7,57,58]. Currently, ECOM only treats the latter case. Neglecting temperature anisotropy, we write p(R, Ψ) = 2n(R, Ψ)T (Ψ), with T (Ψ) = 0.5[T i (Ψ) + T e (Ψ)] a species averaged temperature, and integrate Eq. (57) to find [7]: where p 0 (Ψ) = p(R 0 , Ψ) and p Ω (Ψ) = n 0 (R 0 , Ψ)m i Ω φ (Ψ) 2 R 2 0 /2 is the kinetic pressure due to the ion toroidal flow. ECOM uses the numerical scheme described in Section II to solve the normalized version of Eq. (59), with the normalized pressure term given bȳ In ECOM, the profilep 0 (ψ) is specified in the same way asp(ψ) is in the static case, with the same namelist variable IPTYPE, and the same options. If IPTYPE=1 or IPTYPE=2 and IPTABLE=0, integration is required to obtain p 0 (ψ) from its flux derivative. ECOM uses Chebyshev-Gauss quadrature on the global Chebyshev grid for ψ to compute these integrals. There are several options to specify the kinetic pressure profilep Ω (ψ), with corresponding namelist parameter ITFTYPE. If ITFTYPE=1,p Ω is such that the toroidal Mach number M = 2p Ω /p 0 has the same value at all radii. If ITFTYPE=2,p Ω is given by the explicit formulap Ω (ψ) =p Ω0 (1 − (1 − ψ) pΩin ) pΩout . If ITFTYPE=3,p Ω is given as a numerical table in terms of ψ, as is also done for the pressure and poloidal current profiles.
Among the expressions ECOM uses in postprocessing to evaluate flux functions and figures of merit, only a few need to be modified in the presence of an equilibrium toroidal flow. Equation (20) becomes where dp ∆ /dψ = dp Ω /dψ − (p Ω /p 0 )dp 0 /dψ, and I B1 (ψ) and I B2 (ψ) are defined by As a result, the intermediate step we use to derive Eq. (22) takes a slightly different form, but Eq. (22) itself does not change: the relation I φ (ψ) = I A (ψ)/(λµ 0 ) still holds. Finally, the expression for the volume averaged pressure is now given by Figure 9 shows the flux contours of a stationary equilibrium with toroidal flow computed with ECOM and the flux contours of the corresponding static equilibrium also computed with ECOM. For that example, we chose p Ω so that the toroidal Mach number is uniform with value 1. The flux contours of the stationary equilibrium are the red lines, and the flux contours of the static equilibrium are the black dashed lines. We can clearly see the expected outward shift of the magnetix axis [58,59].

VI. DISCUSSION
ECOM uses conformal mapping from the plasma domain to the unit disk to decouple the numerical issues associated with the plasma geometry from the rest of the problem. Once on the unit disk, ECOM relies on fairly standard high order methods based on the FFT for the angular dependence and a Green's functions formulation for the radial dependence to solve the mapped partial differential equation describing the plasma equilibrium. This approach is not only conceptually elegant, it is also effective. By studying a static equilibrium with Solov'ev profiles, we showed in this article that the solution of the Grad-Shafranov equation as well as its first and second derivatives converge exponentially as grid size is increased. Furthermore, ECOM is much faster than finite element based codes in terms of work per grid point, and requires less memory at equal grid size. Finally, ECOM can be easily parallelized in multi-core system since the solver on the unit disk uses separation of variables and solves an independent radial ODE for each angular grid point.
Because of the crowding effect inherent to the mapping of an elongated shape to a disk, ECOM usually requires a denser grid than FEM solvers to achieve the same accuracy. For tokamak geometries, this weakness is compensated by the speed of the solver: beyond a threshold grid size, ECOM computes the equilibrium quantities that play a key role in wave propagation, stability and transport calculations with more accuracy than CHEASE at equal run time. The threshold grid size depends on the plasma geometry and on the quantity that is calculated. For an ITER-like geometry and quantities that depend on first derivatives of the flux, such as the safety factor and the poloidal magnetic field energy, the threshold grid size is n r = n ϑ 360. For plasma shapes that are less elongated, the threshold grid size is smaller. It is also smaller if the quantity of interest depends on second derivatives of the flux, such as the local magnetic shear for instance, as a direct consequence of the exponential convergence of the second derivatives in ECOM. We also find the threshold grid to be smaller when the quantity of interest is a flux derivative of a flux function, such as the flux averaged magnetic shear for instance. This is because ECOM uses a Chebyshev grid for the flux contours and spectral differentiation for the evaluation of flux derivatives.
ECOM has two important limitations. First, equilibria with κ > 3 require too dense a grid for ECOM to reach a high level of accuracy in a reasonable amount of computing time. In these situations, particularly relevant to FRCs, ECOM is not an attractive option. Second, ECOM can only compute equilibria whose boundaries are smooth. It can therefore not be used for equilibria with a magnetic X-point. A promising idea to address these limitations is to develop a Grad-Shafranov solver based on an integral equation formulation that avoids conformal mapping to treat the geometrical aspects of the problem. Approaches relying on the Fast Multipole Method [60,61] may represent an attractive option, that would lead to high order accuracy for the solution of the G-S equation as well as its derivatives. They are the subject of ongoing research.
Note that there are additional options in ECOM that can be very desirable for certain applications, but are not discussed in this article. For instance, ECOM can handle equilibria that are specified in terms of the parallel current J (ψ) or the safety factor q(ψ) instead of the poloidal current F (ψ). Any one of the three profiles can be used in ECOM along with the specification of the pressure profile. When either the J profile or the q profile constrains the equilibrium, ECOM needs to evaluate the flux functions I A , I B , and I C in Eq. (18) at each iteration. The fast and accurate numerical methods implemented in ECOM to calculate these functions then become a key strength of the solver, leading to fast convergence of the iterations and accurate equilibria. ECOM can also compute equilibria specified by an EFIT g-file [62] containing the pressure profile, the poloidal current profile, and the boundary shape of interest. And conversely, ECOM can print the results of any equilibrium calculation according to the format of an EFIT g-file, which includes the pressure, poloidal current and q profiles, as well as the boundary and ψ(R E , Z E ) where (R E , Z E ) is a uniform grid in the range of R min < R E < R max and Z min < Z E < Z max . For the sake of clarity and conciseness of the presentation, we did not describe these capabilities in the present article. They will be presented in detail in forthcoming articles, in which we explicitly use them to explore properties of tokamak equilibria and to couple ECOM with wave propagation and transport codes.
The entries of the Hessian matrix are , ECOM starts the iterative procedure with the following initial guesses Appendix B: Summary of namelist variables in ECOM Profile specified along with the pressure profile in the G-S equation 0 (default) Poloidal current profile dF 2 /dψ 1 Parallel current profile J 2 Safety factor profile q IPTYPE Specification of the pressure profile 0 µ0dp/dψ = −Cs for the Solov'ev solution given by Eq. (43) Discrete values of dp/dψ or p in terms of ψ or ρ are given by a table in 'file prof ' 3 Discrete values of dp/dψ in terms of ψ is given by the EFIT output 'file efit' Cs = FB(κ + 1/κ)/(R 3 0 q 3 0 ) is determined by the namelist variables 'F0','q0','rkappa', and 'R0' The namelist variables for p 0ψ ,pin and pout are 'p0psi', 'pin' and ,'pout', respectively IFTYPE Specification of the poloidal current profile (activated for IECOM=0) 0 dF 2 /dψ = 0 for the Solov'ev solution given by Eq. (43)  1(default) ( Discrete values of (1/2)(dF 2 /dψ) or F in terms of ψ or ρ are given by a table in 'file prof ' 3 Discrete values of (1/2)(dF 2 /dψ) in terms of ψ given by the EFIT output 'file efit' The namelist variables for F 0ψ ,fin and fout are 'ff0','ffin' and ,'ffout', respectively IJTYPE Specification of the parallel current profile (activated for IECOM=1) 1(default) Discrete values of J in terms of ψ are given by a table in 'file jprof' 3 J as evaluated from ohmic and bootstrap current models at each iteration the namelist variables for J 0 ,jin and jout are 'jpar0','jpin' and ,'jpout', respectively IQTYPE Specification of the q profile (activated for IECOM=2) 1(default) Discrete values of q in terms of ψ are given by a table in 'file qprof' The namelist variables for q0,q f ac qin and qout are 'q0','qfac','qpin' and ,'qpout', respectively IBTYPE Specification of the plasma boundary ∂Ω 0 ∂Ω is given by Eqs. Discrete values ofpΩ in terms of ψ are given by a table in 'file tflow' The namelist variables for M ,pΩ0,pΩin and pΩout are 'mach','ptf0','ptfin' and ,'ptfout', respectively IPTABLE Type of pressure profile table in 'file prof' (activated for IPTYPE=2) 0(default) 1-D arrays of ψ and dp/dψ are given 1 1-D arrays of ψ and p are given 2 1-D arrays of ρ and p are given IRHO Definition for the normalized radius ρ 0(default) ρ(ψ) = (Ro(ψ) − R0)/(Ro(ψ = 0) − R0) 1 ρ(ψ) = (Ro(ψ) − Ri(ψ))/(Ro(ψ = 0) − Ri(ψ = 0)) 2 ρ(ψ) = (Ψ − Ψ0)/(Ψ b − Ψ0). nt1 Number of grid points on the boundary used for the forward conformal mapping (n1) nt2 Number of ϑ grid points in the unit disk (n ϑ ) nt3 Number of θ grid points for contour integrals (n θD = n θE ) nsub Number of radial piecewise Chebyshev intervals in the unit disk (nL) kcheb Number of Chebyshev points in a radial interval (n ch ) nchq Number of flux surfaces in a Chebyshev grid of ψ for postprocessing (n f ) nflx Number of flux surfaces in uniform grid of ρ for postprocessing kLag Order of Lagrange interpolation (kLag) ksamp Oversampling factor for FFT padding used for contour integrals (ksamp) R0 R coordinate of the point that is mapped to the center of D1 by the initial conformal mapping R0 for IBTYPE=0 and Rm0 for IBTYPE=1 Z0 Z coordinate of the point that is mapped to the center of D1 by the initial conformal mapping Z0 for IBTYPE=0 and Zm0 for IBTYPE=1 q0 Value of the safety factor at the magnetic axis q(Ψ = Ψ0) F0 F (Ψ = Ψ0) for IFPOL=0 or F (Ψ = ΨB) for IFPOL=1 ISCALE Control parameter to scale the solutions. See Table I  Printing an Ascii file of Ψ and its derivatives on the grid of the unit disk D1 (0: off, 1: on) IPRINTCON Printing Ascii files of contours of nchq ψ values and nflx ρ values (0: off, 1: on) IPRINTQS Printing an Ascii file of the safety factor and the magnetic shear in terms of nchq ψ and nflx ρ (0: off, 1: on) IPRINTEFIT Printing an Ascii file in EFIT g file format (0: off, 1: on) npsi Number of flux surfaces for EFIT g file format for iprintefit=1 ISTABILITY Evaluation of the Mercier criterion and Troyon limit (0: off, 1: on) IBSCUR Evaluation of the bootstrap and ohmic currents (0: off, 1: on) IJBSMODEL Bootstrap and ohmic current model (1: Hirshman model, 2: Sauter model) nchy Number of pitch angle grid points to evaluate the bootstrap and ohmic currents VLOOP0 Loop voltage giving the ohmic current in unit of volt IFITMIL Fitting the flux surfaces using the Miller parametrization (0: off, 1: on)