Testing of the new JOREK stellarator-capable model in the tokamak limit

In preparation for extending the JOREK nonlinear MHD code to stellarators, a hierarchy of stellarator-capable reduced and full MHD models has been derived and tested. The derivation was presented at the EFTC 2019 conference. Continuing this line of work, we have implemented the reduced MHD model (arXiv:1907.12486) as well as an alternative model which was newly derived using a different set of projection operators for obtaining the scalar momentum equations from the full MHD vector momentum equation. With the new operators, the reduced model matches the standard JOREK reduced models for tokamaks in the tokamak limit and conserves energy exactly, while momentum conservation is less accurate than in the original model whenever field-aligned flow is present.


Introduction
Reduced MHD, a system of approximations introduced in its original form in the 1960s (Greene & Johnson 1961), continues to be used by modern MHD codes, such as JOREK (Franck et al. 2015) and M3D-C 1 (Breslau et al. 2009). While many versions of reduced MHD consisting of different equations and using different methods to derive them have been published, the main idea is the same: the removal of fast magnetosonic waves while retaining relevant physics (Strauss 1997;Jardin et al. 2012). The removal of these waves eliminates the shortest time scale and allows one to use larger time steps due to the Courant condition. Even when implicit time integration methods are used, and the Courant condition is no longer a hard limit, using time steps that are large compared to the shortest time scale can lead to poor accuracy (Jardin et al. 2012;Kruger et al. 1998). An increased time step is especially advantageous in nonlinear simulations, which otherwise tend to be computationally demanding. In particular, stellarator simulations, which are even more computationally expensive than tokamak simulations due to the complex geometry, can benefit from an increased step size.
Originally, reduced MHD was derived via an ordering in a small parameter, often taken to be the inverse aspect ratio. The ordering itself is a system of several approximations and assumptions involving the ordering parameter that allows one to determine the relative order (in terms of the ordering parameter) of any quantity with respect to any other quantity of the same dimension. In this context, terms corresponding to fast magnetosonic waves have a higher order than the terms that one wants to keep, allowing them to be dropped. Different versions of reduced MHD can be derived by using different orderings

Changes made to the model
The original hierarchy of models, as presented in Nikulsin et al. (2019), was derived from the viscoresistive MHD equations: (2.1) The gradient operators parallel and perpendicular to the total magnetic field B are defined as ∇ = B B 2 B · ∇ and ∇ ⊥ = ∇ − ∇ . The usual MHD notation is followed with ρ, p, v and B being density, pressure, velocity and magnetic field, respectively. In addition to that, η is the resistivity, µ is the dynamic viscosity, D ⊥ is the mass diffusion coefficient perpendicular to field lines, κ ⊥ and κ are the thermal conductivity across and along field lines, and S ρ and S e are source terms in the continuity and energy equations, respectively. The ideal gas law p = ρRT is assumed to hold.
We also introduced the following ansatzes for the magnetic field and velocity: and proved that any magnetic field and any velocity can be represented in such a way.
Here, ∇χ is the curl-free vacuum component of the magnetic field, which is generated by the coils, B v = |∇χ| and ∇ ⊥ = ∇ − B −2 v ∇χ(∇χ · ∇). Note that, since the ansatzes do not restrict the magnetic and velocity fields, they can still be used within the context of full MHD. The transition to reduced MHD is done by setting Ω = ζ = 0 and dropping the evolution equations for Ω and ζ (derived below). Further, in the tokamak limit, i.e. when χ = F 0 φ, where φ is the toroidal angle † , the ansatzes reduce to those of the standard JOREK reduced MHD tokamak models, B = F 0 ∇φ+∇ψ×∇φ and v = R 2 ∇u×∇φ+v B, where ψ = F 0 Ψ and u = Φ/F 0 (Hoelzl et al. 2021).
The scalar ψ v is one of the Clebsch potentials for the vacuum field, which can locally be written as (2.3) As we discussed in Nikulsin et al. (2019), one can construct a Clebsch-type coordinate system with the coordinates (ψ v , β v , χ). Note that the variables ψ v and β v are not unique, as can be seen by replacing ψ v with ψ v = ψ v + f (β v ), where f is an arbitrary function, which leaves equation (2.3) unchanged. In general, given some vacuum magnetic field ∇χ, the scalars ψ v and β v are completely determined by fixing some parameterization (ψ v , β v ) of a surface intersected by the vacuum field lines, with the values of ψ v and β v in the rest of space being determined by following the field lines off of the surface while keeping ψ v and β v constant. In a toroidal system, the surface is usually chosen to be a poloidal plane, and, in many cases, a cut needs to be introduced at that poloidal plane to prevent multivaluedness of ψ v and β v when the field lines travel once around the torus and encounter previously labelled points. Finally, β v is determined by integrating 'haeseleer et al. 2012;Stern 1970). Usually, it is advantageous to find a ψ v such that the component of ∇Ω ×∇ψ v perpendicular to ∇χ, F 2 v |∂ Ω| 2 , is minimized at t = 0, so that most of field line bending is captured by the ∇Ψ × ∇χ term.
The variables Ψ , Ω, Φ, v and ζ are the unknown magnetic field and velocity variables, which are solved for. In order to obtain evolution equations for these variables, we projected Faraday's law onto ∇ψ v and ∇χ, and applied the following projection operators † In many tokamaks, the vacuum field also includes a contribution from the poloidal field coils, however the toroidal field is by far the dominant component of the vacuum field, and it makes sense to include only the toroidal field in ∇χ. Throughout this paper, we will refer to χ = F0φ as the tokamak limit.
to the vector momentum equation, which was first divided by ρ Here, e χ = B/B χ is the covariant basis vector in the Clebsch-type coordinate system (α, β, χ) associated with the full magnetic field, B = ∇α×∇β. These projection operators have the property that each one of them is orthogonal to all but one term in the velocity ansatz (2.2), i.e.
However, for the sake of energy conservation, a better set of projection operators would be as we will show in the next section.
To obtain the new evolution equations for the velocity variables, we insert the ansatzes (2.2) into the momentum equation and apply the projection operators (2.5), this time without dividing first by ρ. The momentum equation can be written as where we used the identity ( v · ∇) v = 1 2 ∇v 2 + ω × v and ω is the vorticity: (2.7) Just as in Nikulsin et al. (2019), we do not carry the viscosity term through this derivation, instead adding a generic viscosity term afterwards. To obtain the evolution equation for Φ, we apply to equation (2.6) the first projection operator in (2.5), or its equivalent, −∇ · (B −2 v ∇χ×, when appropriate: We have added the generic viscosity term ∇ · (µ ⊥ ∇ ⊥ ∆ ⊥ Φ), which we have chosen to match the viscosity term in Franck et al. (2015) in the tokamak limit. As discussed in previous papers (Franck et al. 2015;Nikulsin et al. 2019), the viscosity term µ∆ v in the momentum equation (2.1) is only a rough approximation for the divergence of the viscous stress tensor in a plasma, and one does not loose much by using a generic viscosity term in the final equations. Note, however, that the µ ⊥ in the above equation is not the physical viscosity µ in equation (2.1). Indeed, from dimensional analysis it follows that µ ⊥ has units of T 2 · Pa · s. Applying scaling analysis to the term µ∆ v after inserting just the first term of the ansatz (2.2) for v and applying the first projection operator in (2.5), we see that where L ⊥ is the length scale perpendicular to the magnetic field, and L ⊥ L for most magnetic fusion configurations. Applying the same analysis to the generic viscosity term, we get Comparing the two scalings, we see that µ ⊥ scales as µB −2 v ; for the purposes of the scaling, one can take the value of B v at the axis as a typical value and write µ ⊥ ∼ µB −2 v,axis . To get the evolution equation for v , we apply the second operator in (2.5), projecting equation (2.6) on B: (2.9) where we have again added a generic viscosity term with a form analogous to the viscosity term in the previous equation. Here, (f, g) = ∇ ⊥ f · ∇ ⊥ g is an inner product, is a Poisson bracket with respect to ∇ψ v . Finally, to obtain the evolution equation for ζ, we apply the last projection operator in (2.5) (or its equivalent, −∇ · [B −2 v ∇χ × (∇χ×, to some terms) to equation (2.6): (2.10) This concludes the derivation of the scalar momentum equations.
The evolution equations for the magnetic field variables of Nikulsin et al. (2019) are valid, therefore we keep them as is: (2.11) The continuity equation remains unchanged (Nikulsin et al. 2019): while for the pressure, we use the standard internal energy evolution equation, instead of the full energy conservation equation from the system (2.1): (2.13) which now, after using (2.2), reads (2.14) We note that since any arbitrary magnetic field and velocity can be represented in the ansatz form (2.2), the equations derived above still correspond to full MHD, albeit in a unconventional form. We also showed that if the background vacuum field is the dominant part of the magnetic field, then the MHD waves are approximately separated by the terms of the velocity ansatz, with the first term containing Alfvén waves, the second containing slow magnetosonic waves, and the last one containing fast magnetosonic waves (Nikulsin et al. 2019). Thus, to get the reduced MHD system, we set ζ = Ω = 0, eliminating fast magnetosonic waves and field compression, and drop the corresponding evolution equations, namely equation (2.10) and the second equation in (2.11), to avoid having an overconstrained system. We also approximate the electric field in Ohm's Law as E = − v× B +η j , which corresponds to neglecting the components of j perpendicular to ∇χ in the last term of the first equation in (2.11). Alternatively, if we introduce a tensor resistivity, ← → η , this approximation corresponds to neglecting the perpendicular resistivity, η ⊥ ≈ 0, although parallel and perpendicular resistivities are usually of the same order, so it makes more sense to speak about neglecting the perpendicular current (only in the resistive term). This seems to be a fairly good approximation, at least for tokamaks (Pfirsch 1978), and likely also for any configuration with significant bootstrap current. Indeed, we can write ← → η · j = η ⊥ j ⊥ + η j , where subscripts refer to parallel and perpendicular components relative to the total field B, not just ∇χ. When the bootstrap current is high enough, the η j term is dominant, and, since ∇χ is the dominant component of B, we have j ≈ j in the lowest order. As we will show in the next section, the approximation is necessary in order to have energy conservation in reduced MHD. Accordingly, we also have to replace the last term of equation (2.14) with η(j ) 2 to avoid creating extra internal energy, as now only the parallel component of current contributes to the dissipation of magnetic energy. A further reduction would be to also eliminate slow magnetosonic waves by setting v = 0 and dropping equation (2.9).
Finally, we point out that in the tokamak limit, when χ = F 0 φ, the new reduced models derived in this section match the currently implemented JOREK tokamak models (Hoelzl et al. 2021). The Ψ evolution equation, once we set ζ = Ω = 0 and drop the perpendicular components of the current, becomes which can be rewritten as The above equation will be satisfied as long as is satisfied, which exactly matches the currently implemented form of the magnetic field equation in JOREK when χ = F 0 φ (Hoelzl et al. 2021). For all other equations, the fact that they match the corresponding JOREK reduced MHD equations in the tokamak limit is obvious and can be verified by simple substitution of χ = F 0 φ.

Energy conservation
We begin by showing that the new models conserve energy. If we apply the first projection operator in (2.5) to some vector field Q, multiply the result by a test function Φ * and integrate over the plasma volume, we can get the following expression using the identity ∇f · ∇ × U = −∇ · (∇f × U ) and integration by parts: where we let Φ * = 0 on ∂V , so the surface integral term is zero. Doing the same with the third projection operator and test function ζ * , we have (3.2) We can also apply the second projection operator (2.5) to Q, multiply by a test function v * , and integrate over the plasma volume: , as well as Φ * = Φ, v * = v and ζ * = −ζ, and sum up the above three equations, we get where the LHS is zero due to equations (2.8), (2.9) and (2.10). More importantly, if we set ζ = 0 and drop equation (2.10), equation (3.4) still holds, since with ζ * = −ζ = 0 equation (3.2) becomes 0 = 0, and equation (2.10) no longer has to be satisfied in order for the LHS of equation (3.4) to remain zero. Similarly, if we set v = 0 and drop equation (2.9), equation (3.4) continues to hold due to equation (3.3) becoming 0 = 0 and no longer requiring equation (2.9) to be satisfied. This property of equation (3.4) being satisfied even when some terms in the velocity ansatz are missing was not present in the original models, leading to problems which we will describe in detail further in this section. We now proceed to the magnetic field evolution equations. In the full MHD case, both of the equations (2.11) are satisfied. These two equations are simply the vector components of Faraday's Law in the ∇ψ v and ∇χ directions, and as such are equivalent to the original Faraday's Law in vector form. Only two equations are needed to evolve the full magnetic field, as one degree of freedom is eliminated by the ∇ · B = 0 constraint. In the reduced MHD case, Faraday's Law is also satisfied, as one can easily see from multiplying equation (2.15) by ∇χ/B v and taking the curl. Thus, we can take the dot product of Faraday's Law with B and follow the usual procedure to obtain the magnetic energy equation: where we have j * = j in the full MHD case and j * = j in the reduced MHD case. To complete the derivation of the energy conservation law, we rewrite equation (3.4) as where we used the continuity equation. We then divide equation (2.13) by γ − 1 and intgrate it, together with equation (3.5) over the plasma volume. Adding both of the resulting equations to equation (3.6), we obtain: which can be rewritten as Here, we have assumed that the plasma is surrounded by a perfectly conducting wall, which implies the boundary conditions v · n = 0 and B · n = 0, where n is a unit normal vector to the boundary. These conditions, along with Φ = ζ = 0 on ∂V , which we had assumed before to make the boundary integral terms in equations (3.1) and (3.2) zero, are consistent with the boundary conditions that were used in the simulations presented in the following sections.
It is important to point out that, if we replace P by ∂ρ/∂t + ∇ · (ρ v), exact energy conservation continues to hold even when the exact solutions are replaced by Bubnov-Galerkin finite element approximations. Indeed, a Bubnov-Galerkin solution will ensure that the expressions (3.1), (3.2) and (3.3) equal zero whenever the test functions Φ * , ζ * and v * are finite element basis functions. Since the finite element solutions Φ, ζ and v are linear combinations of the basis functions, equation (3.6) will continue to hold. Similarly, since the Ψ and Ω solutions will be linear combinations of the basis functions and the Bubnov-Galerkin method ensures that the weak forms of equations (2.11) is satisfied whenever the test function is a basis function, the integral of equation (3.5) over volume will continue to hold. Finally, since the constant 1 can be represented as a linear combination of basis functions, the integral of equation (2.13) over volume will continue to hold, and so energy will be exactly conserved. Note, however, that this argument does not take into account temporal discretization and Fourier expansion in the toroidal direction, both of which contribute to energy conservation error in practice. However, these errors can be made small, as we discuss in section 4. The practical accuracy of energy conservation in the standard JOREK tokamak model, which is what the new model presented in section 2 reduces to in the tokamak limit, is also considered in Hoelzl et al. (2021).
The necessity of neglecting perpendicular current in the resistive term of the first equation (2.11) when using the reduced MHD model can be clarified by deriving a magnetic field equation for the case when it is not neglected. Let Ψ 1 and Ψ 2 satisfy the equations then Ψ = Ψ 1 + Ψ 2 satisfies the original first equation in (2.11), without the neglect of perpendicular current. The first of the two equations above is equivalent to The evolution equation for Ψ 2 can be written as (∂ B 2 /∂t) ψv = −(∇ × E 2 ) ψv , and so the last two terms in the RHS of the above equation will only cancel partially, leaving behind the nonconservative terms B βv ∂B βv 2 /∂t + B βv (∇ × E 2 ) βv + (∇ × E 2 ) χ . Finally, we show why the projection operators used in our original model allowed for nonphysical generation of kinetic energy. To begin with, the velocity field can be written in a covariant form: as opposed to the contravariant form (2.2). The fact that any vector field can be written in this way is easy to show using the same arguments as in Nikulsin et al. (2019), except using the projection operators (2.5) instead of the operators (2.4). Note that the projection operators (2.5) are orthogonal to all but one term in the covariant velocity ansatz (3.8): Now, applying the first projection operator in (2.4) to a vector field Q and then multiplying by a test function Φ * and integrating, we get: where we again used integration by parts to get the RHS. Similarly, for the third projection operator in (2.4), we can write: and for the second projection operator, we simply multiply by v * and integrate, without doing any transformations: The first three conditions pose no problem in the full MHD case, but if we try to set ζ = 0, it does not imply ζ = 0, and thus when we drop the evolution equation for ζ, equation (3.10) is no longer satisfied, and so the kinetic energy equation (3.6) cannot be satisfied. In addition, because the vector momentum equation was divided by density before the projection operators were applied, unless the density is held constant, the kinetic energy equation (3.6) also cannot be satisfied even when ζ = 0. Thus, both the reduction and spatially varying density can lead to nonphysical generation of kinetic energy.

Numerical examples
In this section, we consider several simulations with the JOREK code for a test case of a tearing mode in a circular cross section tokamak with an aspect ratio of 10 ( Figure 1). We compare the reduced model without parallel flow from the original set of stellaratorcapable models, as presented in Nikulsin et al. (2019), to the standard JOREK tokamak model without parallel flow (Hoelzl et al. 2021). Note that in the tokamak limit χ = F 0 φ, the latter is equivalent to the the model without parallel flow from the set of new stellarator-capable models, as introduced in section 2. In addition, we point out that the standard tokamak reduced MHD model in JOREK is known to accurately reproduce tearing modes when compared to full MHD (Pamela et al. 2020;Haverkort et al. 2016), thus comparing to the standard tokamak model is sufficient.
In the test cases considered below, F 0 = 10 T · m. The initial conditions were set up by solving the Grad-Shafranov equation with F F (ψ n ) = 1.173(1 − ψ n ) in units of T 2 · m, where ψ n = (ψ − ψ axis )/(ψ edge − ψ axis ), and p(ψ n ) = 0, except for the last case. In the tokamak limit, the ψ in the Grad-Shafranov equation is related to the Ψ in the magnetic field ansatz (2.2) by ψ = F 0 Ψ . The density was initialized to be constant at 3.346 · 10 −7 kg/m 3 , corresponding to 10 20 deuterium ions per cubic meter, and Φ was initialized to zero. A viscosity of µ = 5.159 · 10 −8 kg/(m · s) was used in all simulations, which corresponds to 10 −5 in JOREK units for µ t ⊥ in the standard tokamak model. † When using the original stellarator model, the kinematic viscosity was set to be ν = µ/ρ 0 = 0.1542 m 2 /s (10 −7 in JOREK units). Finally, since JOREK discretizes the toroidal direction by approximating all solutions as toroidal Fourier series, in all simulations shown in this section, the Fourier series are truncated after the first mode for the sake of simplicity, keeping only the axisymmetric (n = 0) and the n = 1 terms, unless noted otherwise. Finite elements are used for discretization in the poloidal plane. Temporal discretization is done via the Crank-Nicolson scheme, i.e. an average of the forward and backward Euler schemes. Given equations of the type where U is the set of all unknowns in the model and a i , A i and B are the expressions that comprise the equations and can involve spatial derivatives. The Crank-Nicolson scheme gives: where n in the superscript refers to the already calculated values of the model variables in the current time step, and n + 1 refers to the next time step, yet to be calculated. In addition, a temporal linearization is done around the current time step, leading to the actual numerical scheme that is implemented in the code: The viscosity µ t ⊥ in the standard tokamak model is not the same as the µ ⊥ of equation Original stellarator model Standard tokamak model Table 1. The zero-β forms of the two reduced MHD models that are compared in this section, namely the original stellarator model from Nikulsin et al. (2019) and the standard tokamak model without field-aligned flow from Franck et al. (2015). Note that the new set of equations presented in section 2 reduces to the standard tokamak model when ζ = v = Ω = 0 and χ = F0φ, except for the term containing P (eighth term in the RHS of equation (2.8)). In the cases considered, P = ∇ · (D ⊥ ∇ ⊥ ρ). When using the original stellarator model, we set χ = F0φ and ψv = R, where R is the distance from the central axis of symmetry; in addition a subscript χ means a dot product of the corresponding vector with eχ = B/B 2 v , and f b = −F F ∇Ψ |t=0/R 2 is the force balancing term (see Appendix A). In order for the initial condition to be a true equilibrium, we have introduced j d = j − j0, where j0 is the current at t = 0.
To begin with, we tested the original stellarator model (Nikulsin et al. 2019) in the linear regime by calculating the tearing mode growth rates at different resistivities and comparing them to the growth rates obtained from the standard tokamak model (Hoelzl et al. 2021). In both cases, we did a spatial and temporal resolution scan at each resistivity to ensure that the growth rate that we recorded was converged. We then repeated the simulations with Fourier modes up to and including n = 4 with the same spatial and temporal resolutions at which the growth rates converged for the n = 0, 1 simulations. The results of the n = 0, 1 simulations are shown in Figure 2 a. Figure 2 b shows the growth rates of the nonlinearly driven n = 2 modes in the simulations with n = 0, ..., 4; these modes are not inherently unstable, but are destabilized by the unstable n = 1 mode. As can be seen, both models show decent agreement. At each resistivity, the n = 2 growth rates are roughly twice the n = 1 growth rates, as expected. In both models at higher resistivities, the n = 3 and n = 4 modes are destabilized shortly before the onset of the nonlinear regime, and so the corresponding growth rates (not shown here) do not plateau, but rather peak at values roughly three and four times the value of the n = 1 growth rate and then decline rapidly as nonlinear saturation is reached. All of the growth rates shown in Figure 2 were calculated for β = 0, and the pressure was not evolved. In the case of the the original stellarator model, this amounts to not using the full energy conservation equation (third equation in (2.1)). For clarity, we explicitly list the equations solved in Table 1.
While the original stellarator model (Nikulsin et al. 2019) performs decently in the linear regime and at β = 0, two problems can arise after nonlinear saturation. First, as discussed in the previous section, nonphysical kinetic energy can be generated. The rate at which it is generated is negligible in the linear regime, but can become significant as saturation is reached. This can be seen in Figure 3, where we consider the tearing mode   Table 1.
simulation with a resistivity of 10 −5 JOREK units (1.9382 · 10 −5 Ω · m). In Figure 3, −dE/dt is plotted and compared to the physical energy loss rate. We define E = V EdV , where E is given in (2.1); this is the integrated total energy of the system at a given time. The physical energy loss rate is defined as the sum of the energy fluxes across the plasma boundary and the volume integral of energy sinks due to resistive and viscous dissipation (conversion to internal energy is not accounted for since pressure is not evolved, so dissipated energy is lost). If the inward energy fluxes are greater than outward fluxes plus sinks, the physical loss rate is negative, i.e. energy is gained. As can be seen in Figure  3 a, the difference between −dE/dt and the physical loss rate is negligible in the linear  regime (until approximately 20 ms) and grows rapidly with the onset of saturation, due to the uncontrollable increase in kinetic energy (Figure 4 a). The simulation is stopped at ∼27 ms, as it would crash shortly thereafter if allowed to continue. In order to allow the simulation to continue after saturation, we can introduce an artificial dissipation. Figure  3 b shows −dE/dt and physical energy loss rate for a simulation with artificial dissipation in the form of a hyperviscosity term (with ν h = 1.08·10 −3 m 4 /s, corresponding to 7·10 −10 in JOREK units) in the evolution equation for Φ. In this case, the physical energy loss rate does not include the hyperviscous dissipation. As can be seen, with artificial dissipation, the energy behaves much more reasonably, though it is still not conserved, hence the mismatch between the total rate of energy change and physical energy loss. The reason for such good behavior is that the non-conservation comes mostly from the kinetic energy, which is kept from exploding by the hyperviscosity term. In addition, the hyperviscosity prevents the plasma from responding too quickly to numerical errors in the magnetic field, thus stabilizing it against numerical instabilities that arise from the Ψ evolution equation (discussed below). One can introduce a finite pressure into the original stellarator model by using equation (2.14) without affecting the energy conservation by much, since all of the error comes from the kinetic energy, which is usually not dominant in fusion-relevant situations. However, we find that the original stellarator model cannot correctly reproduce the growth rates when β becomes non-negligible, likely because the pressure terms in the original stellarator model (Nikulsin et al. (2019), not shown in Table 1) contain either ∇ρ or ∂ p, both of which are zero initially, and so they make the pressure terms small compared to the Lorentz force terms.
Finally, for the sake of comparison, we show the same kind of graph for a simulation with the standard tokamak model/new stellarator model without parallel flow (Figure  3 c), both of which are equivalent in the tokamak limit. Note the small discrepancy between −dE/dt and the physical energy loss rate in the last plot; unlike the discrepancy in Figure 3 b, this discrepancy is purely numerical. As the arguments of section 3 continue to apply even after the exact solutions are replaced by finite element approximations, the discrepancy is not due to poloidal discretization, and is independent of the resolution of the finite elements, which we have confirmed by running simulations with higher resolution. Instead, the two sources of error are: a) the temporal discretization, where a small enough time step is used in practice to diminish the deviation, and b) the truncation of the toroidal Fourier series, where a sufficiently large toroidal resolution is used in practice to avoid this effect. We have carried out further simulations (not shown here), where we have confirmed that by decreasing the time step size and increasing the number of Fourier modes kept, one can make the discrepancy arbitrarily small. The second problem is that, while looking benign, the Ψ evolution equation can often produce ill-conditioned time stepping matrices, resulting in numerical instabilities. In fact, the Ψ equation is what destabilizes the kinetic energy, causing it to explode in Figure 4 a. If we replace the Ψ evolution equation in the original stellarator model (see Table 1) by equation (2.15), we can run the simulations well into the post-saturation regime without needing to introduce hyperviscosity, as can be seen in Figure 3 d, albeit with a higher energy conservation error than when hyperviscosity is present. On the other hand, we can replace the ψ evolution equation in the standard tokamak model by the corresponding equation from the original stellarator model (Table 1), thus testing the Ψ equation from the original stellarator model in a setting where energy is conserved in the continuous limit. In this case (Figure 5 a), energy conservation issues begin around the same time that the original model without hyperviscosity would have crashed (Figure  3 a). However, instead of exploding, the energy rapidly decreases, and so this simulation can run a bit longer before crashing. Figure 5 b shows that the kinetic energy plays no part in this numerical instability; all of the energy conservation error comes from the buildup of numerical errors in the magnetic energy. As can be seen, while the Ψ evolution equation in the original stellarator model in Table 1 is equivalent to equation (2.15) in the analytical sense, numerically it behaves very differently. Finally, we verified that replacing the psi evolution equation in all models considered in this section does not have any significant impact on the growth rates.
To conclude this section, we note that we had initially intended to use the full energy conservation equation, as shown in (2.1), to evolve pressure and ensure energy conservation. However, using the full energy conservation equation would make the internal energy into a reservoir, drawing energy from it when there is a nonphysical gain of kinetic energy, or depositing energy into it if there is a nonphysical loss of kinetic energy. This shifts the energy conservation error to pressure while ensuring total energy conservation. While this may be acceptable when the internal energy is much larger than the kinetic energy, using the full energy conservation equation would thus disallow low β simulations. It is much better to use equation (2.14) with the original stellarator model instead of the full energy conservation equation, although, as noted above, the original stellarator model is unable to accurately reproduce growth rates when β is not negligible.

Approximate conservation of momentum
Reduced MHD models usually do not conserve linear momentum exactly, except for some special cases. Locally, this can be seen from the fact that there are only two velocity variables, Φ and v , in reduced MHD, so the three components of the full MHD momentum equation cannot be satisfied simultaneously. However, different reduced MHD models can have different amounts of error in the momentum conservation. In the remainder of this section, we will first compare momentum conservation properties of the original model presented in Nikulsin et al. (2019) to those of the new model derived in section 2, and then show how global momentum behaves in the new model using some numerical examples. When discussing momentum conservation, we do not exclude exchange of momentum with the walls, such as in the case of a vertical displacement event. Exchange of momentum is a physical process, and, while the total momentum of the simulated system can change, no momentum is created or destroyed. In this section, we are concerned with the nonphysical generation of momentum within the plasma due to approximations made in reduced MHD.

General local momentum conservation properties of the reduced models
We begin by reviewing the momentum conservation properties of the original reduced model, which were discussed in detail in Nikulsin et al. (2019). The action of the first projection operator (2.4) on a vector Q can be written as ∇χ · ∇ × [∇χ × ( e χ × Q)] = ∇χ · ∇ × ( e χ Q χ − Q). Thus, if Q is the full MHD momentum equation (2.1), the action of the projection operator corresponds to dropping the two vector components of a vorticitylike equation, ∇×( e χ Q χ − Q), that are perpendicular to ∇χ. If all three components of this equation were satisfied simultaneously (which, in general, is not possible as noted before), then the original full MHD momentum equation would also be satisfied and momentum would be conserved exactly. We can estimate the magnitude of momentum conservation error by considering the components of the vorticity-like equation perpendicular to ∇χ. This equation can be written as (5.1) where the reduced velocity v = −∇χ × ( e χ × v) = ∇Φ × ∇χ/B 2 v and the reduced vorticity ω = ∇ × v were introduced, and the viscosity term is not considered. If the components of this equation perpendicular to ∇χ are identically zero, then there is no approximation in the reduction as nothing is being neglected, and momentum is still conserved. The most general case in which these components are zero is the following: Φ, Ψ, v , p, ρ, P } (5.2) where g ik are the components of the metric tensor of the vacuum field-aligned coordinate system. As can be shown by a simple calculation, in this case both ω and j will be directed strictly along ∇χ. If we allow either g ik , Φ or Ψ to vary along ∇χ, the same calculation will show that ω has nonzero perpendicular components. This will cause ∂ ω/∂t to have nonzero components perpendicular to ∇χ, which cannot be canceled by any other terms since there are no more time derivatives involving Φ in the equation. Similarly, if we let any of the other quantities vary along ∇χ, the last term and the seventh term on the RHS (pressure), the first term and the seventh term on the RHS (density), the last term on the LHS (P ) and third term on the LHS (v ) will have nonzero perpendicular components, which will not be canceled by any other terms. If the conditions (5.2) are met, then only the sixth and eighth terms on the LHS have nonzero perpendicular components. As can be shown by a simple expansion of the eighth term: Now we consider the new reduced model. For the first projection operator (2.5), the vorticity-like equation analogous to equation (5.1) can be written as: (5.3) Equation (2.8) is just the contravariant χ component of this equation. In the above equation, the third and eleventh terms on the LHS will have perpendicular components even when condition (5.2) is met; namely, ∇ ⊥ (ρB −2 v ∂v /∂t) × ∇χ from the third term and ∇ ⊥ (P v B −2 v )×∇χ from the eleventh term. In addition, while the eighth term on the LHS of equation (5.3) is analogous to the sixth term on the LHS of equation (5.1), there is no term in equation (5.3) analogous to the eighth term of equation (5.1) to cancel the nonzero perpendicular components of its eighth term. Thus, a more stringent condition is required to ensure that the perpendicular components of equation (5.3) are identically zero and momentum is exactly conserved in the new reduced MHD model. Namely, in addition to the conditions (5.2), we must also have v = 0. This additional condition causes the third and eleventh terms to vanish, and the eighth term becomes: Clearly, if condition (5.2) is met, the eighth term, as shown above, will not have any perpendicular components. As can be seen, the original reduced model (Nikulsin et al. 2019) has better momentum conservation properties than the new reduced model (section 2). This was the reason why we initially attempted to work with the original model. However, as energy conservation is generally more important than momentum conservation, we will use the new model for future simulations.
To quantify the momentum conservation error, we calculate the total linear momentum in the Cartesian x-and y-directions for the entire plasma. Physically, this momentum should be zero throughout the simulation, however, due to the inaccuracies discussed in the previous subsection, a nonzero momentum tends to appear. Here, we choose the x-direction as the R-direction when φ = 0 and the y-direction as the R-direction when φ = π/2, with the z-axis being the axis of symmetry of the torus. We do not consider momentum in the z-direction since the global z-momentum is conserved in the tokamak limit, as can be seen by setting Φ * = ln R in equation (3.1) and recalling that B v = F 0 /R in the tokamak case.
As can be seen in Figure 7, after about 1.5 ms, the momentum conservation error in the model with v = 0 is worse by more than an order of magnitude than that in the model with v = 0. In both cases, the instability has reached nonlinear saturation, and at that point the momenta no longer grow in amplitude, but tend to oscillate around zero with the amplitude of the oscillations remaining at the same order of magnitude.

Conclusion
In continuation of our previous work (Nikulsin et al. 2019), we have implemented the stellarator-capable reduced MHD model derived in our previous paper and tested it in the tokamak limit. We find that, even in the simple case of a tearing mode in a circular high aspect ratio tokamak, energy conservation can only be approximately satisfied, due to the possibility of nonphysical gain or loss of kinetic energy. Our original plan, as presented in Nikulsin et al. (2019), was to use the full energy conservation equation when evolving pressure, thus ensuring energy conservation. However, that would require the internal energy to be much larger than the kinetic energy, thus excluding low β situations. Even without exact energy conservation, meaningful linear results, such as growth rates, could be obtained without any further modifications, however the lack of energy conservation required artificial dissipation to be introduced to prevent a crash after nonlinear saturation was reached due to nonphysical kinetic energy buildup. Introducing artificial dissipation in the form of a hyperviscosity term kept the kinetic energies in check and allowed the simulations to proceed into the post-saturation regime. Although exact energy conservation could not be achieved, the error in energy conservation remained small because the kinetic energy is small compared to the total energy, and the magnetic energy evolved correctly. Thus, energy should be approximately conserved whenever the kinetic energy is small, which is usually the case in fusion-relevant situations. We also found that the kinetic energy blowup was triggered by numerical errors from the magnetic field. Replacing the magnetic field equation with a different form that is analytically equivalent to the original, but more numerically stable alleviates the need for hyperviscosity, although the energy conservation error in this case is larger than when hyperviscosity is present.
To remedy the lack of energy conservation in our original model, we used a new set of projection operators to obtain scalar momentum equations from the full MHD vector momentum equation. With these new projection operators, the reduced MHD model matches the standard JOREK reduced MHD models for tokamaks (Franck et al. 2015;Hoelzl et al. 2021) in the tokamak limit, i.e. when χ = F 0 φ, where ∇χ is the vacuum magnetic field, and φ is the toroidal angle. When using the new model, energy is conserved almost exactly, with a small amount of error due to temporal discretization and the use of a truncated Fourier decomposition in the toroidal direction. However, the new model also has more error in the momentum conservation than the original model whenever v = 0.

Appendix A.
In most reduced MHD models used for tokamaks, including the tokamak reduced MHD model implemented in JOREK (Franck et al. 2015;Hoelzl et al. 2021), field compression is neglected altogether via the ansatz B = F 0 ∇φ + ∇Ψ × ∇φ, which assumes that the background vacuum field is purely toroidal and only allows field line bending corrections to that. This is possible due to a property of the projection operator used, ∇φ · ∇ × (R 2 * , which allows exact force balance in the reduced system of equations when equilibria satisfy the Grad-Shafranov equation, even if field compression is neglected. Consider the second (poloidal momentum) equation in section 2.6 of Franck et al. (2015) with u = v = 0 †: where R is the distance from the axis of symmetry (major radius coordinate), p is pressure, Ψ is the poloidal flux and j = ∆ * Ψ . One can easily see that this equation can be obtained by applying the projection operator ∇φ · ∇ × (R 2 to the static equilibrium condition j × B = ∇p. Under the assumption of axisymmetry, the derivative in the last term in equation (A 1) is zero. In addition, since pressure is a flux function, we have ∇p = p ∇Ψ . Thus, equation ( and can be obtained by projecting j × B = ∇p onto B. In this equation, the first term is again identically zero due to p being a flux function, and the second term is zero due to axisymmetry. Thus, any solution of the Grad-Shafranov equation is also an exact solution of the reduced MHD system with v = 0, and so this version of reduced MHD has the same set of axisymmetric equilibria as full MHD. As is clear in the proof above, the Grad-Shafranov equation plays a crucial role. Without it, the reduced MHD system would not admit a solution with v = 0 even if we were to find a (Ψ, p) satisfying j × B = ∇p. In other words, for three dimensional equilibria, one should not expect to have force balance in the reduced MHD system. While not proven, it seems unlikely to us that a different choice of projection operator could resolve this problem. At high aspect ratios, lack of force balance becomes less significant due to the neglected field compression becoming smaller, and vanishing completely in the cylindrical limit. However, for an aspect ratio of 10, the lack of force balance is already significant. In order to carry out the simulations presented in section 4 with our original model, we had to introduce a static force balancing term, which was not considered in Nikulsin et al. (2019). For the sake of clarity, we consider the force balancing term before the projection operator is applied, however note that the full MHD vector momentum equation is overconstrained if the velocity ansatz is used but the projection operator is not applied. The modified equation is as follows: where B is the magnetic field given by the ansatz, µ 0 j = ∇ × B and B f and j f are the full MHD magnetic field and current, i.e. the actual field and current, not their ansatzbased approximations. At t = 0, the two ansatz-based Lorentz force terms cancel, and the forces are balanced by the full MHD Lorentz force; the first Lorentz force term then evolves with time while the force balancing term remains frozen. Such an approach should work in most cases, as instabilities tend to not compress the magnetic field due to the large amount of energy required to do so. In the tokamak limit, the force balancing term reduces to −F F ∇Ψ | t=0 /R 2 .
We can show that this extra term does not violate global conservation of momentum. If we ignore the viscosity term, equation (A 5) can be rewritten as already dealt with. Most importantly, if ∂ Ψ = ∂ g ik = 0, then j ∇χ, where g ik are the components of the metric tensor of the vacuum field-aligned coordinate system. The conditions for exact momentum conservation are thus ∂ u = 0, u ∈ {g ik , Φ, Ψ, v , p, ρ, P }.

(A 8)
The only minor difference from condition (5.2) is the replacement of p with p.