2.1. Current discretization

Typically, continuous current profiles are provided by analytical models or after equilibrium reconstruction using experimental data. We now discuss how these profiles can be represented in the framework of MRxMHD. Consider an externally driven current profile, e.g. ECCD, provided as the enclosed toroidal current as a function of the toroidal magnetic flux, i.e. $I_{\phi ,\textrm {ECCD}}(\psi _t)$, and a pressure-driven current profile, e.g. the bootstrap current, provided similarly as the enclosed toroidal current as a function of the toroidal flux, $I_{\phi ,\textrm {BS}}(\psi _t)$. We also assume that the pressure profile, $p(\psi _t)$, the number of volumes, $N_{\textrm {vol}}$, and their enclosed toroidal fluxes, $\{\psi _{t,l}\}_{l=1,\ldots ,N_{\textrm {vol}}}$, are given (see figure 2). The question of how many volumes and where to position their interfaces to best represent a given pressure profile is not addressed in this paper.

Figure 2. Sketch of a pressure profile as a function of the toroidal flux. Blue, continuous pressure profile obtained via experiment or analytical model; red, SPEC discretized pressure profile; black dashed lines, volume interfaces.

A proposed representation of these current density profiles in MRxMHD is achieved as follows. The ECCD current is an externally driven, parallel current and is thus represented as a volume current since it flows parallel to the field lines; on the other hand, the bootstrap current is a pressure-driven, self-generated current and is represented as a surface current, since it is localized at the pressure gradients. The volume currents are obtained by integrating the externally driven current density in each volume (figure 3), which is simply given by the difference

(2.8)\begin{equation} I^v_{l,\phi} = I_{\phi,\textrm{ECCD}}(\psi_{t,l}) - I_{\phi,\textrm{ECCD}}(\psi_{t,l-1}), \end{equation}

and the surface currents are obtained by integrating the pressure driven current density around each interface (figure 4), which is expressed as

(2.9)\begin{equation} I^s_{l,\phi} = I_{\phi,\textrm{BS}}(\psi_{l,\textrm{out}}) - I_{\phi,\textrm{BS}}(\psi_{l,\textrm{in}}), \end{equation}

with

(2.10)\begin{gather} \psi_{l,\textrm{in}} = \begin{cases} 0, & \text{if } l=1,\\ \dfrac{\psi_{t,l-1} + \psi_{t,l}}{2}, & \text{otherwise}, \end{cases} \end{gather}

(2.11)\begin{gather}\psi_{l,\textrm{out}} = \begin{cases} \psi_{a}, & \text{if } l=N_{\textrm{vol}}-1,\\ \dfrac{\psi_{t,l} + \psi_{t,l+1}}{2}, & \text{otherwise}, \end{cases} \end{gather}

with $\psi _{a}$ the total toroidal magnetic flux enclosed by the plasma. In (2.10)–(2.11), care has been taken for the first and last interfaces, where the surface of integration has been extended to include the current density from the magnetic axis and up to the plasma boundary. Note that this difference in the definition of the first and last surface currents vanishes as the number of volumes $N_{\textrm {vol}}$ is increased. Equations (2.8)–(2.11) are only one possible discretization of the continuous current profiles, proposed by the authors for illustration. Advantages of this particular representation are: (i) that the total toroidal current is always preserved; (ii) that the currents are approximately localized at the same location in the discretized than in the continuous case. The following sections do not depend on the particular choice of discretization of the currents.

Figure 3. Sketch of externally driven current density (red curve). Coloured area corresponds to the MRxMHD volume current. Black dashed lines represent volume interfaces.

Figure 4. Sketch of pressure driven current density. Coloured area corresponds to the MRxMHD surface current. Black dashed lines represent volume interfaces.

3.1. The SPEC code

The SPEC code can run in three different geometries, namely in slab (Loizu & Hudson Reference Loizu and Hudson2019), cylindrical (Loizu et al. Reference Loizu, Hudson, Helander, Lazerson and Bhattacharjee2016a) and toroidal geometry (Loizu et al. Reference Loizu, Hudson and Nührenberg2016b). The coordinates used to describe position are

(3.1)\begin{equation} \boldsymbol{x} = \begin{cases} \theta\hat{\boldsymbol{i}} + \phi\hat{\boldsymbol{j}} + R\hat{\boldsymbol{k}}, & \text{in slab geometry},\\ R\cos\theta\hat{\boldsymbol{i}} + R\sin\theta\hat{\boldsymbol{j}} + \phi\hat{\boldsymbol{k}}, & \text{in cylindrical geometry},\\ R\cos\phi\hat{\boldsymbol{i}} + R\sin\phi\hat{\boldsymbol{j}} + Z\hat{\boldsymbol{k}}, & \text{in toroidal geometry}, \end{cases} \end{equation}

where $\{\hat {\boldsymbol {i}},\hat {\boldsymbol {j}},\hat {\boldsymbol {k}}\}$ is the unitary Cartesian basis, $\theta$ is a poloidal angle and $\phi$ is the usual toroidal angle. At each interface $\mathcal {I}_l$, the geometry of the surface is described by decomposing the coordinates $R$ and $Z$ on a Fourier basis,

(3.2)\begin{gather} R_l(\theta,{\phi}) = \sum_{m=0}^{M_{\textrm{pol}}}\sum_{n={-}N_{\textrm{tor}}}^{N_{\textrm{tor}}}R_{l,m,n} \cos(m\theta-nN_p\phi), \end{gather}

(3.3)\begin{gather}Z_l(\theta,{\phi}) = \sum_{m=0}^{M_{\textrm{pol}}}\sum_{n={-}N_{\textrm{tor}}}^{N_{\textrm{tor}}}Z_{l,m,n} \sin(m\theta-nN_p\phi), \end{gather}

where $N_p$ is the number of field periods, $M_{\textrm {pol}}$ and $N_{\textrm {tor}}$ are the poloidal and toroidal mode numbers above which Fourier series are truncated, i.e. $m=\{0,\ldots ,M_{\textrm {pol}}\}$, $n=\{-N_{\textrm {tor}},\ldots ,N_{\textrm {tor}}\}$, and stellarator symmetry has been assumed for simplicity. Between interfaces $l$ and $l+1$, coordinates are constructed by linear interpolation of $(R,Z)_l$ and $(R,Z)_{l+1}$ using a radial-like coordinate $s$. Geometrical degrees of freedom are $\{R_{l,m,n}, Z_{l,m,n}\}\equiv \{x_i\}_{i=1,\ldots ,N}$, with, e.g. in the case of a fixed-boundary, stellarator symmetric equilibrium in toroidal geometry,

(3.4)\begin{equation} N = (N_{\textrm{vol}}-1)\left\{2\left[1+N_{\textrm{tor}} + M_{\textrm{pol}} \left( 2N_{\textrm{tor}} + 1 \right) \right] -1\right\}. \end{equation}

The magnetic field is represented by the magnetic vector potential, $\boldsymbol {B}_l=\boldsymbol {\nabla }\times \boldsymbol {A}_l$, which is described by using a Chebyshev–Fourier basis,

(3.5)\begin{equation} A_{l,i} = \sum_{m,n}\sum_{k=0}^{L_{\textrm{rad}}} A_{l,i,k,m,n} T_k(s)\cos(m\theta-nN_p\phi), \end{equation}

where $L_{\textrm {rad}}$ is the radial resolution, $T_k$ is the Chebyshev polynomial of order $k$ and $i=\{\theta ,\phi \}$. Note that gauge freedom is used to set $A_{l,s}=0\ \forall l$. The coefficients $A_{l,i,k,m,n}$ are the vector potential degrees of freedom, packed in a single array $\boldsymbol {a}_l$ per volume.

Fixed-boundary SPEC can solve the MRxMHD system defined by (2.4) and (2.5) for a given plasma boundary geometry, pressure profile $\{p_l\}_{l=1,\ldots ,N_{\textrm {vol}}}$ and profiles $\{\psi _{t,l}, \psi _{p,l}, \mu _l\}_{l=1,\ldots ,N_{\textrm {vol}}}$. In this case, the Beltrami equation in volume $l$, (2.4), can be cast into a linear system

(3.6)\begin{equation} \boldsymbol{G}_l\boldsymbol{\cdot}\boldsymbol{a}_l = \boldsymbol{C}_l, \end{equation}

where $\boldsymbol {G}_l$ and $\boldsymbol {C}_l$ are matrices that depend only on the geometrical degrees of freedom $\{x_i\}_{i=1,\ldots ,N}$, and linearly on $\psi _{t,l}$, $\psi _{p,l}$ and $\mu _l$.

The constrained quantities $\{\mu _l,\psi _{t,l},\psi _{p,l}\}_{l=1,\ldots ,N_{\textrm {vol}}}$ can be replaced by other independent triplets such as $\{\psi _{t,l},\ {\raise.3pt-\kern-5pt\iota} _l^+,\ {\raise.3pt-\kern-5pt\iota} _l^-\}_{l=1,\ldots ,N_{\textrm {vol}}}$, with ${\raise.3pt-\kern-5pt\iota} _l^\pm$ the rotational transform evaluated on the inner and outer side of interface $\mathcal {I}_l$, respectively. In this case, SPEC iterates on $\{\psi _{p,l}, \mu _l\}_{l=1,\ldots ,N_{\textrm {vol}}}$ until the solution satisfies the input rotational transform profile. In what follows, this method will be referred to as the ‘rotational transform constraint'. Finally, SPEC uses a Newton method to iterate on $\{x_i\}_{i=1,\ldots ,N}$ until the force balance, (2.5), is satisfied. This is achieved by expanding the force imbalance on a Fourier basis,

(3.7)\begin{equation} \left[\left[p+\frac{B^2}{2\mu_0}\right]\right]_l = \sum_{m=0}^{M_{\textrm{pol}}}\sum_{n={-}N_{\textrm{tor}}}^{N_{\textrm{tor}}} F_{l,m,n} \cos(m\theta-nN_p\phi), \end{equation}

and packing the Fourier components $F_{l,m,n}$ in a single array $\boldsymbol {F}$ of size $N$. Force balance is satisfied once all the Fourier modes $F_{l,m,n}$, herein $\{F_j\}_{j=1,\ldots ,N}$, are sufficiently small. Analytical derivatives of $\{F_j\}_{j=1,\ldots ,N}$ with respect to $\{x_i\}_{i=1,\ldots ,N}$ accelerate substantially the code.

It is worth noting that even though the magnetic helicity is not conserved when SPEC iterates on $\{x_i\}_{i=1,\ldots ,N}$ to find an equilibrium that matches a given input $\{\psi _{t,l}, \psi _{p,l}, \mu _l\}_{l=1,\ldots ,N_{\textrm {vol}}}$ or $\{\psi _{t,l},\ {\raise.3pt-\kern-5pt\iota} _l^+,\ {\raise.3pt-\kern-5pt\iota} _l^-\}_{l=1,\ldots ,N_{\textrm {vol}}}$, the final equilibrium satisfies the MRxMHD equilibrium equations ((2.4)–(2.5)). There is a magnetic helicity profile $\{\textrm {K}_l\}_{l=1,\ldots ,N_{\textrm {vol}}}$ corresponding to this equilibrium which is unknown a priori. Thus, the same equilibrium could have been found by minimizing the MRxMHD energy functional while keeping the magnetic helicity profile constant if the initial state had the same magnetic helicity profile (bifurcations are not considered in this paper). This capability is also available in SPEC, and details can be found in the literature (Hudson et al. Reference Hudson, Dewar, Dennis, Hole, McGann, Von Nessi and Lazerson2012a; Dennis et al. Reference Dennis, Hudson, Terranova, Franz, Dewar and Hole2013b).

Lastly, SPEC has been recently extended to allow free-boundary calculations (Hudson et al. Reference Hudson, Loizu, Zhu, Qu, Nührenberg, Lazerson, Smiet and Hole2020), where the plasma is surrounded by a vacuum region bounded by a computational boundary. The solution in the vacuum region, which is a special instance of a Beltrami field, with $\mu =0$, is determined by a different set of scalars, namely the total poloidal current passing through the torus hole, $I_{\textrm {coil}}$, and the total toroidal plasma current, $I_{\textrm {pl}}$. Using Ampere's law, these currents can be written as

(3.8)\begin{gather} \mu_0 I_{\textrm{coil}} = 2{\rm \pi}\tilde{B}^+_{V,\phi}, \end{gather}

(3.9)\begin{gather}\mu_0 I_{\textrm{pl}} = 2{\rm \pi}\tilde{B}^-_{V,\theta}, \end{gather}

with $\tilde {B}^+_{V,\phi }$ the $m=n=0$ Fourier mode of the covariant toroidal magnetic field on the computational boundary and $\tilde {B}^-_{V,\theta }$ the $m=n=0$ Fourier mode of the covariant poloidal magnetic field on the outer side of the plasma boundary. Additionally, the computational boundary is not necessarily a flux surface, thus $\boldsymbol {B}_V\boldsymbol {\cdot }\hat {\boldsymbol {n}}_c \neq 0$ with $\boldsymbol {B}_V$ the magnetic field in the vacuum region and $\hat {\boldsymbol {n}}_c$ a unit vector normal to the computational boundary. The solution satisfying the constraints $I_{\textrm {coil}}$ and $I_{\textrm {pl}}$ is obtained by iteration over the poloidal and toroidal fluxes $\{\psi _{p,V},\psi _{t,V}\}$ in the vacuum volume.

3.2. Constraining the toroidal current profiles in SPEC

The SPEC code has been extended to allow the triplet $\{\psi _{t,l}, I^v_{l,\phi }, I^s_{l,\phi }\}_{l=1,\ldots ,N_{\textrm {vol}}}$ as a constraint, herein ‘current constraint’. In the case of the rotational transform constraint, SPEC finds the solution to the linear system (2.4) volume by volume and iterates on $\{\mu _l, \psi _{p,l}\}$ until the field has the desired rotational transform at the volume interfaces. In the case of the current constraint, the constants $\{\mu _l\}_{l=1,\ldots ,N_{\textrm {vol}}}$ are determined using (2.6), without the need for iterations, and this directly constrains the value of volume currents $\{I^v_{l,\phi }\}_{l=1,\ldots ,N_{\textrm {vol}}}$.

Regarding the poloidal fluxes, it can be shown (Appendix A) that the surface currents depend linearly on the poloidal fluxes, i.e.

(3.10)\begin{equation} \boldsymbol{I}^s=\boldsymbol{M}\boldsymbol{\psi}_p+\boldsymbol{Q}, \end{equation}

where $\boldsymbol {I}^s$ and $\boldsymbol {\psi }_p$ are arrays containing all $\{I^s_l\}_{l=1,\ldots ,N_{\textrm {vol}}-1}$ and all $\{\psi _{p,l}\}_{l=2,\ldots ,N_{\textrm {vol}}}$, respectively, and the matrix $\boldsymbol {M}$ and the array $\boldsymbol {Q}$ depend only on the geometry of the interfaces $\{x_i\}_{i=1,\ldots ,N}$. In this section, we consider the geometry, toroidal fluxes and the constants $\{\mu _l\}_{l=1,\ldots ,N_{\textrm {vol}}}$ to be fixed and seek how the poloidal flux profile $\boldsymbol {\psi }_p$ has to be constrained in order to obtain a surface current profile matching the input profile $\boldsymbol {I}^s$.

The unknown $\boldsymbol {Q}$ is eliminated by subtracting (3.10) evaluated at two different values of $\boldsymbol {\psi }_p$, i.e. evaluated once at $\overline {\boldsymbol {\psi _p}}$ and once at $\boldsymbol {\psi }_p$,

(3.11)\begin{equation} \boldsymbol{M} (\overline{\boldsymbol{\psi}_p} - \boldsymbol{\psi}_p) = \overline{\boldsymbol{I}^{s}} - \boldsymbol{I}^{s}, \end{equation}

where $\overline {\boldsymbol {\psi }_p}$ is an arbitrary choice of poloidal fluxes, and $\overline {\boldsymbol {I}^s}$ is the surface current profile calculated from the Beltrami fields $\{\overline {\boldsymbol {a}_l}\}_{l=1,\ldots ,N_{\textrm {vol}}}$ obtained when the poloidal fluxes are constrained to the values $\overline {\boldsymbol {\psi }_p}$.

The matrix $\boldsymbol {M}$ is evaluated by taking the derivatives of (3.10) with respect to the poloidal fluxes, i.e.

(3.12)\begin{equation} M_{ij} = \frac{\partial I^s_{i,\phi}}{\partial \psi_{p,j}}, \end{equation}

which leads to the following bidiagonal matrix:

(3.13)\begin{equation} \boldsymbol{M} = \frac{2{\rm \pi}}{\mu_0} \begin{bmatrix}\dfrac{\partial \tilde{B}^-_{\theta,2}}{\partial{\psi_{p,2}}} & 0 & \cdots & \cdots & 0 \\ -\dfrac{\partial \tilde{B}^+_{\theta,2}}{\partial{\psi_{p,2}}} & \dfrac{\partial \tilde{B}^-_{\theta,3}}{\partial{\psi_{p,3}}} & 0 & \cdots & 0\\ 0 & -\dfrac{\partial \tilde{B}^+_{\theta,3}}{\partial{\psi_{p,3}}} & \dfrac{\partial \tilde{B}^-_{\theta,4}}{\partial{\psi_{p,4}}} & 0 & 0 \\ \vdots & \ddots & \ddots & \ddots & 0\\ 0 \cdots & \cdots & 0 & -\dfrac{\partial \tilde{B}^+_{\theta,N_{\textrm{vol}}-1}}{\partial{\psi_{p,N_{\textrm{vol}}-1}}} & \dfrac{\partial \tilde{B}^-_{\theta,N_{\textrm{vol}}}}{\partial{\psi_{p,N_{\textrm{vol}}}}}\end{bmatrix}.\end{equation}

Derivatives of $\tilde {B}^\pm _{\theta ,l}$ with respect to the poloidal flux can be easily obtained by applying matrix perturbation theory on the linear system (3.6) (Hudson et al. Reference Hudson, Dewar, Dennis, Hole, McGann, Von Nessi and Lazerson2012a),

(3.14)\begin{equation} \boldsymbol{G}_l \boldsymbol{\cdot} \left[\frac{\partial}{\partial\psi_{p,l}}\boldsymbol{a}_l\right] ={-} \left[\frac{\partial}{\partial\psi_{p,l}}\boldsymbol{G}_l\right] \boldsymbol{\cdot} \boldsymbol{a}_l + \frac{\partial}{\partial\psi_{p,l}}\boldsymbol{C}_l.\end{equation}

Due to the linear nature of (3.10), the coefficients of the matrix $\boldsymbol {M}$, i.e. (3.12), are independent of $\boldsymbol {\psi }_p$ and thus can be evaluated once at any arbitrary value of $\boldsymbol {\psi }_p$. We thus conveniently evaluate them at $\overline {\boldsymbol {\psi }_p}$. Equation (3.11) is then solved to obtain the poloidal flux profile $\boldsymbol {\psi }_p$.

Finally, instead of solving a second time the Beltrami equation, (2.4), at $\boldsymbol {\psi }_p$, we take advantage of the linear dependence of $\boldsymbol {a}$ on $\boldsymbol {\psi }_p$, (i.e. 3.6), and solve

(3.15)\begin{equation} A_{l,i} = \overline{A_{l,i}} - \frac{\partial {A_{l,i}}}{\partial {\psi_{p,l}}} (\overline{\psi_{p,l}} - \psi_{p,l}), \end{equation}

where $\overline {A_{l,i}}$ is one element of $\overline {\boldsymbol {a}_l}$. The algorithm flow is summarized in figure 5.

Figure 5. Flow of the algorithm used to constrain the net toroidal current profiles for a given toroidal flux profile and geometry.

In the case of a free-boundary computation, the toroidal flux in the vacuum region is varied to satisfy the poloidal linking current $I_{\textrm {coil}}$. This slightly modifies the linear system (3.11). Details can be found in Appendix B.

Constraining the toroidal current profile takes away the control of other profiles such as the rotational transform or the magnetic helicity. However, as for the case of the rotational transform constraint, the equilibrium can be accessed by a relaxation process at constant magnetic helicity if the final magnetic helicity is known a priori. The MRxMHD equations are thus still satisfied by an equilibrium obtained by constraining the toroidal current profiles.

The new current constraint has been parallelized with MPI in a similar fashion to the other constraints. Each volume is associated with one central processing unit (CPU); since the solution to the Beltrami equation (2.4) in a volume is independent from other volumes, each CPU can solve the linear system (3.6) in parallel. Finally, the master CPU gathers all required derivatives to construct the matrix $\boldsymbol {M}$ and solves the linear system (3.11), before broadcasting the values of $\{\psi _{p,l}\}_{l=2,\ldots ,N_{\textrm {vol}}}$ and $\{\boldsymbol {a}_l\}_{l=1,\ldots ,N_{\textrm {vol}}}$ to all CPUs.

3.3. Force gradient

The Newton algorithm used in SPEC to iterate on the interfaces' geometry uses analytic derivatives, which is faster than finite differentiation. To keep good performance while using the current constraint, derivatives of the force Fourier coefficients, $\{F_j\}_{j=1,\ldots ,N}$, with respect to the geometrical degrees of freedom, $\{x_i\}_{i=\{1,\ldots ,N\}}$, at constant $\{\psi _{t,l}, I^v_{l,\phi },\ I^s_{l,\phi }\}_{l=1,\ldots ,N_{\textrm {vol}}}$ are provided. Derivatives are first evaluated in real space and then Fourier transformed. Using the chain rule,

(3.16)\begin{gather} \frac{\textrm{d}}{\textrm{d} x_i}\left[\left[p + \frac{B^2}{2\mu_0}\right]\right]_l = \frac{1}{\mu_0}\left( B^-_{l+1}\frac{\textrm{d} B^-_{l+1}}{\textrm{d} x_i} - B^+_l \frac{\textrm{d} B^+_{l}}{\textrm{d} x_i}\right), \end{gather}

(3.17)\begin{gather}\frac{\textrm{d} B^{{\pm}}_l}{\textrm{d} x_i} = \frac{\partial B^\pm_l}{\partial x_i} + \frac{\partial B^\pm_l}{\partial \psi_{p,l}}\frac{\textrm{d} \psi_{p,l}}{\textrm{d} x_i} + \frac{\partial B^\pm_l}{\partial \mu_l}\frac{\textrm{d}\mu_l}{\textrm{d} x_i}+ \frac{\partial B^\pm_l}{\partial \psi_{t,l}}\frac{\textrm{d}\psi_{t,l}}{\textrm{d} x_i}, \end{gather}

where $B^-_l$, $B^+_l$ are the magnetic field strength on the inner and outer side of volume $l$, respectively, and the pressure, $p_l$, is considered constant in each volume with respect to variations in the geometry, $\mu _l$ and $\psi _{p,l}$. Note that all derivatives are taken at constant toroidal flux, volume current and surface current. Enforcing $\textrm {d}\psi _{t,l} / \textrm {d}x_i=0$ and $\textrm {d} I^v_{l,\phi } / \textrm {d}x_i=0$ leads to $\textrm {d} \mu _l/\textrm {d}x_i=0$ using (2.6). The surface current constraint, $\textrm {d} I^s_{l,\phi }/\textrm {d}x_i=0$, leads to a system of coupled equations using (2.7),

(3.18)\begin{equation} \frac{\partial \tilde{B}^-_{l+1,\theta}}{\partial x_i} + \frac{\partial \tilde{B}^-_{l+1,\theta}}{\partial \psi_{p,l+1}} \frac{\partial\psi_{p,l+1}}{\partial x_i} - \frac{\partial \tilde{B}^+_{l,\theta}}{\partial x_i} - \frac{\partial \tilde{B}^+_{l,\theta}}{\partial \psi_{p,l}}\frac{\partial\psi_{p,l}}{\partial x_i} = 0, \end{equation}

which can be written as a linear system using the matrix $\boldsymbol {M}$ defined in (3.13),

(3.19)\begin{equation} \boldsymbol{M} \boldsymbol{\cdot} \begin{bmatrix} \dfrac{\textrm{d}\psi_{p,2}}{\textrm{d}x_i}\\ \vdots\\ \dfrac{\textrm{d}\psi_{p,N_{\textrm{vol}}}}{\textrm{d}x_i} \end{bmatrix} = \frac{2{\rm \pi}}{\mu_0}\begin{bmatrix} \dfrac{\partial \tilde{B}^+_{\theta,1}}{\partial x_i} - \dfrac{\partial \tilde{B}^-_{\theta,2}}{\partial x_i} \\ \vdots \\ \dfrac{\partial \tilde{B}^+_{\theta,N_{\textrm{vol}}-1} } {\partial x_i} - \dfrac{\partial \tilde{B}^-_{\theta,N_{\textrm{vol}}}}{\partial x_i} \end{bmatrix}.\end{equation}

Derivatives of $\tilde {B}_{\theta ,l}$ with respect to $\psi _{p,l}$ and $x_i$ can be obtained by applying matrix perturbation theory to the Beltrami system (3.6). The solution of (3.19), together with (3.17) provides the required derivatives of the force with respect to the geometry. The derivatives of the Fourier components of the force, $\textrm {d}F_j/\textrm {d}x_i$, are obtained by taking the Fourier transform of (3.17) and are packed in a matrix $\boldsymbol {\nabla } F$ of size $N^2$, henceforth named ‘force gradient’. Appendix B provides details on the free-boundary case.

4.1. Verification in cylindrical geometry

We consider a fixed-boundary screw pinch MRxMHD equilibrium that only depends on the radius $R$ and whose solutions can be written analytically (Appendix C). We choose a set of somewhat arbitrary input parameters, i.e. a cylinder of minor radius $a=1$ and length $L=2{\rm \pi}$, $N_{\textrm {vol}}=3$, $p_l=0$ $\forall l\in \{1,2,3\}$, $\boldsymbol {\psi }_t = \{1/9, 4/9, 1\}\,\textrm {Tm}{}^2$, $\mu _0\boldsymbol {I}^v=\{0.2,0.2,0.4\}\,\textrm {Tm}$ and $\mu _0\boldsymbol {I}^s=\{-0.4,0.5\}\,\textrm {Tm}$, which uniquely define the analytical solution. The SPEC code is then run with the same input parameters and the solutions are compared (figure 6). Very good agreement between the analytical solution and the SPEC solution is obtained. Note that by constraining the toroidal current profiles, we lose control on the rotational transform profile, since only two profiles can be constrained in addition to the pressure profile. Hence discontinuities in the magnetic field components arise at the volume interfaces, even when there is no pressure, unless the input parameters are carefully selected.

Figure 6. Magnetic field components as a function of the radius in the case of a screw pinch. Solid and dashed lines, analytical solution as given in Appendix C; circles and triangles, SPEC solution using the current constraint.

The force gradient can also be expressed in terms of Bessel function integrals (see Appendix C). Figure 7 shows the normalized maximum absolute error between the force gradient obtained with SPEC and that obtained analytically as a function of the radial resolution $L_{\textrm {rad}}$. As $L_{\textrm {rad}}$ is increased, exponential convergence is observed up until $10^{-13}$, where the error in the evaluation of the Bessel integrals starts to dominate.

Figure 7. Semilogarithmic plot of the maximal normalized absolute error between the analytical force gradient, $\boldsymbol {\nabla } F_{\textrm {AN}}$, and the force gradient obtained from SPEC, $\boldsymbol {\nabla } F$, as a function of the radial resolution for the screw pinch case.

4.2. Verification in toroidal geometry

A verification is proposed here in the more complex case of a free-boundary, rotating ellipse (also called classical stellarator) equilibrium with five field periods ($N_p=5$), multiple poloidal and toroidal modes ($M_{\textrm {pol}}=4$ and $N_{\textrm {tor}}=2$) and seven plasma volumes ($N_{\textrm {vol}}=7$). The pressure is set to zero and the computational boundary is defined by

(4.1)\begin{gather} R = R_{00} + R_{10} \cos\theta + R_{11}\cos(\theta - N_p\phi), \end{gather}

(4.2)\begin{gather}Z = Z_{00} + Z_{10} \sin\theta + Z_{11}\sin(\theta - N_p\phi), \end{gather}

with $R_{00}=10\ \textrm {m}$, $Z_{00}=0\ \textrm {m}$, $R_{10}=-Z_{10}=1\ \textrm {m}$, $R_{11}=Z_{11}=0.25\ \textrm {m}$ and $N_p = 5$ (see figure 8). We suppose here that some hypothetical coils with a total current of $\mu _0I_{\textrm {coil}}=42.87$ Tm are able to generate a vacuum field without normal component to the computational boundary. The total toroidal magnetic flux in the plasma is set to $\psi _a=0.61\,\text {Tm}^2$ and the toroidal magnetic flux in the vacuum region is set to $\psi _{t,V}=1.39\,\text {Tm}^2$, adding up to a total toroidal magnetic flux enclosed by the computational boundary of $\psi _{t,V} + \psi _a = 2\,\text {Tm}^2$.

Figure 8. Three-dimensional plot of the classical stellarator boundary as described by (4.1)–(4.2). Colours indicate the magnetic field strength.

To the authors’ knowledge, no analytical solution to (2.4) and (2.5) exists in this geometry. The verification is thus carried out as follows. First, a rotational transform constraint case is run with an input ${\raise.3pt-\kern-5pt\iota}$-profile that is chosen to be $10\,\%$ larger than the vacuum rotational transform ${\raise.3pt-\kern-5pt\iota} _{\textrm {vac}}$, i.e. ${\raise.3pt-\kern-5pt\iota} = 1.10\times {\raise.3pt-\kern-5pt\iota} _{\textrm {vac}}$, so that there is a non-zero contribution from the current to the rotational transform. The volume and surface currents are evaluated from the obtained equilibrium and used to run a current constraint calculation to obtain a second equilibrium. The same initial guess for the geometry and the interfaces is used for both calculations. The rotational transform profile $\bar {{\raise.3pt-\kern-5pt\iota} }$ is then extracted from the second equilibrium and compared with the reference ${\raise.3pt-\kern-5pt\iota}$-profile.

The vacuum rotational transform profile, as well as the profiles ${\raise.3pt-\kern-5pt\iota}$ and $\bar {{\raise.3pt-\kern-5pt\iota} }$ are shown in figure 9(a). The toroidal current enclosed by the plasma is mostly contained in the volumes and adds up to a total of ${\sim }2.7\ \textrm {kA}$, see figure 9(b). As expected the surfaces currents $I^s_{\phi ,l}$ remain small ($<10^{-2}\ \textrm {kA}$), since there are no pressure gradients to drive them. The constraint on the rotational transform ${\raise.3pt-\kern-5pt\iota}$ is enforced on each side of the volumes’ interfaces, indicated by grey dashed lines on figure 9(a). The value of ${\raise.3pt-\kern-5pt\iota}$ at the computational boundary is not constrained. Agreement between ${\raise.3pt-\kern-5pt\iota}$ and $\bar {{\raise.3pt-\kern-5pt\iota} }$ up to a relative error of $\max (|{\raise.3pt-\kern-5pt\iota} -\bar {{\raise.3pt-\kern-5pt\iota} }| / |{\raise.3pt-\kern-5pt\iota} |)\sim 10^{-5}$ is observed, showing that the same equilibrium can be obtained using either constraint. The maximum error between both profiles decreases as the numerical resolution is increased (data not shown).

Figure 9. (a) Rotational transform profile versus effective minor radius. Red triangles, ${\raise.3pt-\kern-5pt\iota}$, input profile used in SPEC when run at fixed rotational transform; black dashed line, $\bar {{\raise.3pt-\kern-5pt\iota} }$, the output profile obtained from SPEC when run at fixed toroidal current profile; blue line, $\iota$-profile in vacuum; grey dashed lines, position of volume interfaces. (b) Total toroidal current enclosed by each volume. Surface currents (not plotted), $I^s_{\phi ,l}$, are smaller than $10^{-2}\ (\textrm {kA})$ and are negligible in comparison with the volume current.

To verify the force gradient components, we use a fourth-order centred finite difference formula (Fornberg Reference Fornberg1988),

(4.3)\begin{equation} \frac{\textrm{d} f}{\textrm{d}x} = \frac{f(x-2\Delta x) - 8f(x-\Delta x) + 8f(x+\Delta x) - f(x+2\Delta x)}{12\Delta x} +{ O(\Delta x^4)}, \end{equation}

to obtain $\boldsymbol {\nabla } F_{FD}$, i.e. a finite-difference estimate of the force gradient, and compare it with $\boldsymbol {\nabla } F$, i.e. the force gradient calculated in SPEC by using analytical derivatives. The finite difference estimate is evaluated by perturbing each geometrical degree of freedom $\{x_i\}_{l=1,\ldots ,N}$ by a constant value $\Delta R$. Convergence as $\Delta R\rightarrow 0$ is shown in figure 10. A convergence of order $O(\Delta R^4)$ is observed down to $\sim 10^{-11}$ for $\Delta R\sim 10^{-4}$. For lower values of $\Delta R$, the finite difference approximation error is dominated by round-off error. This shows that the analytical derivatives (the force gradient) are correctly implemented in SPEC.

Figure 10. Normalized maximum absolute error between SPEC force gradient and a finite difference estimate in the case of a rotating ellipse. The dashed line has slope of four.