CREATE-FXB, a Fixed Boundary Code Based on Finite Element Methods for the Solution of Grad-Shafranov Equation and Optimization of Equilibrium Currents

1. Introduction

Numerical models solving the Grad-Shafranov equation in magnetically confined toroidal plasmas are crucial both for the design phase of plasma scenarios [1] and for the development of active magnetic control systems [2]. For these applications, a variety of procedures and associated codes for the calculation of Magnetohydrodynamic (MHD) equilibria are widely used [3,4,5,6,7,8], offering an adequate trade-off between computational efficiency and precision.

In tokamaks, a fixed-boundary equilibrium problem assumes a known plasma shape as input, whereas a free-boundary equilibrium problem simultaneously provides both the plasma shape and internal equilibrium conditions as results. Fixed-boundary problems are computationally simpler and faster, while free-boundary problems are typically employed to test plasma control system performance and analyse experimental data.

The present work focuses on fixed-boundary problems. This approach reduces computational effort, as only the plasma region is meshed, i.e., where the Grad-Shafranov equation needs to be solved. Furthermore, it enables the determination of the plasma response to external poloidal fields without explicitly defining the external coil system. Consequently, it can be used to obtain linearized plasma response to any geometry and number of poloidal field coils, improving flexibility.

This paper presents a new tool called CREATE-FXB, a Finite Element Method (FEM) solver for fixed-boundary problems. The code incorporates several features that make it particularly suited for tokamak applications, including the design of plasma scenarios, high-accuracy derivation of linearized plasma responses, treatment of a wide class of current density profiles (including non-smooth distributions), and improved interoperability with other existing codes. CREATE-FXB employs an innovative nodal numerical formulation, which allows for the treatment of generic current density profiles provided by transport codes [9,10], including profiles with points of non-regularity, in both fixed-boundary problem resolution and subsequent linearization. An optimization tool is also included for estimating the equilibrium currents once the geometry of the magnetic system is fixed.

The paper is organized as follows. Section 2 describes the main assumptions. Section 3 presents the Grad-Shafranov equilibrium equation and its numerical resolution with the formulation used by CREATE-FXB. Section 4 illustrates some applications of the procedure. Finally, Section 5 draws the conclusions and discusses potential future developments compatible with the model.

2. Assumptions

2.1. Basic Equations

At a generic point of the plasma region ${\Omega }_{pl}$ , static (zero plasma velocity v=0) and stationary (∂/∂t=0) plasma equilibria are well described by the following basic equations below, obtained by means of simplifications of the more complex MHD model [11]:

$$\nabla \times \mathit{B}={\mu }_0\mathit{J}$$

(1)

$$\nabla \bullet \mathit{B}=0$$

(2)

$$\mathit{J}\times \mathit{B}=\nabla p$$

(3)

where $\mathit{J}$ is the current density, $\mathit{B}$ is the magnetic field, ${\mu }_0$is the vacuum magnetic permeability assumed in the plasma, and p is the kinetic pressure.

Under the 2D axisymmetric assumption in a cylindrical reference coordinate system $\left(r,\varphi ,z\right)$, using Ampere’s law (1) and the divergence-free condition for the magnetic field (2), both $\mathit{J}$ and $\mathit{B}$ can be expressed in terms of two scalar functions, the magnetic poloidal flux per radian function $\psi \left(r,z\right)$ and poloidal current function $f\left(r,z\right)$:

$$\mathit{B}={\displaystyle \frac{\nabla \psi \times {\mathit{i}}_{\mathit{\varphi }}}{r}}+{\displaystyle \frac{f}{r}}{\mathit{i}}_{\mathit{\varphi }}$$

(4)

$$\mathit{J}={\displaystyle \frac{\nabla \left(\frac{f}{\mu }\right)\times {\mathit{i}}_{\mathit{\varphi }}}{r}}+L\psi {\mathit{i}}_{\mathit{\varphi }}$$

(5)

where $L$ is a 2^nd order elliptic operator defined by:

$$L\psi =-{\displaystyle \frac{\partial }{\partial r}}\left({\displaystyle \frac{1}{{\mu }_{0}r}}{\displaystyle \frac{\partial \psi }{\partial r}}\right)-{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{1}{{\mu }_{0}r}}{\displaystyle \frac{\partial \psi }{\partial z}}\right)$$

(6)

Referring to the plasma region ${\Omega }_{pl}$, by substituting (4) and (5) into the (3) and noting that (5) implies:

$$L\psi =J_{\varphi }$$

(7)

the Grad-Shafranov (G-S) equation is derived:

$$L\psi =r{\displaystyle \frac{dp}{d\psi }}+{\displaystyle \frac{1}{{\mu }_{0}r}}{\displaystyle \frac{d\left(\frac{f^2}{2}\right)}{d\psi }}$$

(8)

The functions ${\displaystyle \frac{dp}{d\psi }}$ and ${\displaystyle \frac{d\left(\frac{f^2}{2}\right)}{d\psi }}$are named free functions.

In the fixed boundary problems, the differential problem is closed by coupling the PDE to some Dirichlet boundary conditions at the plasma boundary $\partial {\Omega }_{pl}$:

$$\psi \left(r,z\right)\left.\right|\partial {\Omega }_{pl}={\psi }_b, p\left(r,z\right)\left.\right|\partial {\Omega }_{pl}=p_b, f\left(r,z\right)\left.\right|\partial {\Omega }_{pl}=f_b$$

(9)

The problem resolution requires the expression of $J_{\varphi }$. Thus, the free functions must be assigned or derived in some way.

2.2. Current Density Modelling

The current density modelling passes through the derivation of the free functions $dp/d\psi$ and $d\left(f^2/2\right)/d\psi$. Usually, they are given as function of the non-dimensional variable $\overline{\psi }\in \left[\mathrm{0,1}\right]$ defined as follows:

$$\overline{\psi }={\displaystyle \frac{\psi -{\psi }_a}{{\psi }_b-{\psi }_a}}$$

(10)

where ${\psi }_a$ and ${\psi }_b$ are the values of the flux per radian at the magnetic axis and at the plasma boundary, respectively.

These functions can be obtained by coupling a simplified one-dimensional version of the Grad-Shafranov equilibrium to the transport equations [9,10], adopting iterative schemes to solve the non-linear system and derive the plasma pressure $p\left(\psi ,t\right)$ and poloidal current function $f\left(\psi ,t\right)$. Afterwards, these profiles are given to specific 2D axisymmetric codes coupled to circuit equations like CREATE-NL+ [5] to get the time evolution of the poloidal flux map.

An alternative procedure takes $dp/d\psi$ and $d\left(f^2/2\right)/d\psi$ from suitable family of functions $g\left(\overline{\psi },\mathit{\alpha }\right)$. The $\mathit{\alpha }$ variable is a set of shape parameters necessary to model the features of the current density profile as shown in [12]. They represent the degrees of freedom of the parameterization $g$, linked to the main plasma quantities, namely plasma current $I_{pl}$, poloidal beta ${\beta }_{pol}$, and internal inductance $l_i$, which can be defined as follows:

$$I_{pl}={\int }_{ {\Omega }_{pl}} J_{\varphi }d\Omega , {\beta }_{pol}=4{\int }_{ V_{pl}}{\displaystyle \frac{p}{{\mu }_0{R}_0{I_{pl}}^2}}dV, l_i=4{\int }_{ V_{pl}}{\displaystyle \frac{B_{pol}^2}{{{2\mu }_0}^2{R}_0{I_{pl}}^2}}dV$$

(11)

where ${\mathit{B}}_{pol}$is the magnetic field component in the poloidal plane $rz$.

However, different definitions of ${\beta }_{pol}$ and $l_i$ can also be employed [13]. The choice of this parameterization is motivated by the fact that the flux surfaces outside the plasma boundary are nearly identical for two equilibria with different current density profiles but the same values of $I_{pl}$, ${\beta }_{pol}$ and $l_i$ [14]. Therefore, the magnetic coupling of the plasma with the external conductors can be modelled by means of a suitable choice of the $\mathit{\alpha }$ parameters.

These $\mathit{\alpha }$ parameters are calculated using iterative methods, i.e., by solving (7)-(9) until the desired values of plasma quantities (11) are obtained. The computational effort depends on the number of degrees of freedom (DoFs). Therefore, the choice of $g$ should be made so as to ensure coverage of the widest possible range of plasma quantities with a limited amount of computational effort.

Several methods have been proposed to model the current density profile in terms of $\overline{\psi }$. For instance, power functions can be employed for both $dp/d\psi$ and $d\left(f^2/2\right)/d\psi$, as in the CREATE-L code [5], which usually implements a formulation featuring three DoFs:

$$J_{\varphi }=\lambda \left[\beta {\displaystyle \frac{r}{R_0}}+\left(1-\beta \right){\displaystyle \frac{R_0}{r}}\right]{\left(1-\overline{\psi }\right)}^{{\alpha }_m}$$

(12)

which corresponds to bell-shaped current density profiles due to the term ${\left(1-\overline{\psi }\right)}^{{\alpha }_m}.$This bell-shaped expression is relatively simple, but it cannot be used for plasma equilibria characterized by low values of $l_i$, where the current density typically exhibits a hollow profile and tends to concentrate near the boundary. A new family of functions must be determined to be non-monotonic with $\overline{\psi }$for specific values of α correspondent to a low internal inductance. Moreover, it should be as simple as possible, to guarantee convergence with sufficient rate.

The CREATE-FXB code can also use the following parameterization in terms of the four parameters ${[\alpha }_0$,${\beta }_0, {\gamma }_0, {\delta }_0]$ derived by modifying the McCarthy class of profiles [15]:

$$J_{\phi }=\left[r{\beta }_0\left({\displaystyle \frac{1}{{\mu }_0}}+{\displaystyle \frac{{\gamma }_0}{{\mu }_0}}\overline{\psi }\right)+{\displaystyle \frac{\left(1-{\beta }_0\right)}{{\mu }_{0}r}}\left({\alpha }_0\overline{\psi }+1\right)\right]{\psi }_f, {\psi }_f=\left\{\begin{array}{l}\begin{array}{c}1 \overline{\psi }\le {\delta }_0 \end{array}\\ \begin{array}{c} \\ 1-{\displaystyle \frac{{\left(\overline{\psi }-{\delta }_0\right)}^2}{{\left(1-{\delta }_0\right)}^2}} \overline{\psi }>{\delta }_0\end{array}\end{array}\right.$$

(13)

which enables the treatment of a wide class of plasmas, as shown in Section 4. It is worth noticing that the filter function ${\psi }_f$ plays a key role in shaping the plasma current density.

In some cases, the profiles both $dp/d\psi$ and $d\left(f^2/2\right)/d\psi$ provided by transport codes may be highly rough and irregular. Examples of treatment of this kind of profiles are also shown in Section 4, where an innovative technique based on information from the mesh nodes has been adopted.

All these features represent a challenge more for the linearization procedure than for the numerical resolution of the fixed boundary problem: indeed, the coupling between standard FEM with linear triangular elements with iterative fixed-point scheme (Picard iterations) is quite robust to face the non-linearity in the computation of the plasma equilibrium quantities. On the other hand, the derivation of the linearised model becomes challenging, because the linearization involves the numerical calculation of profiles derivatives, which cannot be rigorously evaluated in the non-regularity points. For these reasons, adequate meshing processes and suitable numerical techniques to treat the irregularities of the profile are crucial for derivation of the linear model. In this context, the so-called nodal formulation, a numerical technique formalized in the subsequent paragraphs, which represents a suitable way to overcome the limitations of conventional methods. Additionally, a procedure of adaptive refinement of mesh can also be used.

3. Finite Element Method Formulation in CREATE-FXB

3.1. Picard Iterations

Equation (8) is a non-linear elliptical partial differential equation, perfectly suited for the resolution using the standard FEM. The first step is to convert it into the weak formulation, by multiplying all the members by a weight function $w\left(r,z\right)\in H_0^1\left({\mathfrak{R}}^2\right)$ and applying integration by parts:

$${\int }_{{\Omega }_{pl}}{\displaystyle \frac{\nabla w\bullet \nabla \psi }{{\mu }_{0}r}} d\Omega -{\int }_{\partial {\Omega }_{pl}}{\displaystyle \frac{w}{{\mu }_{0}r}}\bullet {\displaystyle \frac{\partial \psi }{\partial \mathit{n}}} dS={\int }_{{\Omega }_{pl}}w\bullet J_{\varphi }\left(r,\psi \right) d\Omega$$

(14)

The non-linearity appears in the right-hand side, since the toroidal current density $J_{\varphi }$ depends on the unknown $\psi$.

The plasma domain is discretized into a set of triangular elements with N nodes, corresponding to the vertices of the triangles. The unknown solution $\psi \left(r,z\right)$is approximated via a linear combination of shape functions ${\varphi }_l\left(r,z\right), l=1,\dots ,N$:

$$\psi \left(r,z\right)={\sum }_{l=1}^N{\psi }_l{\varphi }_l\left(r,z\right)$$

(15)

where ${\varphi }_l\left(r,z\right)=1$at the $l$-th node, is zero at the other nodes, and is piecewise linear within the elements containing the l-th node.

The unknown variable can thus be represented as a numerical vector of length N rather than as a function, so that $N$ independent equations are required to find the solution. These are obtained from weak formulation (14) by applying Galerkin’s approach, which selects $w_l\equiv {\varphi }_l, \forall l=1,\dots ,N$:

$${\sum }_{m=1}^N\left({\int }_{{\Omega }_{pl}}\nabla {\varphi }_l\bullet {\displaystyle \frac{1}{{\mu }_{0}r}}\nabla {\varphi }_m d\Omega \right)d\Omega ={\int }_{{\Omega }_{pl}}{\varphi }_l\bullet J_{\varphi }\left(r,\psi \right) d\Omega +{\int }_{\partial {\Omega }_{pl}}{\displaystyle \frac{{\varphi }_l}{{\mu }_{0}r}}\bullet {\displaystyle \frac{\partial \psi }{\partial \mathit{n}}} dS$$

(16)

which, in a compact form, become the matrix relation:

$$K\mathit{\psi }=\mathit{F}\left(\mathit{\psi }\right)+\mathit{r}$$

(17)

where:

$$K_{lm}={\int }_{{\Omega }_{pl}}{\displaystyle \frac{\nabla {\varphi }_l\bullet \nabla {\varphi }_m}{{\mu }_{0}r}} d\Omega , F_l={\int }_{{\Omega }_{pl}}{\varphi }_l\bullet J_{\varphi }\left(r,\psi \right) d\Omega , r_l={\int }_{\partial {\Omega }_{pl}}{\displaystyle \frac{{\varphi }_l}{{\mu }_{0}r}}\bullet {\displaystyle \frac{\partial \psi }{\partial \mathit{n}}} dS$$

(18)

At this point, Dirichlet’s boundary conditions must be imposed on (17), according to (9). In this case the last term $r_l$, which would be needed in case of Neumann boundary conditions, becomes redundant. To this end, the set of nodes is partitioned into two subsets: ${\mathit{N}}_i$ corresponding to internal nodes, and $N_b$ of internal and boundary nodes, marked by indices ii and ib, respectively. The same operation is done for the nodal values of $\psi$. In system (17) K is then rewritten in block format:

$$\left[\begin{array}{cc}{K}_{ii,ii}& K_{ii,ib}\\ K_{ib,ii}& K_{ib,ib}\end{array}\right]\left[\begin{array}{c}{\mathit{\psi }}_{\mathit{i}\mathit{i}}\\ {\mathit{\psi }}_{\mathit{i}\mathit{b}}\end{array}\right]=\left[\begin{array}{c}{\mathit{F}}_{\mathit{i}\mathit{i}}\\ {\mathit{F}}_{\mathit{i}\mathit{b}}\end{array}\right]+\left[\begin{array}{c}{\mathit{r}}_{\mathit{i}\mathit{i}}\\ {\mathit{r}}_{\mathit{i}\mathit{b}}\end{array}\right]$$

(19)

where clearly ${\mathit{r}}_{\mathit{i}\mathit{i}}=0$ due to properties of ${\varphi }_l$ functions, ${\mathit{r}}_{\mathit{i}\mathit{b}}$ can be ignored since only the first subset of $ii$ equations is required for the next step, ${\mathit{\psi }}_{\mathit{i}\mathit{b}}$ is known (Dirichlet’s boundary condition) and $\mathit{F}$ depends on $\mathit{\psi }$.

Therefore, a Picard scheme is adopted to solve the non-linearity, with the unknown vector ${\mathit{\psi }}_{\mathit{i}\mathit{i}}$ uniquely determined by solving, at each step $k$:

$${\mathit{\psi }}_{\mathit{i}\mathit{i}}^{\left(\mathit{k}\right)}=K_{ii,ii}^{-1}\bullet \left({\mathit{F}}_{\mathit{i}\mathit{i}}\left({\mathit{\psi }}^{(k-1)}\right)-K_{ii,ib}\bullet {\mathit{\psi }}_{\mathit{i}\mathit{b}}\right)$$

(20)

until a convergence criterion is satisfied.

Additional information on the well-posedness of system (17) can be found in [16].

3.2. Nodal Formulation for $J_{\varphi }$

A critical aspect in solving the Grad–Shafranov equation is the accurate evaluation of the current density $J_{\varphi }$, especially in the presence of steep gradients or discontinuities in the equilibrium profiles. Standard approaches use linear triangular elements and Gauss integration, e.g., in the CREATE-NL+ code [5], where the term ${\mathit{F}}_{\mathit{i}\mathit{i}}\left({\mathit{\psi }}^{(k-1)}\right)$ in (18) is computed iteratively from the values of $J_{\varphi }(r_g^{\left(e\right)},z_g^{\left(e\right)})$ and the shape functions ${\varphi }_l(r_g^{\left(e\right)},z_g^{\left(e\right)})$ in each triangle barycentre $e=1,\dots ,N_e$ as well as the $\overline{{\psi }_g}$ map is determined in these points, while the free function values are read from the diagrams ${\displaystyle \frac{dp}{d\psi }}\left(\overline{\psi }\right)$ and $f{\displaystyle \frac{df}{d\psi }}\left(\overline{\psi }\right)$. However, the main issue arises when the profiles present discontinuities. In such cases, if the $\overline{\psi }$ level at which $J_{\varphi }$ is discontinuous falls into a bend of elements, the integrals are done with a wrong value of $J_{\varphi }\left({\overline{\psi }}_g\right)$. To address this limitation, a nodal formulation is introduced, in which the current density is reformulated in terms of nodal quantities, ensuring a smoother and more robust treatment of discontinuities.

The main issue arises when the profiles present discontinuities. In such cases, if the $\overline{\psi }$ level at which $J_{\varphi }$ is discontinuous falls into a bend of elements, the integrals are calculated with a wrong value of $J_{\varphi }\left({\overline{\psi }}_g\right)$.

In order to properly handle the high steep current density regions within the plasma domain, as well as the possible discontinuities of $J_{\varphi }$ respect to $\overline{\psi }$, the current density is reformulated as follows:

$$L\psi =r{\displaystyle \frac{\nabla p\bullet \nabla \psi }{{\left|\nabla \psi \right|}^2}}+{\displaystyle \frac{1}{{\mu }_{0}r}}{\displaystyle \frac{\nabla \left(\frac{f^2}{2}\right)\bullet \nabla \psi }{{\left|\nabla \psi \right|}^2}}, \mathrm{with} p=p\left(\overline{\psi },\alpha \right), f=f\left(\overline{\psi },\alpha \right)$$

(21)

where $\nabla p$ and $\nabla \left({\displaystyle \frac{f^2}{2}}\right)$, are computed from $p$ and $f$ defined at the mesh nodes according to their respective profile functions. In this way, spatial discontinuities resulting from discontinuities in ${\displaystyle \frac{dp}{d\psi }}$ and ${\displaystyle \frac{d\left(\frac{f^2}{2}\right)}{d\psi }}$ are smoothed out. As a result, the computation of $J_{\varphi }$ (the flux source) becomes more accurate, and consequently, so does the evaluation of the flux itself. The singularity in the right-hand side of (21) wherever $\left|\nabla \psi \right|=0$ can easily be handled. Indeed, the flux gradient is typically zero for well defined values of $\overline{\psi }$, namely at the magnetic axis ($\overline{\psi }=0$) and at the null points at the plasma boundary ($\overline{\psi }=1$).

The nodal formulation is particularly suited in two cases that are beyond the scope of this paper, namely:

• the free boundary problem in cases where the current density is non-zero at the plasma edge;
• the derivation of the linearized plasma response to external currents and profile variations.

3.3. Accurate Calculation of the Magnetic Flux and Field

The above MHD equilibrium calculation with linear triangles provides a piecewise constant value of the flux gradient and yield magnetic field discontinuities across elements.

To tackle these issues, especially when the derivative of the magnetic field is needed, two different actions are undertaken. The first action is a mesh refinement wherever the discontinuities in the flux gradient and/or the gradients of the current density exceed prescribed thresholds.

The other method is a post-processing algorithm that requires very limited computational effort while ensuring the regularity of field derivatives, which is crucial for certain applications. This method is based on Helmholtz’s theorem [17] and employs the finite element method (FEM) with linear shape functions to calculate $\partial \psi /\partial r$ and $\partial \psi /\partial z$ in the second order calculation. With the technique proposed in [18], the poloidal magnetic field can be determined by post-processing the standard FEM flux map as the solution to the minimization problem:

$$min{\displaystyle \frac{1}{2}}\left\{{\int }_{\Omega }\left[{{\alpha }_D\left(\nabla \bullet {\mathit{B}}_{pol}\right)}^2+{{\alpha }_C\left(\nabla \times {\mathit{B}}_{pol}-{\mu }_0{J}_{\mathit{\varphi }}{i}_{\mathit{\varphi }}\right)}^2\right]r^{2}d\Omega +{\int }_{\partial \Omega }{\alpha }_B{\left(\mathit{n}\bullet {\mathit{B}}_{pol}-{\displaystyle \frac{\partial /\partial t}{r}}\right)}^{2}rd\Gamma \right\}$$

(22)

where, from (4), ${\mathit{B}}_{\mathit{p}\mathit{o}\mathit{l}}={\displaystyle \frac{\nabla \psi \times {\mathit{i}}_{\mathit{\varphi }}}{r}}$, while ${\alpha }_D, {\alpha }_C and {\alpha }_B$ are suitable non-dimensional constants.

The quantities $J_{\mathit{\varphi }}$, and $/ t$, which is the tangential component of $\nabla$, are obtained from the first order FEM solution of the G-S equation. This formulation ensures continuity and achieves a convergence rate of order O(h²) also for ${\mathit{B}}_{\mathit{p}\mathit{o}\mathit{l}}$, where h is the mesh size. These properties are highly useful, offering a high level of accuracy and enabling the derivation of reliable linearized models with limited computational effort, especially when compared to other methods, such as [19].

3.4. Optimization of Poloidal Field Coil Currents

The optimization of the number, position and dimension of the Poloidal Field (PF) coils in a tokamak is a challenging task, due to linear and nonlinear constraints related to maintenance requirements, and technological limits on electric power, current density, magnetic fields, and vertical forces [8,9,10,11,20]. Here, the problem of the optimization of the currents of the PF coils in air core tokamaks while keeping their positions fixed has been addressed [4,5,6,7]. It is based on the linearization of the GS equation around the desired plasma equilibrium. An iterative optimization problem, subject to both linear and quadratic constraints, is then solved.

The total flux per radian ${\psi }_b$ at the plasma boundary is the sum of two terms:

$${\psi }_b={\psi }_{bpl}+{\psi }_{bPF}$$

(23)

Here ${\psi }_{bPF}$, the flux pre radian due to the PF coils, can readily be computed using Green functions based on elliptic integrals [21]. The calculation of the plasma contribution ${\psi }_{bpl}$ using elliptic integrals may be cumbersome, since the boundary is in contact with the source currents. A solution, having the same degree of accuracy as FEM can be obtained by surrounding the fixed mesh boundary (i.e., the plasma contour) with a halo region ${\Omega }_h$ made of a few layers of elements, and solving:

$$L{\psi }_{pl}=J_{\varphi }\mathrm{in} {\Omega }_{pl} {\Omega }_h, \mathrm{with} {\psi }_{pl}={\psi }_{p{l}_G}\mathrm{on} \partial {(\Omega }_{pl}\mathrm{ }{\Omega }_h)$$

(24)

where $J_{\varphi }$ is given by the solution of Section 2.1 in ${\Omega }_{pl}$, and vanishes in ${\Omega }_h$.

The boundary condition ${\psi }_{p{l}_G}$ on $\partial {(\Omega }_{pl}{\Omega }_h)$ is computed using Green functions from $J_{\varphi }$. This method gets rid of the singularities, since the minimum distance between the points of the boundary and the sources is the thickness of the layer.

The PF coil currents are then optimized using one of the methods described in [8,9,10,11,20]. Since the optimized flux per radian is hardly equal to the difference between ${\psi }_b$ and ${\psi }_{bpl}$, it is necessary to check whether the free boundary equilibrium is sufficiently close to the desired fixed boundary configuration. This task is fulfilled by a specific version of the CREATE-NL+ code, solving the free boundary equilibrium code in a semicircle ${\Omega }_{CNL+}$ of the rz plane in the first and fourth quadrants, adding further elements to the mesh used for solving the fixed boundary problem. In this version of the free boundary code, the PF equilibrium currents are treated as a set of filaments at the same locations used for Gauss integration of elliptic integrals.

3.5. Linearization

The formulation described in this paper can effectively be used to derive a linearized plasma response model. The description of the details of the linearization procedure is outside the scope of the present paper. However, it is worth noticing that linearized models obtained with this procedure are extremely accurate and can efficiently be employed for plasma identification, magnetic control, and optimization of PF coil current waveforms.

4. Results and Discussion

This section shows the flexibility of parameterization (13) and demonstrates the reliability of the CREATE-FXB code with a number of benchmark cases, obtaining a very good agreement with available analytical solutions and other numerical results in the literature.

4.1. Envelopes for Current Density Parameterization

This section compares the limits of validity of three parameterization choices for the current density in the ${\beta }_{pol}$-$l_i$ plane for an elliptic plasma (Figure 1a).

The CREATE-FXB parameterization proposed here in (13) is able to obtain plasmas with a wide range of ${\beta }_{pol}$ and $l_i$ values. Figure 1 shows that the upper bounds (more than 20 for $l_i$ and more than 4 for ${\beta }_{pol}$) are clearly larger than those obtained with the other parameterizations. As stated in Section 2, the CREATE-FXB parameterization (13) is also able to obtain low values of $l_i$, i.e., down to about 0.2 like McCarthy’s family, whereas the corresponding lowest value for CREATE-L is about 0.5.

4.2. Circular Plasma with a Large Aspect Ratio and Discontinuous Current Density Profile

The robustness of the nodal formulation is shown by comparing the standard and new numerical solutions with the analytical one, in a case with a very large aspect ratio and a discontinuous current density profile. The test case is hence a circular plasma with a minor radius $a=1m$, a major radius $R_0\gg a$, and current density profile defined by:

$$I_{pl}=1 MA, {\displaystyle \frac{d\left(\frac{f^2}{2}\right)}{d\psi }}=0, {\displaystyle \frac{dp}{d\psi }}=\left\{\begin{array}{l}\begin{array}{c}{\displaystyle \frac{J_0}{R_0}} \overline{\psi }\le 0.4191 \end{array}\\ \begin{array}{c} \\ 0 \overline{\psi }>0.4191\end{array}\end{array} \right.$$

(25)

This 1D problem, which has a straightforward analytical solution, is solved as a 2D problem using FEM. Figure 2 shows the azimuthal field in the plasma region evaluated both using the standard and the nodal formulation of the current density $J_{\varphi }$ in the numerical solver. Moreover, for each formulation, the 2^nd order procedure has been applied.

Two aspects are highlighted in Table 1, which reports the values of the numerical error. Firstly, the nodal formulation provides better results when the discretization is coarse. Secondly, higher accuracy and smoothness of the magnetic field is obtained by post-processing the FEM solution with the method introduced in [18] and summarized in Section 3.3, no matter about the involved $J_{\varphi }$ formulation.

4.3. Non-Circular Low Aspect Ratio Plasma Equilibrium

The procedure has been applied also to a low-aspect ratio plasma equilibrium with an available analytical solution. The reference case is described in [22], with

$${\displaystyle \frac{dp}{d\psi }}={\displaystyle \frac{8{\psi }_0}{{\mu }_0{R}_0^4}}\left(1+{\alpha }^2\right), f{\displaystyle \frac{df}{d\psi }}=0$$

where ${\psi }_0$ is the poloidal flux per radian at the magnetic axis, $R_0$ is the radial coordinate of the magnetic axis, and $\alpha$ is a non-dimensional constant (there is a typo in [22] in the exponent of $R_0)$.

The analytic solution of the Grad-Shafranov equation is given by:

$$\psi \left(R,Z\right)={\displaystyle \frac{{\psi }_0{R}^2}{R_0^4}}\left(2R_0^2-R^2-4{\alpha }^2{Z}^2\right)$$

By fixing the previous parameters as follows:

$${\psi }_0=5{\displaystyle \frac{Vs}{rad}}, R_0=2m, \alpha =1$$

and choosing as boundary flux per radian:

$${\psi }_b=1{\displaystyle \frac{Vs}{rad}}$$

the $\left(R,Z\right)$ coordinates of the separatrix $\psi \left(R,Z\right)=$ ${\psi }_b$are shown in Figure 3a.

In this case the nodal and the standard formulations yield exactly the same results, since the current density is linear in each triangular element. Figure 3b shows the trend of the numerical error on the flux with the infinity norm as a function of the maximum element size $h_{max}$, which scales with $h_{max}^2$, as shown in [16].

4.4. Optimization of PF Equilibrium Currents in DTT

The fixed boundary solver has been used to determine the set of PF current needed to keep in equilibrium a single plasma in DTT.

Figure 4 shows the target plasma boundary at 27 s in the DTT plasma scenario with $I_{pl}$=5.5 MA, ${\psi }_b=0.70 Vs/rad$, ${\beta }_{pol}=0.10$, and $l_i=0.80$ [23]. The CREATE-FXB code has firstly calculated the fixed boundary equilibrium in ${\Omega }_{pl}$, and then the flux per radian due to the plasma ${\psi }_{pl}$, obtained in the region ${\Omega }_{pl} {\Omega }_h$, using the procedure described in Section 3.4.

Afterwards, the CREATE-FXB code minimizes ${‖\Delta {\psi }_{PF}‖}_{\infty }$, which is the difference between ${\psi }_{PF}$ (flux per radian due to the PF system at the plasma boundary) and the desired equilibrium value, namely ${\psi }_b$ - ${\psi }_{bpl}$. The PF system includes 12 circuits and 24 coils. The details are given in [23]. Figure 5a shows the ${‖\Delta {\psi }_{PF}‖}_{\infty }$ as a function of the maximum current in the circuits ${‖I_{PF}‖}_{\infty }$. Figure 5b shows the optimal currents in comparison with those reported in [23]. The slight improvement of CREATE-FXB is due to a different choice of the figure of merit, since the procedure used in [23] optimizes the entire scenario adding further constrains.

The validity and the effectiveness of the CREATE-FXB procedure is in any case demonstrated in Figure 5c, where the free boundary equilibrium obtained by the CREATE-NL+ code using the PF circuit currents obtained by CREATE-FXB provides a plasma boundary almost coincident with the desired shape, since the maximum distance is less than 9 mm.

5. Conclusions

In this work, CREATE-FXB is presented as a new fixed-boundary solver for the Grad–Shafranov equation, based on the Finite Element Method with an innovative technique. The code provides an efficient and accurate tool for plasma equilibrium analysis in tokamaks, offering robust handling of non-smooth current density profiles, second-order accuracy in the evaluation of magnetic flux and field, and an integrated framework for PF coil current optimization.

Validation against the free-boundary code CREATE-NL+ has confirmed the reliability of the approach for scenario design and equilibrium current estimation. Future extensions will focus on the systematic use of the derived linearized models for plasma identification and control, on the integration of plasma response-based scenario optimization strategies.

Overall, CREATE-FXB constitutes a flexible and accurate complement to the CREATE code family, enhancing its capability to address fixed-boundary problems in tokamak applications.