JOREK3D: An extension of the JOREK nonlinear MHD code to stellarators

The stellarator, having been proposed by Lyman Spitzer in 1951, is one of the oldest plasma confinement concepts potentially applicable as a fusion power plant. However, early stellarators were plagued with problems stemming from neoclassical transport losses, leading to them being largely phased out in favor of tokamaks by the 1970s.^2–4 However, improved mathematical models and increased computational power, which became available by the late 1980s, allowed us to overcome the main challenges faced by the stellarator concept. Moreover, the revival of stellarators brought with it a new strategy for fusion research, where numerical modeling drives the development of future machines, as opposed to the traditional strategy, where smaller-scale machines had to be built and experimented on before advancing to larger-scale machines. The creation of Wendelstein 7-X is one example of the successful application of this new strategy. The advantages are clear: not only is it more cost-effective, but it also allows one to consider a much wider range of potential machine designs in a much shorter amount of time.² 

Most of the computational developments mentioned above focused on the optimization of stellarator equilibria. While there has been work on nonlinear magnetohydrodynamic (MHD) simulations of stellarators since the 1970s,⁵ this area is not as developed as stellarator optimization. Most early studies applied the straight stellarator approximation by neglecting toroidicity,^5–7 and it was not until the 2000s that fully 3D geometries were simulated.^8,9 At present, several well-known MHD codes exist, with a few of them, including M3D-C¹,¹⁰ M3D,⁸ and MIPS,¹¹ having been extended to stellarators. All three of these codes use full MHD on flux-surface-aligned grids, except for MIPS, which uses a cylindrical grid. NIMROD, another major tokamak code, is still in the process of being extended to stellarators,¹² although the tokamak capabilities of NIMROD are already enough to simulate a stellarator with an axisymmetric vacuum vessel.¹³ Also, the FLUXO nonlinear MHD code¹⁴ is applicable to stellarator geometries. However, the stellarator capabilities of these codes have not been used much so far. This study reports on a similar extension of JOREK, one of the leading nonlinear MHD codes for tokamaks,^1,15,16 to model stellarators. The work consists of two parts: first, a reduced MHD model compatible with three-dimensional geometries was derived by generalizing the ideas of Breslau et al., Izzo et al., and Strauss;^17–19 then, this model is implemented in the JOREK code and tested on a simple stellarator. The first part of the work has already been published in previous studies,^20,21 and so this study will present the results of the second part of this effort.

As discussed in our previous studies,^20,21 the magnetic field ansatz and equations in stellarator-capable reduced MHD involve a magnetic scalar potential χ, which represents the part of the magnetic field that is generated by the coils. In the tokamak limit, this potential reduces to $χ=F0ϕ$⁠, where $ϕ$ is the toroidal angle; however, a stellarator-capable code using this model will need to allow arbitrary scalar potentials. Fortunately, it is possible to analytically represent an arbitrary χ: since $∇χ$ is a magnetic field, it must be divergence-free, so $Δχ=0$⁠. One then needs to find a general solution to the Laplace equation in the toroidal coordinate system $(R,z,ϕ)$⁠. This was done by Dommaschk,²² who provides his solution as a sum over harmonics, where any particular solution is determined by the coefficients of these harmonics. Naturally, each harmonic individually satisfies the Laplace equation. In order to determine the coefficients for a particular equilibrium, one needs to first calculate the vacuum field on an $(R,z,ϕ)$ grid, which we do using the EXTENDER_P code.²³ 

Using the Dommaschk potential formulation for χ in conjunction with nonaxisymmetric flux-surface-aligned grids allows one to relatively efficiently simulate stellarators. The steps to run a stellarator simulation can then be summarized as follows:

1. Calculate an equilibrium for the stellarator in question using the GVEC code.²⁴ 

2. Use the output of GVEC to calculate the contribution to the stellarator's magnetic field from the coils (i.e., the curl-free/vacuum field) with the EXTENDER_P code.

3. Calculate the coefficients for the Dommaschk representation of the scalar potential from the output of EXTENDER_P.

4. Build a flux-surface-aligned grid from the geometry data in the GVEC solution and import it into JOREK.

5. Calculate the initial values for the reduced MHD variables from the GVEC solution.

6. Evolve the system implicitly in time using the stellarator reduced MHD equations in JOREK.

The GVEC²⁴ code mentioned above is a new fixed-boundary 3D MHD equilibrium solver, which follows the ideas of the well-established VMEC code,^25,26 assuming nested flux surfaces and using a constraint minimization of the MHD energy. In contrast to VMEC, the radial discretization is based on nonuniform B-splines of arbitrary order, allowing smooth representation of equilibrium quantities.

The rest of this paper is organized as follows. Section II states the reduced MHD equations that will be used throughout the rest of the paper. The same section also discusses the compatibility of full MHD equilibria with reduced MHD and shows that the error introduced by reduced MHD is small. Section III explains how the coefficients of the Dommaschk representation can be calculated from EXTENDER_P output, while Sec. IV explains how the initial conditions for the reduced MHD variables can be calculated from the GVEC equilibrium. In Sec. V, several simulations of stable equilibria in the Wendelstein 7-A stellarator²⁷ are presented to show that spurious instabilities and problems with maintaining equilibrium do not appear. Finally, Sec. VI shows tearing mode simulations in Wendelstein 7-A and benchmarks the growth rates against the CASTOR3D linear MHD code.^28,29 Section VII presents a similar benchmark for ballooning modes in the same device. In addition,  Appendix A briefly discusses the new nonaxisymmetric grid feature in JOREK, which allows the stellarator simulations to have flux surface-aligned grids. In  Appendix B, we show that the linearized ideal version of the model presented in Sec. II has a self-adjoint operator, and thus, ideal perturbations will have real eigenvalues.

Throughout this paper, we will use the reduced MHD model that was derived in Ref. 21. The advantage of reduced MHD is that it allows one to use a larger time step than full MHD by eliminating the shortest timescale in the system; even if implicit time stepping is used, accuracy will deteriorate if the time step is too much larger than the shortest timescale. In the tokamak limit, reduced MHD is well tested and can accurately model tearing and ballooning modes in a wide range of betas and resistivities.³⁰ As discussed in Refs. 20 and 21, the full MHD magnetic field can be written as (no approximations)

where the first term is the part of the magnetic field generated by the coils, the second is field line bending, and the last term mostly corresponds to field compression but also adds a small correction to field line bending. Since the vacuum field $∇χ$ must be divergence-free, it can be written in terms of Clebsch potentials: $∇χ=∇ψv×∇βv$⁠. As discussed in Ref. 20, one can construct a Clebsch-type coordinate system with coordinates $(ψv,βv,χ)$⁠; this coordinate system will be used further in this section. Expression (1) can be seen as just the plasma current-induced magnetic field (whose vector potential is $Ψ∇χ+Ω∇ψv$⁠, with the $∇βV$ component removed by a gauge transform) added to the coil field. The full MHD velocity can be written as (no approximations)

The terms approximately separate the MHD waves, with the first term containing Alfvén waves, the second containing slow magnetosonic waves, and the last one containing fast magnetosonic waves. Here, $Bv=|∇χ|$⁠. The reduced model is obtained by setting ζ = 0 and Ω = 0, and dropping the component of current perpendicular to $∇χ$ in Ohm's law. In addition to that, we perform a further reduction in this paper by setting $v∥=0$⁠. The removal of $v∥$ decreases the accuracy of the model and narrows the range of scenarios where it is valid^31,32 (note that the model tested in Ref. 30 has $v∥≠0$⁠.); however, the simplified model without $v∥$ is still applicable to the cases considered here, as will be seen. The resulting model consists of Eqs. (2.9), (2.13), (2.15), and (2.18) from Ref. 21, with ζ, Ω, and $v∥$ zeroed out,

In addition, $[f,g]=Bv−1∇χ·(∇f×∇g)$ is the Poisson bracket of two scalar functions f and g, and $(f,g)=∇⊥f·∇⊥g$ is the inner product of their gradients. In the reduced model, the ansatzes (1) and (2) become

In Eq. (3), η is the resistivity, $μ⊥$ is a viscosity-like parameter (it is not the same as physical dynamic viscosity, as discussed in Ref. 21), μ_h is the hyperviscosity, η_h is the hyper-resistivity (artificial dissipation parameters in the $Φ$ and $Ψ$ equations, respectively, that can help with numerical stabilization), $D⊥$ is mass diffusion across field lines, $κ⊥$ and $κ∥$ are the heat conduction coefficients across and along field lines, respectively, and $Sρ$ and S_e are the mass and energy sources, respectively. The Ohmic resistivity $ηOhm$ is a separate parameter, which allows one to neglect part of or all of the resistive contribution to internal energy and simply remove that energy from the system, which can be useful if the resistivity is artificially high.

Finally, there are two auxiliary variables: the normalized current in the $∇χ$ direction $j̃=−∇χ·j→/Bv2=−∇χ·∇×B→/(μ0Bv2)$ and the contravariant χ component of vorticity $ω̃=−∇χ·ω→=−∇χ·∇×v→$⁠, both taken with the opposite sign. The definition equations [Eqs. (3e) and (3f)] can be obtained by taking the dot product of $∇χ$ with the curl of $B→$ and $v→$⁠, respectively, and then using the identity $∇f·∇×Q→=−∇·(∇f×Q→)$⁠, where f is an arbitrary scalar field and $Q→$ is an arbitrary vector field. Instead of simply substituting the definition equations [Eqs. (3e) and (3f)] into the rest of the Eq. (3), $j̃$ and $ω̃$ were treated as separate variables, each with their own degrees of freedom on the finite element grid. These degrees of freedom are simultaneously evaluated at each time step with the degrees of freedom for the other variables by using the definition equations [Eqs. (3e) and (3f)] alongside with the rest of Eq. (3). This approach, used in combination with transforming the equations to weak form, allows one to avoid second-order derivatives in Eq. (3), except for the hyperviscosity term in Eq. (3a). Second-order derivatives can have discontinuities, as the finite elements in JOREK presently only have G¹ continuity.¹ This also means that if one tries to represent χ in the finite element basis instead of using the Dommaschk analytical form, one will not be able to avoid discontinuities in the first term on the RHS of Eq. (3a) even after applying integration by parts. In our experience, having discontinuities in the advective terms can decrease numerical accuracy or may lead to numerical instabilities, which was one of the reasons for the existence of $ω̃$ as a separate variable. A recent development in JOREK allows one to use basis functions with higher-order Gⁿ continuity, where n is a user-selected parameter.³³ This will allow one to eliminate $j̃$ and $ω̃$ as separate variables; however, it has not been ported to JOREK3D yet.

It is important to note that the $Φ$ evolution equation [Eq. (3a)] was obtained in Ref. 21 by applying a projection operator to the MHD momentum equation

and then inserting the ansatzes (4). The viscosity and hyperviscosity terms are separately added after the projection operator is applied (see Ref. 21 for more details). The projection operator that produces Eq. (3a) is

If the $v∥$ variable is kept in the reduced model, then the following projection operator produces the $v∥$ evolution equation (not shown here) when applied to Eq. (5):

The $v∥$ evolution equation is not included in the model (3) and was not used in any of the simulations presented in this paper, but it will be considered in the equilibrium error discussion, which comprises the remainder of this section.

A natural question that arises when considering reduced MHD is whether or not the reduction preserves equilibria. In other words, if a particular equilibrium solution to the full MHD equations is known, will it also be a solution to the reduced MHD equilibrium equations? As shown in Refs. 1 and 21, a simple argument involving the Grad–Shafranov equation shows that this is indeed the case in the tokamak limit, where the reduced MHD model (3) reduces to the model that was already used in the tokamak version of JOREK. However, this does not work for a general stellarator, where a full MHD equilibrium does not exactly satisfy the reduced MHD equilibrium equations, and a residual force arises and contributes to Eq. (3a). However, it can be shown that this contribution is small using an ordering argument.

Let $L⊥$ be the length scale perpendicular to $∇χ$ and $L∥$ be the length scale along $∇χ$⁠. Then, defining $λ≡L⊥/L∥$ as the ordering parameter, the spatial derivatives must satisfy $|∂∥|∼λ|∇⊥|$⁠. The terms in the full magnetic field (1) are ordered as follows:

and

where $Fv=|∇ψv|$⁠. Identifying $L⊥$⁠, B_v, and F_v as zeroth-order quantities, $L⊥=O(1), Bv=O(1), Fv=O(1)$⁠, it follows that $L∥=O(λ−1), ∇⊥=O(1), ∂∥=O(λ), Ψ=O(λ)$⁠, and $Ω=O(λ2)$⁠. Meanwhile, the residual force due to the reduction is $f→res=∇p−j→×B→=j→f×B→f−j→×B→$⁠, where $B→$ is the reduced MHD magnetic field (4), $B→f$ is the full magnetic field (1), and $j→$ and $j→f$ are the curls of the corresponding field divided by μ₀. Note that the residual force is just the difference between the full and reduced MHD Lorentz forces. It arises due to the neglect of the last term of (1) in reduced MHD and will be present even in the zero β limit. Inserting the ansatzes, the following expression is obtained for the residual force:

After some algebra, the reduced MHD current can be written as

where the Einstein summation convention is used, with $k,n∈{ψv,βv,χ}$ and $i∈{ψv,βv}$⁠. Here, g is the metric tensor of the Clebsch-type coordinate system aligned to $∇χ$⁠, which was introduced in Ref. 20, qⁱ represents the actual coordinates: $qi∈{ψv,βv}$⁠, and $e→ i$ are the contravariant basis vectors: $e→ i∈{∇ψv,∇βv}$⁠. With this, the first term in the residual force (8) expands to

Since $j→⊥=O(λ2)$⁠, it is easy to see that the first two terms above are $O(λ4)$ and the third term is $O(λ5)$⁠.

The second term in the residual force (8) can be expanded as

Note that $∇Ω×∇ψv=−(∂Ω/∂βv)∇χ+(∂Ω/∂χ)e→βv/J$⁠, where $e→βv=J∇χ×∇ψv$ is the covariant basis vector in the β_v direction in the Clebsch-type coordinate system aligned to $∇χ$⁠, and $J=[(∇ψv×∇βv)·∇χ]−1=1/Bv2$ is the Jacobian. Furthermore, since $Bv∂/∂χ=∂∥$⁠, one has $(∇Ω×∇ψv)ψv=gψvβvBv∂∥Ω$ and $(∇Ω×∇ψv)βv=gβvβvBv∂∥Ω$⁠. Thus, the first two terms in (10) are $O(λ4)$ and the third term is $O(λ2)$⁠. However, it is easy to see that the third term in (10) is in the kernel of the projection operator (6), and so this term will not contribute to the residual force in the $Φ$ evolution equation [Eq. (3a)]. On the other hand, the third term in (10) is not in the kernel of the projection operator (7), but its image under the operator, which can be written as $Bv3[∂Ω/∂βv,Ψ]$⁠, will be canceled by the image of another term, as will be shown below.

The curl of the last term of the full magnetic field (1) can be written as

Note that the first term in (11) is $O(λ2)$⁠, and the other two terms are $O(λ3)$⁠. Using (11), the third term in the residual force (8) becomes

Terms of the order $λ4$ are not written out explicitly in (12), since there is no need to consider them, as it is already established that there is at least an $O(λ4)$ contribution to equation (3a). As can be seen, the $O(λ3)$ term in (12) will be canceled by the projection operator (6) and as such will not contribute to the $Φ$ evolution equation [Eq. (3a)]; however, it will not be canceled by the projection operator (7). Indeed, its image under the operator (7) will be $Bv3[Ψ,∂Ω/∂βv]$⁠. Note that this image is equal to the negative image of the third term in (10), which was discussed above. These two images will cancel, and thus, the lowest order in which the residual force will contribute to the $v∥$ evolution equation is $λ4$⁠.

Finally, the last term in the residual force (8) is clearly of order $λ4$ or higher, so there is no need to consider it in detail like the other terms. In order to compare the residual force contributions to the other terms in Eq. (3a) and the $v∥$ evolution equation, some more ordering needs to be done. Consider that the shortest timescale in the reduced system is the Alfvén time $τA≡L∥/cA$⁠, where the parallel length scale is used because the Alfvén wave travels along field lines, and so the time derivative is ordered as $|∂/∂t|∼1/τA$⁠. As such, $∂/∂t=O(λ)$⁠. The $Φ$ and $v∥$ terms in the velocity ansatz are then ordered as

Assuming that the partial and convective terms in the material derivative are of the same order, $|∂/∂t|∼|v→·∇|$⁠, one has $Φ,v∥=O(λ)$⁠. After identifying $ρ=O(1)$ and $p=O(λ2)$⁠, it is clear that the lowest-order terms in Eq. (5) are $O(λ2)$⁠, and both projection operators (6) and (7) are O(1). As such, the lowest-order terms in both Eq. (3a) and the $v∥$ evolution equation are $O(λ2)$⁠. Thus, the residual force contribution to the reduced MHD equations is at least two orders of λ higher than the lowest-order terms. In Sec. V, it will be confirmed with numerical simulations that the residual force is indeed small.

Since χ is a solution of the Laplace equation in a torus, it can be represented as a summation over toroidal harmonics

where $F0ϕ$ corresponds to a tokamak-like toroidal field, n is the toroidal mode number, m is the poloidal mode number, and each harmonic individually satisfies the Laplace equation: $Δχn,m=0$⁠. Dommaschk gives a more explicit representation for χ,²² 

where a tilde denotes normalization: $χ=F0χ̃, R=R0R̃$ and $z=R0z̃$⁠; the normalization factor R₀ is the toroidally averaged radial position of the magnetic axis of the vacuum field. The functions $Dn,m$ and $Nn,m$ are defined as

and

The coefficients α_i, β_i, and γ_i are defined as

Although not written out explicitly, it can be seen that the coefficients also depend on n, the toroidal mode number of the D or N function that is being evaluated. The expressions above are only well defined if the following conditions on i and n are met: $i≥0$ for α_i and $αi*, i≥0$ and n > i for β_i and $βi*$⁠, and i > 0 for γ_i and $γi*$⁠. Otherwise, the corresponding coefficient and its starred version are zero. Finally, the coefficients $an,m, bn,m, cn,m$⁠, and $dn,m$ in Eq. (14) are what determine a particular configuration and must be calculated from the EXTENDER_P output.

Note that, since the harmonics $χn,m$ are analytically given, the property that $Δχn,m=0$ is exactly satisfied. This is an important advantage of using the Dommaschk representation for χ instead of the finite element representation (see  Appendix A), as it guarantees that the divergence-free condition on the magnetic field will be satisfied to machine precision. The second advantage is that χ and its derivatives are smooth.

EXTENDER_P provides the values of the three cylindrical components of the vacuum magnetic field, which will be referred to as $B→E$⁠, on an $(R,z,ϕ)$ grid. Setting $∇χ=B→E$ and considering the $ϕ$ component $BE,ϕ=ϕ̂·B→E=R−1∂χ/∂ϕ$⁠, one has the following:

We will make use of the properties [from Ref. 22, Eqs. (10) and (11)] $Dn,m|R̃=1=z̃m/m!$ and $Nn,m|R̃=1=0$⁠. Evaluating Eq. (18) at $R̃=1$ gives

If one also evaluates at $z̃=0$ and integrates over $ϕ$⁠, F₀ can be calculated

To calculate the coefficients $an,m$ and $bn,m$⁠, one must first multiply by either $sin nϕ$ or $cos nϕ$ and then use the orthogonality property of trigonometric functions

The number of terms M in the summations over m in Eq. (21) that is necessary to accurately represent the magnetic field is usually less than the number of poloidal modes used in the GVEC equilibrium. In practice, it is best to scan through different values of M, starting with the number of poloidal modes and decreasing from there while trying to minimize the error in $∇χ$ as compared to $B→E$⁠. Note that using higher values of M than necessary can lead to higher errors away from the $R̃=1$ surface due to overfitting, as the integration in Eqs. (22), (23), (26), and (27), which will be derived shortly, is only over the $R̃=1$ surface. Figure 1 shows the volume-averaged relative squared error of the Dommaschk potential representation as a function of the number M of poloidal modes kept in a Wendelstein 7-A equilibrium with $β=2.3×10−3 %$ (see Sec. V for more details about this equilibrium).

One can convert equations (21) into two linear algebraic systems with triangular matrices by changing the variable to $z′=z̃/Z$ and, after multiplying both equations by a Legendre polynomial $Pi(z′)$⁠, integrating from −1 to 1. Here, Z is determined as follows. In each poloidal plane at $R̃=1, z̃∈[−z̃−(ϕ),z̃+(ϕ)]$⁠, so $Z<minϕ{z̃−(ϕ),z̃+(ϕ)}$⁠. The value of Z is chosen to be slightly smaller than the minimum to avoid using the components of $B→E$ close to the boundary, where the output of EXTENDER_P can be less accurate. There is some freedom in choosing the specific value of Z, and it may take some trial and error to find the best value. As an example, Fig. 2 shows a segment of the surface of integration in one field period of the Wendelstein 7-A equilibria as described in Sec. V.

The Legendre polynomials are orthogonal to monomials of a lower order than the polynomial, since a monomial z^m can be expanded exactly in the Legendre polynomial basis of the same order m

Using the orthogonality property discussed above, one has, starting with i = M and descending to i = 0, the following linear algebraic system for $an,m$⁠:

where $⟨zi,Pj(z)⟩=∫−11ziPj(z)dz$⁠. Similarly, for the coefficient $bn,m$⁠, one has the following linear algebraic system:

As can be seen, both of these systems of equations have triangular matrices.

At this point, the equations for the coefficients $cn,m$ and $dn,m$ have yet to be determined. Consider now the R component of $B→E$⁠: $BE,R=R̂·B→E=∂χ/∂R$⁠. One has the following:

Again, evaluating at $R̃=1$ and using the properties $∂Dn,m/∂R̃|R̃=1=0$ and $∂Nn,m/∂R̃|R̃=1=z̃m/m!$ [also from Ref. 22, Eqs. (10) and (11)], one has the following:

From here, it is straightforward to follow the same steps as for $an,m$ and $bn,m$⁠, obtaining the following linear algebraic systems for $cn,m$⁠:

and for $dn,m$⁠:

Note that there are only M equations in each system for the unknowns $cn,1,…,cn,M$ and $dn,1,…,dn,M$ because $Nn,−1$ is not defined, and so terms with $cn,0$ and $dn,0$ are not included in the sum (14).

The only coefficients for which a system of equations has not yet been obtained are $a0,m$ (there are no $b0,m$ coefficients since $sin 0=0$⁠). These coefficients cannot be obtained from the system (22) since the matrices of this system are singular when n = 0. To get a solvable system, one must use the z component of $B→E$ as follows:

Evaluating at $R̃=1$ using the properties [from Ref. 22, Eqs. (10) and (11)] $Dn,m|R̃=1=z̃m/m!$ and $Nn,m|R̃=1=0$ after differentiating by $z̃$ gives

Integrating over $ϕ$ leaves only the n = 0 term in the sum, as all others are harmonic,

To finalize the derivation, we multiply the equation by a Legendre polynomial $Pi(z′)$ and integrate from -1 to 1. Starting from $i=M−1$ and descending to i = 0, the system of equations is

Just as in the case of systems (26) and (27), there are only M equations for the unknowns $a0,1,…,a0,M$⁠. This is because $D0,0=1$⁠, and so $a0,0$ is an additive constant in the scalar potential, which has no effect on the vacuum magnetic field.²² 

The linear algebraic systems of Eqs. (22), (23), (26), (27), and (28) are solved in a Python script using the NumPy library.³⁶ The solution is then written out to a Fortran namelist file, which can be read by JOREK. When evaluating the Dommaschk potential and its derivatives at any particular point, JOREK will then use the analytical representation (14).

As was mentioned in the Sec. II, although $j̃$ and $Ψ$ are related by $j̃=Δ*Ψ, j̃$ is stored as a separate variable in finite element representation for numerical purpose. It makes sense to first calculate the initial condition for $j̃$ from the GVEC data and then calculate $Ψ0$ from $j̃0$ using Eq. (3e) at t = 0,

Here, the subscript 0 refers to the fact that $Ψ0$ and $j̃0$ are the initial values of $Ψ$ and $j̃$⁠.

The equilibrium magnetic field provided by the GVEC solution will be referred to as $B→GVEC$⁠. Since GVEC works with full MHD, one needs to consider the full MHD ansatz, as given by (21), when working with $B→GVEC$⁠,

Taking the curl of the above equation and dotting it with $∇χ$⁠, after some algebra, one has the following:

Using the same ordering as in Sec. II, where $Bv=O(1), Ψ=O(λ),Ω=O(λ2)$ and $∂∥=O(λ)$⁠, it can be seen that the first term is $O(λ)$ and the second term is $O(λ3)$⁠. Thus, the second term can be neglected, due to being two orders of λ higher than the first term. This significantly simplifies the calculation, as now one can just set $j̃0=−jGVECχ/Bv2=−∇χ·∇×B→GVEC/Bv2$⁠.

Having determined $j̃0$⁠, it remains to solve Eq. (32) for $Ψ0$⁠. First, however, one needs to determine the boundary condition on $Ψ$⁠. Note that $Ψ$ is not constant on flux surfaces in stellarator geometry, unlike the tokamak situation. When running a fixed boundary simulation, as done in this study, it is usually assumed that the plasma is surrounded by a perfect conductor, so the magnetic field at the boundary does not have a normal component: $n→·B→=0$⁠. In the reduced MHD model, this means that $Ψ$ has to satisfy $n→·(∇Ψ×∇χ)=−n→·∇χ$ at all times. This is a nonhomogeneous linear differential equation, which must be solved on the boundary of the torus; the solution to this differential equation then provides a nonhomogeneous Dirichlet boundary condition for Eq. (32). Note that the kernel of the differential operator in the boundary equation is quite large, consisting of all functions $f(χ)$⁠. Using a flux-surface-aligned coordinate system $(ψ,θ,ϕ)$⁠, where ψ is a flux surface label, θ is the GVEC poloidal angle, which, like in VMEC, is constructed for each particular equilibrium in the course of minimization, and $ϕ$ is the geometric toroidal angle, the boundary equation becomes

where $J=[∇ψ·(∇θ×∇ϕ)]−1$ is the Jacobian. Note that $(ψ,θ,ϕ)$ is not a straight field line coordinate system. However, solving the equation in this form is numerically difficult because one cannot easily separate the kernel and remove it from the solution space. To do so, one must switch to a coordinate system where χ is one of the coordinates. It is best to switch out $ϕ$ for χ, since a stellarator must have a nonvanishing toroidal component to its vacuum field [c.f. the $F0ϕ$ term in Eq. (13)], so $∂χ/∂ϕ$ is nonvanishing and the Jacobian of the new coordinates is nowhere singular. The boundary equation in $(ψ,θ,χ)$ coordinates is

where $J′=[∇ψ·(∇θ×∇χ)]−1$ is the new Jacobian. It is easy to solve this equation in JOREK. Due to the JOREK grid for our applications here being flux-surface-aligned, the element local coordinates s and t (see  Appendix A) can be related to the coordinates ψ and θ as $ψ=s, θ=2π(t+ibnd_elm)/Nbnd_elm$⁠, where $ibnd_elm$ is the zero-based index of the current boundary element, and $Nbnd_elm$ is the total number of boundary elements. Finally, the χ coordinate is given by the Dommaschk representation (14). The solution space in which the solution to Eq. (34) is searched for can now be represented as

where $Ncp$ is defined in  Appendix A, and $χ̃=χ/F0$⁠. Excluding the m = 0 mode removes the kernel of the differential operator of Eq. (34) from $Vsol$⁠, and the equation can then be solved using the standard Fourier–Galerkin method. The solution obtained this way is then projected back onto the JOREK finite element basis and written to the boundary nodes. Finally, Eq. (32) is solved by splitting $Ψ0=Ψ0,i+Ψb$⁠, where $Ψb$ is the solution to Eq. (34) and thus satisfies the nonhomogeneous Dirichlet boundary condition, while $Ψ0,i$ is an unknown function, which is zero at the boundary. The solution $Ψ0,i$ is then found using the standard JOREK solver with homogeneous Dirichlet boundary conditions. When $Ψ$ is evolved in time, JOREK will solve for the increment $δΨ$ at each time step, which must also be zero at the boundary (and thus can also be obtained using the standard solver with homogeneous Dirichlet boundary conditions), so that the nonhomogeneous boundary condition continues to be satisfied for the total $Ψ$⁠.

The last step is determining an initial condition for temperature, which is almost trivial. The GVEC solution provides a pressure profile $pGVEC$⁠, which must be simply converted to JOREK units and divided by the initial density profile ρ₀. In all of the stellarator simulations presented in this paper, the initial density is taken to be constant for simplicity, which corresponds to $ρ0=1$ in JOREK units.