Although the basic concept of a stellarator was known since the early days of fusion research, advances in computational technology have enabled the modeling of increasingly complicated devices, leading up to the construction of Wendelstein 7-X, which has recently shown promising results. This recent success has revived interest in the nonlinear 3D MHD modeling of stellarators in order to better understand their performance and operational limits. This study reports on the extension of the JOREK code to 3D geometries and on the first stellarator simulations carried out with it. The first simple simulations shown here address the classic Wendelstein 7-A stellarator using a reduced MHD model previously derived by us. The results demonstrate that stable full MHD equilibria are preserved in the reduced model: the flux surfaces do not move throughout the simulation and closely match the flux surfaces of the full MHD equilibrium. Furthermore, both tearing and ballooning modes were simulated, and the linear growth rates measured in JOREK are in reasonable agreement with the growth rates from the CASTOR3D linear MHD code.

## I. INTRODUCTION

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.^{2}

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,^{5} 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*^{1},^{10} M3D,^{8} and MIPS,^{11} 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,^{12} although the tokamak capabilities of NIMROD are already enough to simulate a stellarator with an axisymmetric vacuum vessel.^{13} Also, the FLUXO nonlinear MHD code^{14} 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 $\chi =F0\varphi $, where $\varphi $ 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 $\u2207\chi $ is a magnetic field, it must be divergence-free, so $\Delta \chi =0$. One then needs to find a general solution to the Laplace equation in the toroidal coordinate system $(R,z,\varphi )$. This was done by Dommaschk,^{22} 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,\varphi )$ grid, which we do using the EXTENDER_P code.^{23}

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:

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

^{24}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.

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

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

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

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

The GVEC^{24} 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^{27} 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.

## II. REDUCED MHD AND THE RESIDUAL FORCE

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.^{30} 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 $\u2207\chi $ must be divergence-free, it can be written in terms of Clebsch potentials: $\u2207\chi =\u2207\psi v\xd7\u2207\beta v$. As discussed in Ref. 20, one can construct a Clebsch-type coordinate system with coordinates $(\psi v,\beta v,\chi )$; 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 $\Psi \u2207\chi +\Omega \u2207\psi v$, with the $\u2207\beta 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=|\u2207\chi |$. The reduced model is obtained by setting *ζ* = 0 and Ω = 0, and dropping the component of current perpendicular to $\u2207\chi $ in Ohm's law. In addition to that, we perform a further reduction in this paper by setting $v\u2225=0$. The removal of $v\u2225$ 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\u2225\u22600$.); however, the simplified model without $v\u2225$ 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\u2225$ zeroed out,

In addition, $[f,g]=Bv\u22121\u2207\chi \xb7(\u2207f\xd7\u2207g)$ is the Poisson bracket of two scalar functions *f* and *g*, and $(f,g)=\u2207\u22a5f\xb7\u2207\u22a5g$ is the inner product of their gradients. In the reduced model, the ansatzes (1) and (2) become

In Eq. (3), *η* is the resistivity, $\mu \u22a5$ is a viscosity-like parameter (it is not the same as physical dynamic viscosity, as discussed in Ref. 21), *μ _{h}* is the hyperviscosity,

*η*is the hyper-resistivity (artificial dissipation parameters in the $\Phi $ and $\Psi $ equations, respectively, that can help with numerical stabilization), $D\u22a5$ is mass diffusion across field lines, $\kappa \u22a5$ and $\kappa \u2225$ are the heat conduction coefficients across and along field lines, respectively, and $S\rho $ and

_{h}*S*are the mass and energy sources, respectively. The Ohmic resistivity $\eta 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.

_{e}Finally, there are two auxiliary variables: the normalized current in the $\u2207\chi $ direction $j\u0303=\u2212\u2207\chi \xb7j\u2192/Bv2=\u2212\u2207\chi \xb7\u2207\xd7B\u2192/(\mu 0Bv2)$ and the contravariant *χ* component of vorticity $\omega \u0303=\u2212\u2207\chi \xb7\omega \u2192=\u2212\u2207\chi \xb7\u2207\xd7v\u2192$, both taken with the opposite sign. The definition equations [Eqs. (3e) and (3f)] can be obtained by taking the dot product of $\u2207\chi $ with the curl of $B\u2192$ and $v\u2192$, respectively, and then using the identity $\u2207f\xb7\u2207\xd7Q\u2192=\u2212\u2207\xb7(\u2207f\xd7Q\u2192)$, where *f* is an arbitrary scalar field and $Q\u2192$ is an arbitrary vector field. Instead of simply substituting the definition equations [Eqs. (3e) and (3f)] into the rest of the Eq. (3), $j\u0303$ and $\omega \u0303$ 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*^{1} continuity.^{1} 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 $\omega \u0303$ as a separate variable. A recent development in JOREK allows one to use basis functions with higher-order *G ^{n}* continuity, where

*n*is a user-selected parameter.

^{33}This will allow one to eliminate $j\u0303$ and $\omega \u0303$ as separate variables; however, it has not been ported to JOREK3D yet.

It is important to note that the $\Phi $ 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\u2225$ variable is kept in the reduced model, then the following projection operator produces the $v\u2225$ evolution equation (not shown here) when applied to Eq. (5):

The $v\u2225$ 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\u22a5$ be the length scale perpendicular to $\u2207\chi $ and $L\u2225$ be the length scale along $\u2207\chi $. Then, defining $\lambda \u2261L\u22a5/L\u2225$ as the ordering parameter, the spatial derivatives must satisfy $|\u2202\u2225|\u223c\lambda |\u2207\u22a5|$. The terms in the full magnetic field (1) are ordered as follows:

and

where $Fv=|\u2207\psi v|$. Identifying $L\u22a5$, *B _{v}*, and

*F*as zeroth-order quantities, $L\u22a5=O(1),\u2009Bv=O(1),\u2009Fv=O(1)$, it follows that $L\u2225=O(\lambda \u22121),\u2009\u2207\u22a5=O(1),\u2009\u2202\u2225=O(\lambda ),\u2009\Psi =O(\lambda )$, and $\Omega =O(\lambda 2)$. Meanwhile, the residual force due to the reduction is $f\u2192res=\u2207p\u2212j\u2192\xd7B\u2192=j\u2192f\xd7B\u2192f\u2212j\u2192\xd7B\u2192$, where $B\u2192$ is the reduced MHD magnetic field (4), $B\u2192f$ is the full magnetic field (1), and $j\u2192$ and $j\u2192f$ are the curls of the corresponding field divided by

_{v}*μ*

_{0}. 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\u2208{\psi v,\beta v,\chi}$ and $i\u2208{\psi v,\beta v}$. Here, *g* is the metric tensor of the Clebsch-type coordinate system aligned to $\u2207\chi $, which was introduced in Ref. 20, *q ^{i}* represents the actual coordinates: $qi\u2208{\psi v,\beta v}$, and $e\u2192\u2009i$ are the contravariant basis vectors: $e\u2192\u2009i\u2208{\u2207\psi v,\u2207\beta v}$. With this, the first term in the residual force (8) expands to

Since $j\u2192\u22a5=O(\lambda 2)$, it is easy to see that the first two terms above are $O(\lambda 4)$ and the third term is $O(\lambda 5)$.

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

Note that $\u2207\Omega \xd7\u2207\psi v=\u2212(\u2202\Omega /\u2202\beta v)\u2207\chi +(\u2202\Omega /\u2202\chi )e\u2192\beta v/J$, where $e\u2192\beta v=J\u2207\chi \xd7\u2207\psi v$ is the covariant basis vector in the *β _{v}* direction in the Clebsch-type coordinate system aligned to $\u2207\chi $, and $J=[(\u2207\psi v\xd7\u2207\beta v)\xb7\u2207\chi ]\u22121=1/Bv2$ is the Jacobian. Furthermore, since $Bv\u2202/\u2202\chi =\u2202\u2225$, one has $(\u2207\Omega \xd7\u2207\psi v)\psi v=g\psi v\beta vBv\u2202\u2225\Omega $ and $(\u2207\Omega \xd7\u2207\psi v)\beta v=g\beta v\beta vBv\u2202\u2225\Omega $. Thus, the first two terms in (10) are $O(\lambda 4)$ and the third term is $O(\lambda 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 $\Phi $ 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[\u2202\Omega /\u2202\beta v,\Psi ]$, 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(\lambda 2)$, and the other two terms are $O(\lambda 3)$. Using (11), the third term in the residual force (8) becomes

Terms of the order $\lambda 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(\lambda 4)$ contribution to equation (3a). As can be seen, the $O(\lambda 3)$ term in (12) will be canceled by the projection operator (6) and as such will not contribute to the $\Phi $ 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[\Psi ,\u2202\Omega /\u2202\beta 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\u2225$ evolution equation is $\lambda 4$.

Finally, the last term in the residual force (8) is clearly of order $\lambda 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\u2225$ evolution equation, some more ordering needs to be done. Consider that the shortest timescale in the reduced system is the Alfvén time $\tau A\u2261L\u2225/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 $|\u2202/\u2202t|\u223c1/\tau A$. As such, $\u2202/\u2202t=O(\lambda )$. The $\Phi $ and $v\u2225$ 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, $|\u2202/\u2202t|\u223c|v\u2192\xb7\u2207|$, one has $\Phi ,v\u2225=O(\lambda )$. After identifying $\rho =O(1)$ and $p=O(\lambda 2)$, it is clear that the lowest-order terms in Eq. (5) are $O(\lambda 2)$, and both projection operators (6) and (7) are *O*(1). As such, the lowest-order terms in both Eq. (3a) and the $v\u2225$ evolution equation are $O(\lambda 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.

## III. FINDING THE DOMMASCHK REPRESENTATION OF A SCALAR POTENTIAL

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

where $F0\varphi $ 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: $\Delta \chi n,m=0$. Dommaschk gives a more explicit representation for *χ*,^{22}

where a tilde denotes normalization: $\chi =F0\chi \u0303,\u2009R=R0R\u0303$ and $z=R0z\u0303$; the normalization factor *R*_{0} 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}*,

*β*, and

_{i}*γ*are defined as

_{i}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\u22650$ for *α _{i}* and $\alpha i*,\u2009i\u22650$ and

*n*>

*i*for

*β*and $\beta i*$, and

_{i}*i*> 0 for

*γ*and $\gamma i*$. Otherwise, the corresponding coefficient and its starred version are zero. Finally, the coefficients $an,m,\u2009bn,m,\u2009cn,m$, and $dn,m$ in Eq. (14) are what determine a particular configuration and must be calculated from the EXTENDER_P output.

_{i}Note that, since the harmonics $\chi n,m$ are analytically given, the property that $\Delta \chi 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\u2192E$, on an $(R,z,\varphi )$ grid. Setting $\u2207\chi =B\u2192E$ and considering the $\varphi $ component $BE,\varphi =\varphi \u0302\xb7B\u2192E=R\u22121\u2202\chi /\u2202\varphi $, one has the following:

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

If one also evaluates at $z\u0303=0$ and integrates over $\varphi $, *F*_{0} can be calculated

To calculate the coefficients $an,m$ and $bn,m$, one must first multiply by either $sin\u2009n\varphi $ or $cos\u2009n\varphi $ 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 $\u2207\chi $ as compared to $B\u2192E$. Note that using higher values of *M* than necessary can lead to higher errors away from the $R\u0303=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\u0303=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 $\beta =2.3\xd710\u22123\u2009%$ (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\u2032=z\u0303/Z$ and, after multiplying both equations by a Legendre polynomial $Pi(z\u2032)$, integrating from −1 to 1. Here, *Z* is determined as follows. In each poloidal plane at $R\u0303=1,\u2009z\u0303\u2208[\u2212z\u0303\u2212(\varphi ),z\u0303+(\varphi )]$, so $Z<min\varphi {z\u0303\u2212(\varphi ),z\u0303+(\varphi )}$. The value of *Z* is chosen to be slightly smaller than the minimum to avoid using the components of $B\u2192E$ 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 $\u27e8zi,Pj(z)\u27e9=\u222b\u221211ziPj(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\u2192E$: $BE,R=R\u0302\xb7B\u2192E=\u2202\chi /\u2202R$. One has the following:

Again, evaluating at $R\u0303=1$ and using the properties $\u2202Dn,m/\u2202R\u0303|R\u0303=1=0$ and $\u2202Nn,m/\u2202R\u0303|R\u0303=1=z\u0303m/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,\u2026,cn,M$ and $dn,1,\u2026,dn,M$ because $Nn,\u22121$ 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\u20090=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\u2192E$ as follows:

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

Integrating over $\varphi $ 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\u2032)$ and integrate from -1 to 1. Starting from $i=M\u22121$ 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,\u2026,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.^{22}

The linear algebraic systems of Eqs. (22), (23), (26), (27), and (28) are solved in a Python script using the NumPy library.^{36} 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).

## IV. DETERMINING INITIAL CONDITIONS FROM THE GVEC SOLUTION

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

Here, the subscript 0 refers to the fact that $\Psi 0$ and $j\u03030$ are the initial values of $\Psi $ and $j\u0303$.

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

Taking the curl of the above equation and dotting it with $\u2207\chi $, after some algebra, one has the following:

Using the same ordering as in Sec. II, where $Bv=O(1),\u2009\Psi =O(\lambda ),\Omega =O(\lambda 2)$ and $\u2202\u2225=O(\lambda )$, it can be seen that the first term is $O(\lambda )$ and the second term is $O(\lambda 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\u03030=\u2212jGVEC\chi /Bv2=\u2212\u2207\chi \xb7\u2207\xd7B\u2192GVEC/Bv2$.

Having determined $j\u03030$, it remains to solve Eq. (32) for $\Psi 0$. First, however, one needs to determine the boundary condition on $\Psi $. Note that $\Psi $ 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\u2192\xb7B\u2192=0$. In the reduced MHD model, this means that $\Psi $ has to satisfy $n\u2192\xb7(\u2207\Psi \xd7\u2207\chi )=\u2212n\u2192\xb7\u2207\chi $ 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(\chi )$. Using a flux-surface-aligned coordinate system $(\psi ,\theta ,\varphi )$, 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 $\varphi $ is the geometric toroidal angle, the boundary equation becomes

where $J=[\u2207\psi \xb7(\u2207\theta \xd7\u2207\varphi )]\u22121$ is the Jacobian. Note that $(\psi ,\theta ,\varphi )$ 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 $\varphi $ for *χ*, since a stellarator must have a nonvanishing toroidal component to its vacuum field [c.f. the $F0\varphi $ term in Eq. (13)], so $\u2202\chi /\u2202\varphi $ is nonvanishing and the Jacobian of the new coordinates is nowhere singular. The boundary equation in $(\psi ,\theta ,\chi )$ coordinates is

where $J\u2032=[\u2207\psi \xb7(\u2207\theta \xd7\u2207\chi )]\u22121$ 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 $\psi =s,\u2009\theta =2\pi (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 $\chi \u0303=\chi /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 $\Psi 0=\Psi 0,i+\Psi b$, where $\Psi b$ is the solution to Eq. (34) and thus satisfies the nonhomogeneous Dirichlet boundary condition, while $\Psi 0,i$ is an unknown function, which is zero at the boundary. The solution $\Psi 0,i$ is then found using the standard JOREK solver with homogeneous Dirichlet boundary conditions. When $\Psi $ is evolved in time, JOREK will solve for the increment $\delta \Psi $ 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 $\Psi $.

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 *ρ*_{0}. In all of the stellarator simulations presented in this paper, the initial density is taken to be constant for simplicity, which corresponds to $\rho 0=1$ in JOREK units.

## V. A CONSISTENCY CHECK FOR THE STELLARATOR MODEL

After having derived and implemented the stellarator model, it remains to validate it for stellarators, showing that it does work. However, before proceeding to more complicated cases, a set of initial tests must be performed using stable equilibria to demonstrate that the model is indeed consistent, the error due to the neglect of fourth-order terms in the Lorentz force, which was discussed in Sec. II, is small, and no significant change is observed in the stable cases after simulating them for some time.

The consistency checks were done using four equilibria based on the historic Wendelstein 7-A stellarator^{27} with different values of *β*. Note that the *β* values here and throughout the rest of this paper are volume-averaged. These equilibria were intended to be unstable to the (2,1) tearing mode; however, since Wendelstein 7-A had five field periods, the simulations can be performed with fivefold periodicity, excluding the unstable *n* = 1 Fourier mode and its mode family. Thus, in the computational setting used in this section, there are no physical instabilities. The equilibria were first calculated with NEMEC,^{37} and then, GVEC was used to refine them. Poloidal modes $m=0,\u2026,12$ and toroidal modes $n=0,5,10,\u2026,50$, which correspond to $Nctor=21$ and $Ncp=5$ in JOREK (see Appendix A), were used to calculate the equilibrium. All of the equilibria have the same boundary: a rotating ellipse with a minor axis of 0.09131 m and a major axis of 0.1178 m; the major radius of the torus is 1.99 m. The normalized toroidal current profile was also the same for all equilibria,

where $\psi tn$ is the toroidal flux normalized so that $\psi tn=0$ at the axis and $\psi tn=1$ at the boundary. $In(\psi tn)$, which represents the toroidal current enclosed by the flux surface $\psi tn$, is normalized by the total toroidal current, which was 17.5 kA in the cases considered, such that $In(1)=1$. The total toroidal magnetic flux through a poloidal plane was 0.08 Wb. The pressures at the axis were 1, 100, and 500 Pa, and 1 kPa, which correspond to the *β*-values of $2.3\xd710\u22125\u2009%,\u20092.3\xd710\u22123\u2009%$, 0.011%, and 0.022%, respectively. The pressure profiles are given by

where *p* is the pressure in pascals, $pa$ is the pressure at the axis, and $pb$ is the pressure at the boundary. For the $\beta =2.3\xd710\u22125\u2009%$ ($pa=1$ Pa) case, $pb=0.01$ Pa, while for the other three cases, $pb=1$ Pa. When finding the initial conditions from the GVEC equilibrium, $Ntor=9$ and $Np=5$ (see Appendix A) were used for the variables, which correspond to Fourier modes $n=0,5,\u2026,20$. The profile of the rotational transform *ι*, which was the same for all equilibria considered in this study, is shown in Fig. 3.

All of the simulations were run with a spatially constant resistivity $\eta =1.938\u200910\u22129\u2009\u2009\Omega \u2009m$ and viscosity $\mu =2.90\xd710\u22129\u2009kg/(m\u2009s)$. In addition, a hyperviscosity $\mu h=2.90\xd710\u221212\u2009\u2009kg\u2009m/s$ was applied. The radial resolution of the finite elements was 41 nodes, and the poloidal resolution was 48 nodes, that is, moving from the axis to the edge, 41 grid nodes will be passed, counting the nodes at the axis and edge, and 48 grid nodes will be passed in one poloidal turn. The first part of the simulation was run using the implicit Euler time-stepping scheme to damp out small oscillations that were present due to the neglect of fourth-order terms in the Lorentz force and different discrete representations in JOREK and GVEC (see Fig. 4). This consisted of 20 time steps of length $6.484\xd710\u22124\u2009ms$ (1 in JOREK units), followed by 20 time steps of length $6.484\xd710\u22123\u2009ms$ (10 in JOREK units), followed by 10 time steps of length $6.484\xd710\u22122\u2009ms$ (100 in JOREK units). For the $\beta =0.022%$ case, but not for the others, this was followed by another 10 time steps of length $6.484\xb710\u22122\u2009ms$. In the second part of the simulation, the Crank–Nicolson time-stepping scheme was used,^{38} and all four cases were simulated for 6.484 ms (10 000 in JOREK units). The $\beta =2.3\xd710\u22125\u2009%$ and $\beta =2.3\xd710\u22123\u2009%$ cases used time steps of length $6.484\xd710\u22122\u2009ms$ in the second part; however, the $\beta =0.011\u2009%$ and $\beta =0.022\u2009%$ required shorter time steps ($3.242\xd710\u22122$ and $1.621\xd710\u22122\u2009ms$, respectively) for numerical stability. When evaluating the integrals in the weak form of Eq. (3), the toroidal integration was performed by summing over 40 poloidal planes evenly spread over one period.

As expected, no large-scale motion was observed in any of the four simulations, confirming in fact that equilibrium is maintained in the reduced MHD model and its implementation. This can be seen in Fig. 5, where the *R* coordinate of the magnetic axis is plotted as a function of time for each of the four simulations, along with the error bars. The straight line of the same color outside the error bars represents the axis position in the full MHD equilibrium as calculated in GVEC. The difference between the JOREK and GVEC axis positions is due to the magnetic field being approximated by the reduced MHD ansatz. The axis was determined by making an initial guess for its (*R*, *z*) position in the $\varphi =0$ poloidal plane and then tracing the field line at that position for ten toroidal turns, after which the tolerance $T=0.1(max\u2009Ri\u2212min\u2009Ri)(max\u2009zi\u2212min\u2009zi)$, where $i=1,\u2026,10$, is calculated. If this tolerance is smaller than the cutoff, which was set to $5\xd710\u22125\u2009m$, then the axis is considered found: the axis position at $\varphi =0$ is $(Rc,zc)=((max\u2009Ri+min\u2009Ri)/2,(max\u2009zi+min\u2009zi)/2)$. If not, then the field line tracing is restarted at $(Rc,zc)$, and the process is repeated until the tolerance is less than the cutoff. The error in the *R*-coordinate was estimated as $ER=0.1(max\u2009Ri\u2212min\u2009Ri)$ and plotted as error bars in Fig. 5.

To demonstrate that there is no significant motion even away from the axis, the Poincare plots for the $\beta =0.022\u2009%$ case are shown in Fig. 6, both before and after the simulation, along with the flux surfaces of the GVEC equilibrium. As can be seen, the flux surfaces in JOREK coincide with the GVEC flux surfaces, so the error introduced by using the reduced MHD ansatz for the magnetic field has no noticeable effect on the flux surfaces. Moreover, the flux surfaces do not move during the simulation, preserving the stable equilibrium as expected.

## VI. TEARING MODE BENCHMARK

Having demonstrated that basic stellarator simulations can be run with the correct equilibrium in the newly implemented model in JOREK, the next step is to simulate instabilities and benchmark them against known results. Tearing modes in the Wendelstein 7-A stellarator were used for this purpose. Three cases at different values of *β* were considered: $2.3\xd710\u22125\u2009%,\u20092.3\xd710\u22124\u2009%$, and $2.3\xd710\u22123\u2009%$. For reference, Fig. 7 shows the velocity stream function $\Phi $ without the *n* = 0, 5, 10 Fourier modes, which do not contribute to the tearing mode, on the $\varphi =0$ poloidal plane during the presaturation (linear) phase of the full-torus simulation (see below). The characteristic (2,1) structure of the mode is clearly visible. The $\beta =2.3\xd710\u22125\u2009%$ and $\beta =2.3\xd710\u22123\u2009%$ are the same equilibria that were used in the previous section, with the $\beta =2.3\xd710\u22124\u2009%$ being a new equilibrium with the same boundary and current profile as the other two and an intermediate value of *β*. In this new intermediate equilibrium, $pa=10$ Pa and $pb=0.1$ Pa. When finding the initial conditions from the GVEC equilibrium, $Ntor=5$ and $Np=5$ (see Appendix A) were used for the variables, which correspond to Fourier modes *n* = 0, 5, and 10.

Just as before, the stellarator simulations were run with the implicit Euler time-stepping scheme and fivefold periodicity to damp out oscillations. This consisted of 20 time steps of length $6.484\xd710\u22124\u2009\u2009ms$, followed by 20 time steps of length $6.484\xd710\u22123\u2009\u2009ms$, followed by 5 time steps of length $6.484\xd710\u22122\u2009ms$. The resistivity was set to $\eta =1.938\xd710\u22126\u2009\u2009\Omega \u2009m$, and the viscosity was zero. The hyperviscosity was $\mu h=2.90\xd710\u221215\u2009\u2009kg\u2009m/s$ for the $\beta =2.3\xd710\u22125\u2009%$, and the $\beta =2.3\xd710\u22124\u2009%$ cases, and $\mu h=7.25\xd710\u221215\u2009\u2009kg\u2009m/s$ for the $\beta =2.3\xd710\u22123\u2009%$ case. The finite element resolution was 41 nodes radially and 48 nodes poloidally, just as before. Both heat conduction and mass diffusion were set to zero. It should be noted that, when using $Ntor=5$, anisotropic transport cannot be properly modeled, as field lines tend to slightly drift from flux surfaces after many toroidal turns, which leads to parallel transport contributing to perpendicular transport after enough time steps. This problem can be remedied by including more toroidal modes.

In the second part of the simulation, the domain was extended to the full torus, taking now into account all of the $n=0,\u2026,10$ Fourier modes, corresponding to $Ntor=21$ and $Np=1$ (see Appendix A). The Crank–Nicolson scheme was used with time steps of length $1.621\xd710\u22122\u2009ms$. The toroidal integration in the weak form of equations (3) was performed by summing over 40 poloidal planes evenly spread throughout the full torus. The number of Fourier modes, the number of poloidal planes, and the values of hyperviscosity, resolution, and time step size were chosen after scanning over several values for each parameter and choosing the value at which the growth rate of the tearing mode converged. For the present purposes, convergence is considered to be achieved when halving the time step size or hyperviscosity, or doubling the resolution, the number of modes, or the number of planes leads to a change in the growth rate of less than 1.5%. The convergence test was performed for the $\beta =2.3\xd710\u22123\u2009%$ and $\beta =2.3\xd710\u22125\u2009%$ cases, resulting in all of the parameters converging to the same values, except for hyperviscosity, which converged to $\mu h=7.25\xd710\u221215\u2009kg\u2009m/s$ for the $\beta =2.3\xd710\u22123\u2009%$ case and $\mu h=2.90\xd710\u221215\u2009kg\xb7m/s$ for the $\beta =2.3\xd710\u22125\u2009%$ case. The $\beta =2.3\xd710\u22124\u2009%$ case was then run using the lower value of hyperviscosity. Figure 8(a) shows the values of the growth rates from JOREK alongside the values calculated in the linear MHD code CASTOR3D.^{28,29} The maximum deviation between the two codes is 13.2% and occurs at $\beta =2.3\xd710\u22124\u2009%$.

For reference, Fig. 9 shows the magnetic energies of each individual Fourier mode in the $\beta =2.3\xd710\u22125\u2009%$ case, except for modes that belong to the *n* = 0 mode family (*n* = 0, 5, 10). The time axis starts slightly before 0.5 ms, as that is where the full-torus simulation begins and the Fourier modes shown are initialized. As expected, the *n* = 1 Fourier mode drives the instability, with the rest of its mode family (*n* = 4, 6, and 9) growing with it due to linear mode coupling. Around 2 ms, the *n* = 1 mode begins to drive the *n* = 2 mode via nonlinear coupling, which in turn drives the rest of its mode family (*n* = 3, 7, and 8) via linear coupling. Finally, the mode saturates around 3–3.5 ms. The saturated magnetic island structure is shown in the Poincare plot in Fig. 10. Note that, aside from the dominant (2,1) island chain, there is a secondary (3,2) island chain toward the interior of the plasma, which is nonlinearly excited by the mode.

Two more simulations were performed with the $\beta =2.3\xd710\u22125\u2009%$ case, this time using resistivities of $\eta =1.938\xd710\u22127$ and $\eta =1.938\xd710\u22125\u2009\u2009\Omega \u2009m$, while all of the other parameters were kept the same as before. For the $\eta =1.938\xd710\u22127\u2009\u2009\Omega \u2009m$ case, a hyper-resistivity of $\eta h=9.691\xd710\u221214\u2009\u2009\Omega \u2009m2$ ($5\xd710\u221214$ in JOREK units) had to be introduced in order for the iterative solver to converge in a reasonable amount of time. However, it was first confirmed that introducing this amount of hyper-resistivity in the $\eta =1.938\xd710\u22126\u2009\u2009\Omega \u2009m$ case, which could be run with or without hyper-resistivity, changes the growth rate by less than 1.5%. Figure 8(b) shows the growth rates for the $\beta =2.3\xd710\u22125\u2009%$ case at different values of resistivity alongside the growth rates calculated by CASTOR3D. The maximum deviation between the two codes is 26.2%, occurring at $\eta =1.938\xd710\u22125\u2009\Omega \u2009m$. This is most likely due to the neglect of $v\u2225$ by the model used in these simulations, as $v\u2225$ can be large within the resistive layer, and the size of the resistive layer increases with resistivity. A similar effect is known to exist for quasi-interchange modes. Neglecting the parallel velocity leads to an overestimation of the quasi-interchange growth rates by a factor of $1+2qs2$, where *q _{s}* is the safety factor of the flux surface where the mode appears.

^{39}In general, the agreement on the growth rates for the (2,1) tearing mode looks convincing, with deviations on the order of 10% from CASTOR3D, which solves the linearized full MHD equations.

## VII. BALLOONING MODE BENCHMARK

A similar benchmark with CASTOR3D for ballooning mode growth rates in Wendelstein 7-A using equilibria with $\beta =0.11%$ and $\beta =0.21%$ (corresponding to axis pressures of 5 and 10 kPa) has been performed at resistivities of $\eta =1.938\xd710\u22127\u2009\u2009\Omega \u2009m$ and $\eta =5.814\xd710\u22127\u2009\Omega \u2009m$. For reference, the velocity stream function $\Phi $ is shown in Fig. 11 during the linear phase of the $\beta =0.21%,\u2009\eta =1.938\xd710\u22127\u2009\Omega \u2009m$ simulation. These equilibria have the same toroidal current and pressure profiles as the equilibria in Sec. V, with $pb$ set to 1 Pa and 100 Pa, respectively. As before, the simulations were initially run with the implicit Euler scheme to damp out oscillations. This first phase of the simulation consisted of 30 time steps of length $6.484\xd710\u22124\u2009ms$. In the second part of the simulation, the Crank–Nicolson scheme was used; however, both parts were run with a fivefold periodicity, since ballooning modes can be simulated with just one period.

Based on linear ballooning mode theory in the tokamak limit, ballooning mode growth rates are known to diverge as $n\u2192\u221e$ if there are no background flows present in the initial equilibrium.^{4} In order to realistically model ballooning modes, one would have to include an equilibrium flow, such as diamagnetic drift, which will stabilize modes with $n>nmax$.^{40} However, since we only want to do a benchmark, we do not use equilibrium flows in either JOREK or CASTOR3D, but simply cut off the number of Fourier modes in the JOREK simulations, and then limit the CASTOR3D run to the same number of modes. We chose to keep only the *n* = 0, 5, 10 modes in the simulations considered here; the highest mode (*n* = 10) is dominant in this case, as it grows the fastest, and its energy quickly exceeds that of the other modes. We have also tried including the *n* = 15 and *n* = 20 modes in case the prediction about increasing growth rates from the tokamak limit is no longer valid for Wendelstein 7-A; however, we found that the highest-*n* mode is the fastest-growing one in those simulations as well.

Holding the number of modes fixed at $Ntor=5$, a convergence test was performed for the $\beta =0.21%$ equilibrium while using a resistivity of $\eta =1.938\xd710\u22127\u2009\u2009\Omega \u2009m$. The growth rate converged at a finite element resolution of 61 nodes radially and 72 nodes poloidally, time step size of $6.484\xd710\u22124\u2009\u2009ms$, and a hyperviscosity of $\mu h=7.25\xd710\u221215\u2009\u2009kg\u2009m/s$. The other three simulations were then run with these parameters. Both heat conductivity and mass diffusion parameters were set to zero in all four simulations. Figure 12 shows the growth rates measured in the four cases, along with the growth rates calculated in CASTOR3D for the same four cases. The maximum deviation is 6.4% and occurs at $\beta =0.11%,\u2009\eta =1.938\xd710\u22127\u2009\u2009\Omega \u2009m$, with the other three deviations all being less than 3%.

Cutting off the Fourier series without stabilizing the high-*n* modes could lead to spectral blocking,^{41} where existing modes couple to higher modes that are not included in the Fourier series, and these coupling contributions are aliased back onto the existing modes, leading to inaccurate results. To check whether our results are affected by this numerical error, we first repeated the CASTOR3D run for the $\beta =0.21%,\u2009\eta =1.938\xd710\u22127\u2009\Omega \u2009m$ case with the *n* = 15 mode included and compared the individual growth rate of the *n* = 10 Fourier mode in this simulation to that in the simulation without the *n* = 15 mode also included. We found that the *n* = 10 growth rate changed by less than 0.003%, indicating that the CASTOR3D result is trustworthy. We then took a time step from the *n* = 0, 5, and 10 JOREK simulation of the $\beta =0.21%,\u2009\eta =1.938\xd710\u22127\u2009\Omega \u2009m$ case where the ballooning mode had already emerged and was linearly growing, and restarted it with the *n* = 15 mode also included. We found that the *n* = 10 growth rate decreased by about 5%–6% (depending on the exact time step at which it was restarted), becoming much closer to the CASTOR3D value. Thus, in JOREK simulations spectral blocking has more of an effect and it accounts for most of the discrepancy between JOREK and CASTOR3D results. Generally, in both codes, the growth rates of the n = 10 dominated mode are only weakly affected by the choice of including or excluding the n = 15 mode in the simulation.

## VIII. CONCLUSIONS

Continuing the work of previous studies,^{20,21} we implement a stellarator-capable reduced MHD model in JOREK and run several test cases based on the simple geometry of the Wendelstein 7-A stellarator. This paper presents the results, starting with Sec. II, which shows that the reduced model introduces an error into Lorentz force, but the error is negligible. The implementation is discussed in Secs. III and IV. In order to guarantee $\u2207\xb7B\u2192=0$ to machine precision, an analytical representation of the vacuum magnetic field (i.e., the curl-free component), as derived by Dommaschk,^{22} was used. This representation is compatible with arbitrary vacuum fields in a toroidal device. In order to run a simulation, the GVEC code is used to calculate an equilibrium, which is then used as an initial condition for the JOREK run. The actual Wendelstein 7-A simulations are presented in Secs. V–VII. Stable full MHD equilibria are preserved in the reduced model: the flux surfaces do not move throughout the simulation and closely match the flux surfaces calculated in GVEC, just as one would expect from the ordering argument in Sec. II. Furthermore, both tearing and ballooning modes were simulated, and the linear growth rates measured in JOREK are in decent agreement with the growth rates calculated by the CASTOR3D linear full MHD code.

Already in its current form, JOREK is capable of handling more complicated machines, such as Wendelstein 7-X and LHD. Benchmarks involving instabilities in these advanced devices are in progress, and the results will be reported in a future publication. We also plan to look for pressure-induced islands in high-*β* simulations to see how well their widths match the theoretical predictions of Cary and Kotschenreuther.^{42} Studies of scenarios relevant to ongoing experiments are also planned. Of particular interest are the current-driven sawtooth-like crashes observed in Wendelstein 7-X.^{43} Previous studies, which included both linear fully three-dimensional simulations with CASTOR3D^{44} and nonlinear simulations in a simplified cylindrical geometry with the TM1 code,^{45} have found that the corresponding Wendelstein 7-X equilibria are unstable to single and double modes, as well as resistive kink modes, and that the coupling of double-tearing modes with kink modes produces the sawtooth-like crashes. While the family of reduced MHD models used in JOREK, including the models derived in our previous work, cannot accurately reproduce kink modes at higher *β,*^{1} similar sawtooth-like crashes have also been simulated in TM1 at zero *β.*^{45} Using JOREK will allow us to nonlinearly simulate these modes in a fully three-dimensional geometry.

Another line of work that we intend to pursue is implementing more advanced models for stellarators than the one studied here. The most immediate improvement to the model used here would be to add $v\u2225$; other improvements include implementing separate temperatures for electrons and ions, adding a neutral density, and other model extensions already implemented for the tokamak models.^{1} Going further, we also intend to implement a full MHD model for stellarators. This will most likely involve extending the full MHD model described in Ref. 46, which uses the standard MHD variables ${A\u2192,v\u2192,\rho ,T}$ and does not involve any ansatzes or projections, to stellarators. This may require that the vector components of $v\u2192$ and $A\u2192$ are stored in a flux-surface-aligned coordinate system instead of the cylindrical $(R,z,\varphi )$ coordinate system as in Ref. 46. However, implementing the full MHD model with the ${\Psi ,\Omega ,\Phi ,v\u2225,\zeta ,\rho ,T}$ variables, which was derived in Ref. 21 and can be seen as a direct extension of the model used in this study, would also be interesting for comparison purposes. Finally, although we do not expect any issues, it remains to be seen if the model used in this study will hold up for $\beta \u223c3\u22125%$. Further modifications to the model may be needed if future work does not produce satisfactory results for this range of *β*.

## ACKNOWLEDGMENTS

The authors thank Alessandro Zocco and Carolin Nuehrenberg for fruitful discussions, and Michael Drevlak for providing access to his EXTENDER_P code and patiently answering questions about code use.

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (grant agreement no. 101052200—EUROfusion). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: THE JOREK NONAXISYMMETRIC GRID

The spatial discretization in JOREK is performed via two-dimensional quadrilateral finite elements in the poloidal plane and a toroidal Fourier expansion. The finite element discretization has *G*^{1} continuity, meaning that any discretized functions and their first derivatives are continuous across element boundaries, but second derivatives can jump.

In each element, an element-local coordinate system $(s,t,\varphi ),s,t\u2208[0,1]$ is set up, where $(s,t)=(0,0),(0,1),(1,0),(1,1)$ corresponds to the four vertices of the element, and $\varphi $ is the geometric toroidal angle, identical to the $\varphi $ coordinate of the cylindrical coordinate system $(R,z,\varphi )$. In general, *s* and *t* can have arbitrary orientations in the poloidal plane; however, in most configurations without an X-point, *s* is the radial coordinate and *t* is the poloidal coordinate. All quantities, including the cylindrical coordinates *R* and *z*, are expressed in terms of the element-local coordinates. Expressing *R* and *z* in terms of element-local coordinates allows one to adjust the positions of the vertices of an element, which is normally used to build a flux-surface-aligned grid. Previously, *R* and *z* could only depend on *s* and *t*, but not $\varphi $;^{1} however, this constraint was removed as part of the JOREK stellarator effort. Now, the cylindrical coordinates inside a particular element are represented as

where *i* sums over the four vertices of the element, *j* sums over the degrees of freedom at each vertex, and *n* sums over the toroidal Fourier modes, with $Nctor$ being an adjustable parameter. In addition, $Bij(s,t)$ are Bezier basis functions, and

where $Ncp$ is the periodicity of the underlying geometry. Allowing *R* and *z* to depend on $\varphi $ makes it possible to build a flux-surface-aligned grid in a stellarator configuration. The physical quantities, such as density, temperature, and *ψ*, are represented in a similar way

Note that $Ntor$ and $Nctor$ are distinct parameters; on a flux surface-aligned grid, less modes are needed to reasonably represent the physical quantities than the geometry. The Fourier basis function $Zn(\varphi )$ is defined in a similar way to $Znc(\varphi )$, with the difference that $Ncp$ is replaced by $Np$; this allows running full-torus simulations without having to add unnecessary modes to the geometry.

### APPENDIX B: SELF-ADJOINTNESS OF THE LINEARIZED REDUCED MHD OPERATOR

If Eq. (3) is linearized and the nonideal terms are dropped, then the operator of the resulting linear equation will be self-adjoint. In the ideal case, Eq. (3) can be linearized as follows:

where a subscript of 0 denotes the equilibrium value of the corresponding quantity, and perturbations do not have any decorations. Now, analogous to the standard textbook derivation, let $\Phi =\u2202\xi \u2212/\u2202t$ integrate Eqs. (B1b)–(B1d) over time and insert them back into Eq. (B1a). The following linear equation is obtained as follows:

The linear operator $L(\xi \u2212)$ is analogous to the force operator $F\u2192(\xi \u2192)$ of linear full MHD but has a different dimensionality and cannot be interpreted as the force.

To demonstrate self-adjointness, one has to show that for any two scalar functions $\xi \u2212$ and $\eta \u2212$, which satisfy $\xi \u2212=0$ and $\eta \u2212=0$ on the boundary, one has $\u222b\eta \u2212L(\xi \u2212)dV=\u222b\xi \u2212L(\eta \u2212)dV$. Since $L(\xi \u2212)$ can be obtained by applying the projection operator (6) to the negative force operator $\u2212F\u2192(\u2207\xi \u2212\xd7\u2207\chi /Bv2)$, one can use the identity $\u2207f\xb7\u2207\xd7U\u2192=\u2212\u2207\xb7(\u2207f\xd7U\u2192)$ and integration by parts to write the following:

as in Ref. 21. Here, $B\u21920=\u2207\chi +\u2207\Psi 0\xd7\u2207\chi $, which allows one to express the magnetic field and pressure perturbations in a familiar way

Now identify $\xi \u2192=\u2207\xi \u2212\xd7\u2207\chi /Bv2$ and $\eta \u2192=\u2207\eta \u2212\xd7\u2207\chi /Bv2$, and note that the boundary conditions on $\xi \u2212$ and $\eta \u2212$ imply that the vector fields $\xi \u2192$ and $\eta \u2192$ satisfy the usual boundary condition for displacement in a plasma surrounded by a wall: $\xi \u2192\xb7n\u2192=\eta \u2192\xb7n\u2192=0$. Finally, one can apply the self-adjointness property of the force operator $F\u2192(\xi \u2192)$, which allows us to write

This concludes the self-adjointness proof for the linear operator $L(\xi \u2212)$.