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

Although the basic concept of a stellarator was known since the early days of fusion research, advances in computational technology have enabled the modelling 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 modelling of stellarators in order to better understand their performance and operational limits. This paper reports on the extension of the JOREK code to 3D geometries and on the ﬁrst stellarator simulations carried out with it. The ﬁrst 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 ﬂux surfaces do not move throughout the simulation, and closely match the ﬂux surfaces of the full MHD equilibrium. Further, 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][3][4] .However, improved mathematical models and increased computational power, which became available by the late 1980s, allowed 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 modelling 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 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][6][7] , and it wasn't 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 a) Electronic mail: nikita.nikulsin@ipp.mpg.deb) See author list of Ref 1 for full list of team members.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 paper reports on a similar extension of JOREK, one of the leading nonlinear MHD codes for tokamaks 1,15,16 , to 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][18][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 papers 20,21 , and so this paper will present the results of the second part of this effort.
As discussed in our previous papers 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 χ = F 0 φ, 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 represent an arbitrary χ analytically: 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 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, φ) grid, which we do using the EXTENDER P code 23 .
Using the Dommaschk potential formulation for χ in conjunction with non-axisymmetric flux surface aligned grids allows one to simulate stellarators relatively efficiently.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 24 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 24 code mentioned above is a new fixedboundary 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 non-uniform 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 section IV explains how the initial conditions for the reduced MHD variables can be calculated from the GVEC equilibrium.In section 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, section VI shows tearing mode simulations in Wendelstein 7-A and benchmarks the growth rates against the CAS-TOR3D 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 non-axisymmetric 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 section II has a selfadjoint 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 time scale.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 corresponds mostly to field compression, but also adds a small correction to field line bending.Since the vacuum field ∇χ must be divergencefree, it can be written in terms of Clebsch potentials: 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, B v = |∇χ|.The reduced model is obtained by setting ζ = 0 and Ω = 0, as well as 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 where the ideal gas law p = ρRT applies, and the short- The operators are defined as follows: 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 the equations (3), η is the resistivity, µ ⊥ is a viscositylike parameter (it is not the same as physical dynamic viscosity, as discussed in Ref 21), µ h is the hyperviscosity, η h is the hyperresistivity (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 auxilliary variables: the normalized current in the ∇χ ) and the contravariant χ component of vorticity ω = −∇χ • ω = −∇χ • ∇ × v, both taken with the opposite sign.The definition equations (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 where f is an arbitrary scalar field and Q is an arbitrary vector field.Instead of simply substituting the definition equations (3e) and (3f) into the rest of the equations (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 evaluated at each time step simultaneously with the degrees of freedom for the other variables by using the definition equations (3e) and (3f) alongside with the rest of the equations (3).This approach, used in combination with transforming the equations to weak form, allows one to avoid second order derivatives in equations (3), except for the hyperviscosity term in equation (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 equation (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 n continuity, where n is a user-selected parameter 33 .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 (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 added separately after the projection operator is applied (see Ref 21 for more details).
The projection operator that produces equation (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 equation ( 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 equation (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 and Ω = O(λ 2 ).Meanwhile, the residual force due to the reduction is 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 µ 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 Here, g is the metric tensor of the Clebsch-type coordinate system aligned to ∇χ, which was introduced in Ref 20, q i represents the actual coordinates: q i ∈ {ψ v , β v }, and e i are the contravariant basis vectors: 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 , where e βv = J∇χ × ∇ψ v is the covariant basis vector in the β v direction in the Clebsch-type coordinate system aligned to ∇χ, and 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 (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 B 3 v [∂Ω/∂β v , Ψ], will be cancelled 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 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 cancelled by the projection operator ( 6) and as such will not contribute to the Φ equation (3a), however it will not be cancelled by the projection operator (7).Indeed, its image under the operator (7) will be 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 equation (3a) and the v evolution equation, some more ordering needs to be done.Consider that the shortest time scale in the reduced system is the Alfvén time τ A ≡ L /c A , 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, , it is clear that the lowest order terms in equations ( 5) are O(λ 2 ), and both projection operators ( 6) and ( 7) are O(1).As such, the lowest order terms in both equation (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 section 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 F 0 φ corresponds to a tokamak-like toroidal field, n is the toroidal mode number, m is the poloidal mode number, and each harmonic satisfies the Laplace equation individually: ∆χ n,m = 0. Dommaschk gives a more explicit representation for χ 22 : where a tilde denotes normalization: χ = F 0 χ, R = R 0 R and z = R 0 z; the normalization factor R 0 is the toroidally averaged radial position of the magnetic axis of the vacuum field.The functions D n,m and N n,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 a n,m , b n,m , c n,m and d n,m in equation ( 14) are what determines a particular configuration and must be calculated from the EXTENDER P output.
Note that, since the harmonics χ n,m are given analytically, the property that ∆χ n,m = 0 is satisfied exactly.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, (18) We will make use of the properties (from Ref 22, equations (10) and ( 11) (19) If one also evaluates at z = 0 and integrates over φ, F 0 can be calculated: To calculate the coefficients a n,m and b n,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 equations ( 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 equations ( 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 section 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 P i (z ), integrating from -1 to 1. Here, Z is determined as follows.In each poloidal plane at 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, Figure 2 shows a segment of the surface of integration in one field period of the Wendelstein 7-A equilibria described in section V.
The Legendre polynomials are orthogonal to monomials of 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 a n,m : where z i , P j (z) = 1 −1 z i P j (z)dz.Similarly, for the coefficients b n,m , one has the following linear algebraic system: (24) Again, evaluating at R = 1 and using the properties ∂D n,m /∂ R| R=1 = 0 and ∂N n,m /∂ R| R=1 = z m /m! (also from Ref 22, equations ( 10) and ( 11)), one has: From here, it is straightforward to follow the same steps as for a n,m and b n,m , obtaining the following linear algebraic systems for c n,m : and for d n,m : Note that there are only M equations in each system for the unknowns c n,1 , . . ., c n,M and d n,1 , . . ., d n,M because N n,−1 is not defined, and so terms with c n,0 and d n,0 are not included in the sum (14).
The only coefficients for which a system of equations has not yet been obtained are a 0,m (there are no b 0,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 : To finalize the derivation, we multiply the equation by a Legendre polynomial P i (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 a 0,1 , . . ., a 0,M .This is because D 0,0 = 1, and so a 0,0 is an additive constant in the scalar potential, which has no effect on the vacuum magnetic field 22 .
The linear algebraic systems of equations ( 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 section II, although j and Ψ are related by j = ∆ * Ψ, j is stored as a separate variable in finite element representation for numerical purposes.It makes sense to first calculate the initial condition for j from the GVEC data, and then calculate Ψ 0 from j 0 using equation (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 ∇χ, one has, after some algebra: Using the same ordering as in section II, where B v = O(1), Ψ = O(λ), Ω = O(λ 2 ) and ∂ = O(λ), it can be seen that the first term is O(λ) and 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 Having determined j 0 , it remains to solve the differential equation (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 paper, 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 equation (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 F 0 φ term in equation ( 13)), so ∂χ/∂φ is nonvanishing and the Jacobian of the new coordinates is nowhere singular.The boundary equation in (ψ, θ, χ) coordinates is where 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 + i bnd elm )/N bnd elm , where i bnd elm is the zero-based index of the current boundary element and N bnd 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 equation ( 34) is searched for can now be represented as: where N cp is defined in Appendix A and χ = χ/F 0 .Excluding the m = 0 mode removes the kernel of the differential operator of equation ( 34) from V sol , 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, equation ( 32) is solved by splitting Ψ 0 = Ψ 0,i + Ψ b , where Ψ b is the solution to equation (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 p GVEC , which must simply be 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 ρ 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 done using stable equilibria to demonstrate that the model is indeed consistent, the error due to neglecting of fourth-order terms in the Lorentz force, which was discussed in section 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 done with five-fold 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, ..., 12 and toroidal modes n = 0, 5, 10, ..., 50, which corresponds to N ctor = 21 and N cp = 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 ψ tn is the toroidal flux normalized so that ψ tn = 0 at the axis and ψ tn = 1 at the boundary.I n (ψ tn ), which represents the toroidal current enclosed by the flux surface ψ tn , is normalized by the total toroidal current, which was 17.5 kA in the cases considered, such that I n (1) = 1.The total toroidal magnetic flux through a poloidal plane was 0.08 Wb.The pressures at the axis were 1 Pa, 100 Pa, 500 Pa and 1 kPa, which corresponds to β-values of 2.3 • 10 −5 %, 2.3 • 10 −3 %, 0.011 % and 0.022 %, respectively.The pressure profiles are given by where p is the pressure in pascals, p a is the pressure at the axis and p b is the pressure at the boundary.For the β = 2.3 • 10 −5 % (p a = 1 Pa) case, p b = 0.01 Pa, while for the other three cases, p b = 1 Pa.When finding the initial conditions from the GVEC equilibrium, N tor = 9 and N p = 5 (see appendix A) was used for the variables, which corresponds to Fourier modes n = 0, 5, ..., 20.The profile of the rotational transform ι, which was the same for all equilibria considered in this paper, is shown in Figure 3.All of the simulations were run with a spatially constant resistivity η = 1.938 • 10 −9 Ω•m and viscosity µ = 2.90 • 10 −9 kg/(m • s).In addition, a hyperviscosity µ h = 2.90 • 10 −12 kg • m/s was applied.The radial resolution of the finite elements was 41 nodes, and the poloidal resolution was 48 nodes, i.e. 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 Figure 4).This consisted of 20 time steps of length 6.484 • 10 −4 ms (1 in JOREK units), followed by 20 time steps of length 6.484 • 10 −3 ms (10 in JOREK units), followed by 10 time steps of length 6.484•10 −2 ms (100 in JOREK units).For the β = 0.022% case, but not for the others, this was followed by another 10 time steps of length 6.484 • 10 −2 ms.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 (10000 in JOREK units).The for numerical stability.When evaluating the integrals in the weak form of the equations (3), the toroidal integration was done by summing over 40 poloidal planes spread evenly 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 Figure 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 φ = 0 poloidal plane, and then tracing the field line at that position for ten toroidal turns, after which the tolerance T = 0.1 (max R i − min R i )(max z i − min z i ), where i = 1, ..., 10, is calculated.If this tolerance is smaller than the cutoff, which was set to 5 • 10 −5 m, then the axis is considered found: the axis position at If not, then the field line tracing is restarted at (R c , z c ), and the process is repeated until the tolerance is less than the cutoff.The error in the R-coordinate was was estimated as E R = 0.1(max R i − min R i ) and plotted as error bars in Figure 5.
To demonstrate that there is no significant motion even away from the axis, the Poincare plots for the β = 0.022 % case are shown in Figure 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•10 −5 %, 2.3•10 −4 % and 2.3•10 −3 %.For reference, Figure 7 shows the velocity stream function Φ without the n = 0, 5, 10 Fourier modes, which do not contribute to the tearing mode, on the φ = 0 poloidal plane during the pre-saturation (linear) phase of the full torus simulation (see below).The characteristic (2,1) structure of the mode is clearly visible.The β = 2.3 • 10 −5 % and β = 2.3 • 10 −3 % are the same equilibria that were used in the previous section, with the β = 2.3 • 10 −4 % 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, p a = 10 Pa and p b = 0.1 Pa.When finding the initial conditions from the GVEC equilibrium, N tor = 5 and N p = 5 (see appendix A) was used for the variables, which corresponds to Fourier modes n = 0, 5, 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 • 10 −4 ms, followed by 20 time steps of length 6.484 • 10 −3 ms, followed by 5 time steps of length 6.484The 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 N tor = 5, anisotropic transport cannot be properly modelled, 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, ..., 10 Fourier modes, corresponding to N tor = 21 and N p = 1 (see appendix A).The Crank-Nicolson scheme was used with time steps of length   28,29 .The maximum deviation between the two codes is 13.2%, and occurs at β = 2.3 • 10 −4 %.
For reference, Figure 9 shows the magnetic energies of each individual Fourier mode in the β = 2.3•10 −5 % 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, 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, 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 Figure 10.Note that, aside from the dominant (2,1) island chain, there is a secondary (3,2)  island chain towards the interior of the plasma, which is nonlinearly excited by the mode.
Two more simulations were done with the β = 2.3 • 10 −5 % case, this time using resistivities of η = 1.938 • 10 −7 Ω•m and η = 1.938 • 10 −5 Ω•m, while all of the other parameters were kept the same as before.For the η = 1.938 • 10 −7 Ω•m case, a hyperresistivity of η h = 9.691 • 10 −14 Ω•m 2 (5 • 10 −14 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 hyperresistivity in the η = 1.938 • 10 −6 Ω•m case, which could be run with or without hyperresistivity, changes the growth rate by less than 1.5%.Figure 8 b shows the growth rates for the β = 2.3•10 −5 % case at different values of resistivity alongside the growth rates calculated by CASTOR3D.The maximum deviation between the two codes is 26.2%, occuring at η = 1.938 • 10 −5 Ω•m.This is most likely due to the neglect of v by the model used in these simulations, as v 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 quasiinterchange modes.Neglecting the parallel velocity leads to overestimation of the quasi-interchange growth rates by a factor of 1 + 2q 2 s , 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 β = 0.11% and β = 0.21% (corresponding to axis pressures of 5 kPa and 10 kPa) has been done at resistivities of η = 1.938 • 10 −7 Ω•m and η = 5.814 • 10 −7 Ω•m.For reference, the velocity stream function Φ is shown in Figure 11 during the linear phase of the β = 0.21%, η = 1.938 • 10 −7 Ω•m simulation.These equilibria have the same toroidal current and pressure profiles as the equilibria in section V, with p b 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 • 10 −4 ms.In the second part of the simulation, the Crank-Nicolson scheme was used, however both parts were run with a five-fold 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 → ∞ 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 > n max 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 N tor = 5, a convergence test was done for the β = 0.21% equilibrium while using a resistivity of η = 1.938 • 10 −7 Ω•m.The growth rate converged at a finite element resolution of 61 nodes radially and 72 nodes poloidally, time step size of 6.484 • 10 −4 ms and a hyperviscosity of µ h = 7.25 • 10 −15 kg • m/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 β = 0.11%, η = 1.938 • 10 −7 Ω • m, 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 if our results are affected by this numerical error, we first repeated the CASTOR3D run for the β = 0.21%, η = 1.938 • 10 −7 Ω • m 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.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, 10 JOREK simulation of the β = 0.21%, η = 1.938•10 −7 Ω•m case where the ballooning mode had already emerged and was growing linearly, 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 discrepancy between JOREK and CASTOR3D results.Generally, in both codes, the growth rates of the n=10 dominated mode are affected only weakly by the choice of including or excluding the n=15 mode in the simulation.

VIII. CONCLUSION
Continuing the work of previous papers 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 section II, which shows that the reduced model introduces an error into Lorentz force, but the error is negligible.The implementation is discussed in sections III and IV.In order to guarantee ∇ • B = 0 to machine precision, an analytical representation of the vacuum magnetic field (i.e. the curlfree 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 sections V, VI and 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 section II.Further, 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 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 as well as 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 tearing 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 to simulate these modes nonlinearly 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 ; 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, v, ρ, T } and doesn't involve any ansatzes or projections, to stellarators.This may require that the vector components of v and A are stored in a flux-surface-aligned coordinate system instead of the cylindrical (R, z, φ) coordinate system as in Ref 46.However, implementing the full MHD model with the {Ψ, Ω, Φ, v , ζ, ρ, T } variables, which was derived in Ref 21 and can be seen as a direct extension of the model used in this paper, 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 paper will hold up for β ∼ 3 − 5%.Further modifications to the model may be needed if future work does not produce satisfactory results for this range of β.

FIG. 1 .
FIG. 1.The volume-averaged squared relative error of the Dommaschk potential representation (∇χ − BE) 2 /B 2 E as a function of the number of poloidal modes M .The values shown in this plot were calculated using a Python implementation of Dommaschk potentials based on the one written by Paul Huslage for the BOUT++ code 34,35 .
) As can be seen, both of these systems of equations have triangular matrices.At this point, the equations for the coefficients c n,m and d n,m have yet to be determined.Consider now the R component of B E : B E,R = R • B E = ∂χ/∂R.One has: ,m cos nφ + b n,m sin nφ) ∂D n,m ∂ R + (c n,m cos nφ + d n,m sin nφ) ∂N n,m−1 ∂ R .

FIG. 3 .
FIG.3.The ι profile as a function of ψtn.This profile is the same for all configurations considered in this paper.

FIG. 4 .
FIG. 4. The total kinetic energy of the plasma in the β = 0.022 % case during the first 0.144 ms of the simulation (20 time steps of 1 and 20 time steps of 10 JOREK time units) showing the damping out of motion due to the neglect of fourth order terms in the Lorentz force.

FIG. 6 .
FIG.6.The Poincare plots (black points) for the β = 0.022 % case at t = 0 and φ = 0 (a), at φ = 0 after the simulation is over (t = 7.275 ms) (b), and at t = 0 and φ = 3π/10 (3/4 of the way through one period) (c).The flux surfaces in the full MHD equilibrium as calculated in GVEC are shown by the orange lines.
1.621•10 −2 ms.The toroidal integration in the weak form of equations (3) was done by summing over 40 poloidal planes spread evenly throughout the full torus.The number of Fourier modes, 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, number of modes or number of planes leads to a change in the growth rate of less than 1.5%.The convergence test was done for the β = 2.3 • 10 −3 % and β = 2.3 • 10 −5 % cases, resulting in all of the parameters converging to the same values, except for hyperviscosity, which converged to µ h = 7.25 • 10 −15 kg • m/s for the β = 2.3 • 10 −3 % case and µ h = 2.90 • 10 −15 kg • m/s for the β = 2.3 • 10 −5 % case.The β = 2.3 • 10 −4 % 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
.5.The R coordinate of the magnetic axis as a function of time for the four different β cases.The axis position in the full MHD equilibrium as calculated in GVEC is shown by the black dot at t = 0. FIG