## 1. Introduction

Reduced magnetohydrodynamics (MHD), a system of approximations introduced in its original form in the 1960s (Greene & Johnson Reference Greene and Johnson1961), continues to be used by modern MHD codes, such as JOREK (Franck *et al.* Reference Franck, Hölzl, Lessig and Sonnendrücker2015) and M3D-$C^1$ (Breslau, Ferraro & Jardin Reference Breslau, Ferraro and Jardin2009). 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 Reference Strauss1997; Jardin *et al.* Reference Jardin, Ferraro, Breslau and Chen2012). 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 (Kruger, Hegna & Callen Reference Kruger, Hegna and Callen1998; Jardin *et al.* Reference Jardin, Ferraro, Breslau and Chen2012). 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 using different orderings that keep different physical effects (e.g. low-$\beta$ versus high-$\beta$ orderings) (Kruger *et al.* Reference Kruger, Hegna and Callen1998; Strauss Reference Strauss1976, Reference Strauss1977, Reference Strauss1980, Reference Strauss1997). An alternative ansatz-based approach was introduced by Park *et al.* (Reference Park, Monticello, White and Jardin1980), where an ansatz form that eliminates fast magnetosonic waves is used for the velocity and terms of all orders are kept. Later papers also adopt an ansatz for the magnetic field that eliminates field compression (Huysmans & Czarny Reference Huysmans and Czarny2007; Breslau *et al.* Reference Breslau, Ferraro and Jardin2009; Franck *et al.* Reference Franck, Hölzl, Lessig and Sonnendrücker2015). The ansatz-based approach makes fewer assumptions and keeps more physical effects, while generally resulting in more complicated equations.

An ordering-based reduced MHD model for stellarators was recently derived in Zocco, Helander & Weitzner (Reference Zocco, Helander and Weitzner2021) with $\beta \sim \epsilon \ll 1$, where $\epsilon$ is the inverse aspect ratio. Here, as in our previous paper (Nikulsin *et al.* Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019), we adopt an ansatz-based approach. Previously, we derived a hierarchy of models, which we intended to use for stellarator simulations. We started by introducing ansatzes for the magnetic field and velocity and proving that any arbitrary magnetic field and velocity can be represented by those ansatzes. We also showed that if the background vacuum field is strong enough, then the three terms of velocity ansatz approximately separate the three MHD waves. Dropping the fast magnetosonic wave term from the velocity ansatz and the field compression term from the magnetic field ansatz, we obtain reduced MHD (Nikulsin *et al.* Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019). However, we later found that the projection operators used had good momentum conservation properties, but not as good energy conservation properties. We thus adopt a different set of projection operators, which makes our new stellarator models a direct generalization of the currently used JOREK tokamak models to three-dimensional geometries; the stellarator models reduce back to the tokamak models in the tokamak limit.

The rest of this paper is structured as follows. In § 2, we briefly review the original set of models, introduce the changes that we made and derive the new set of models. In § 3, we show how the deficiencies of our original model are fixed in the new model. In § 4, we provide some test runs, which show how the deficiencies of the original model manifest themselves in practice by comparing it with the new stellarator model/standard tokamak model (both are the same in this situation) in the case of a tearing mode in a tokamak. In § 5, we consider analytically the local momentum conservation properties of both the original and new models, and then compare the global momentum conservation of the new stellarator model with and without field-aligned flow in the case of a ballooning mode. Finally, in Appendix A, a technicality in the implementation of the original model is discussed.

## 2. Changes made to the model

The original hierarchy of models, as presented in Nikulsin *et al.* (Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019), was derived from the viscoresistive MHD equations:

The gradient operators parallel and perpendicular to the total magnetic field $\boldsymbol {B}$ are defined as $\boldsymbol {\nabla }_\parallel = ({\boldsymbol {B}}/{B^2})\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }$ and $\boldsymbol {\nabla }_\perp = \boldsymbol {\nabla } - \boldsymbol {\nabla }_\parallel$. The usual MHD notation is followed with $\rho$, $p$, $\boldsymbol {v}$ and $\boldsymbol {B}$ being density, pressure, velocity and magnetic field, respectively. In addition to that, $\eta$ is the resistivity, $\mu$ is the dynamic viscosity, $D_\perp$ is the mass diffusion coefficient perpendicular to field lines, $\kappa _\perp$ and $\kappa _\parallel$ are the thermal conductivity across and along field lines, respectively, and $S_\rho$ and $S_e$ are source terms in the continuity and energy equations, respectively. The ideal gas law $p = \rho 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, $\boldsymbol {\nabla }\chi$ is the curl-free vacuum component of the magnetic field, which is generated by the coils, $B_v = |\boldsymbol {\nabla }\chi |$ and $\boldsymbol {\nabla }^{\perp } = \boldsymbol {\nabla } - B_v^{-2}\boldsymbol {\nabla }\chi (\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla })$. 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 $\varOmega = \zeta = 0$ and dropping the evolution equations for $\varOmega$ and $\zeta$ (derived below). Further, in the tokamak limit, i.e. when $\chi = F_0\phi$, where $\phi$ is the toroidal angle,Footnote ^{1} the ansatzes reduce to those of the standard JOREK reduced MHD tokamak models, $\boldsymbol {B} = F_0\boldsymbol {\nabla }\phi + \boldsymbol {\nabla }\psi \times \boldsymbol {\nabla }\phi$ and $\boldsymbol {v} = R^2\boldsymbol {\nabla } u\times \boldsymbol {\nabla }\phi + v_{\parallel }\boldsymbol {B}$, where $\psi = F_0\varPsi$ and $u = \varPhi /F_0$ (Hoelzl *et al.* Reference Hoelzl, Huijsmans, Pamela, Becoulet, Nardon, Artola, Nkonga, Atanasiu, Bandaru and Bhole2021).

The scalar $\psi _v$ is one of the Clebsch potentials for the vacuum field, which can locally be written as

As we discussed in Nikulsin *et al.* (Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019), one can construct a Clebsch-type coordinate system with the coordinates $(\psi _v,\beta _v,\chi )$. Note that the variables $\psi _v$ and $\beta _v$ are not unique, as can be seen by replacing $\psi _v$ with $\psi _v' = \psi _v + f(\beta _v)$, where $f$ is an arbitrary function, which leaves (2.3) unchanged. In general, given some vacuum magnetic field $\boldsymbol {\nabla }\chi$, the scalars $\psi _v$ and $\beta _v$ are completely determined by fixing some parameterization $(\psi _v,\tilde {\beta }_v)$ of a surface intersected by the vacuum field lines, with the values of $\psi _v$ and $\tilde {\beta }_v$ in the rest of space being determined by following the field lines off of the surface while keeping $\psi _v$ and $\tilde {\beta }_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 $\psi _v$ and $\tilde {\beta }_v$ when the field lines travel once around the torus and encounter previously labelled points. Finally, $\beta _v$ is determined by integrating $\partial \beta _v/\partial \tilde {\beta }_v = B_v/|\boldsymbol {\nabla }\psi _v\times \boldsymbol {\nabla }\tilde {\beta }_v|$ (Stern Reference Stern1970; D'haeseleer *et al.* Reference D'haeseleer, Hitchon, Callen and Shohet2012). Usually, it is advantageous to find a $\psi _v$ such that the component of $\boldsymbol {\nabla }\varOmega \times \boldsymbol {\nabla }\psi _v$ perpendicular to $\boldsymbol {\nabla }\chi$, $F_v^2|\partial ^{\parallel }\varOmega |^2$, is minimized at $t = 0$, so that most of field line bending is captured by the $\boldsymbol {\nabla }\varPsi \times \boldsymbol {\nabla }\chi$ term.

The variables $\varPsi$, $\varOmega$, $\varPhi$, $v_{\parallel }$ and $\zeta$ 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 $\boldsymbol {\nabla }\psi _v$ and $\boldsymbol {\nabla }\chi$, and applied the following projection operators to the vector momentum equation, which was first divided by $\rho$:

Here, $\boldsymbol {e}_\chi = \boldsymbol {B}/B^\chi$ is the covariant basis vector in the Clebsch-type coordinate system $(\alpha ,\beta ,\chi )$ associated with the full magnetic field, $\boldsymbol {B} = \boldsymbol {\nabla }\alpha \times \boldsymbol {\nabla }\beta$. 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.

where $\varDelta ^{\perp } = \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\nabla }^{\perp }$. However, for the sake of energy conservation, a better set of projection operators would be

as we 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.6), this time without dividing first by $\rho$. The momentum equation can be written as

where we used the identity $(\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v} = \frac {1}{2}\boldsymbol {\nabla } v^2 + \boldsymbol {\omega }\times \boldsymbol {v}$ and $\boldsymbol {\omega }$ is the vorticity:

Just as in Nikulsin *et al.* (Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019), we do not carry the viscosity term through this derivation, instead adding a generic viscosity term afterwards. To obtain the evolution equation for $\varPhi$, we apply to (2.7) the first projection operator in (2.6), or its equivalent, $-\boldsymbol {\nabla }\boldsymbol {\cdot }(B_v^{-2}\boldsymbol {\nabla }\chi \times$, when appropriate:

where $[f,g] = B_v^{-1}\boldsymbol {\nabla }\chi \boldsymbol {\cdot }(\boldsymbol {\nabla } f\times \boldsymbol {\nabla } g)$ is a Poisson bracket, $\partial ^{\parallel } = B_v^{-1}\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla }$, superscripts indicate contravariant vector components and, for any vector $\boldsymbol {U}$, $\boldsymbol {U}^\perp = \boldsymbol {U} - U^\chi \boldsymbol {\nabla }\chi /B_v^2$. We have added the generic viscosity term ${\boldsymbol {\nabla }\boldsymbol {\cdot }(\mu _\perp \boldsymbol {\nabla }^{\perp }\varDelta ^{\perp }\varPhi )}$, which we have chosen to match the viscosity term in Franck *et al.* (Reference Franck, Hölzl, Lessig and Sonnendrücker2015) in the tokamak limit. As discussed in previous papers (Franck *et al.* Reference Franck, Hölzl, Lessig and Sonnendrücker2015; Nikulsin *et al.* Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019), the viscosity term $\mu \Delta \boldsymbol {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 lose much by using a generic viscosity term in the final equations. Note, however, that the $\mu _\perp$ in the above equation is not the physical viscosity $\mu$ in (2.1). Indeed, from dimensional analysis it follows that $\mu _\perp$ has units of T$^2$ Pa s. Applying scaling analysis to the term $\mu \Delta \boldsymbol {v}$ after inserting just the first term of the ansatz (2.2) for $\boldsymbol {v}$ and applying the first projection operator in (2.6), we see that $\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla }\times [B_v^{-2}\mu \varDelta (\boldsymbol {\nabla }\varPhi \times \boldsymbol {\nabla }\chi /B_v^2)] \sim \mu B_v^{-2}\varPhi /L_\perp ^4$, where $L_\perp$ is the length scale perpendicular to the magnetic field, and $L_\perp \ll L_\parallel$ for most magnetic fusion configurations. Applying the same analysis to the generic viscosity term, we get $\boldsymbol {\nabla }\boldsymbol {\cdot }(\mu _\perp \boldsymbol {\nabla }^{\perp }\varDelta ^{\perp }\varPhi ) \sim \mu _\perp \varPhi /L_\perp ^4$. Comparing the two scalings, we see that $\mu _\perp$ scales as $\mu B_v^{-2}$; for the purposes of the scaling, one can take the value of $B_v$ at the axis as a typical value and write $\mu _\perp \sim \mu B_{v,\textrm {axis}}^{-2}$.

To get the evolution equation for $v_{\parallel }$, we apply the second operator in (2.6), projecting equation (2.7) on $\boldsymbol {B}$:

where we have again added a generic viscosity term with a form analogous to the viscosity term in the previous equation. Here, $(f,g) = \boldsymbol {\nabla }^{\perp } f\boldsymbol {\cdot }\boldsymbol {\nabla }^{\perp } g$ is an inner product, $F_v = |\boldsymbol {\nabla }\psi _v|$, $\partial _{\psi _v} = F_v^{-1}\boldsymbol {\nabla }\psi _v\boldsymbol {\cdot }\boldsymbol {\nabla }$ and $[f,g]_{\psi _v} = F_v^{-1}\boldsymbol {\nabla }\psi_v \boldsymbol {\cdot }(\boldsymbol {\nabla } f\times \boldsymbol {\nabla } g)$ is a Poisson bracket with respect to $\boldsymbol {\nabla }\psi _v$. Finally, to obtain the evolution equation for $\zeta$, we apply the last projection operator in (2.6) (or its equivalent, ${-\boldsymbol {\nabla }\boldsymbol {\cdot }[B_v^{-2}\boldsymbol {\nabla }\chi \times (\boldsymbol {\nabla }\chi \times }$, to some terms) to equation (2.7):

This concludes the derivation of the scalar momentum equations.

The evolution equations for the magnetic field variables of Nikulsin *et al.* (Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019) are valid; therefore we keep them as they are:

The continuity equation remains unchanged (Nikulsin *et al.* Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019):

while for the pressure, we use the standard internal energy evolution equation, instead of the full energy conservation equation from the system (2.1):

which now, after using (2.2), reads

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 term containing slow magnetosonic waves and the last term containing fast magnetosonic waves (Nikulsin *et al.* Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019). Thus, to get the reduced MHD system, we set $\zeta = \varOmega = 0$, eliminating fast magnetosonic waves and field compression, and drop the corresponding evolution equations, namely (2.11) and the second equation in (2.12), to avoid having an overconstrained system. We also approximate the electric field in Ohm's law as $\boldsymbol {E} = -\boldsymbol {v}\times \boldsymbol {B} + \eta \,\boldsymbol {j}^\parallel$, which corresponds to neglecting the components of $\boldsymbol {j}$ perpendicular to $\boldsymbol {\nabla }\chi$ in the last term of the first equation in (2.12). Alternatively, if we introduce a tensor resistivity, $\boldsymbol {\eta }$, this approximation corresponds to neglecting the perpendicular resistivity, $\eta ^\perp \approx 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 Reference Pfirsch1978), and likely also for any configuration with significant bootstrap current. Indeed, we can write $\boldsymbol {\eta }\boldsymbol {\cdot }\boldsymbol {j} = \eta _\perp \boldsymbol {j}_\perp + \eta _\parallel \,\boldsymbol {j}_\parallel$, where subscripts refer to parallel and perpendicular components relative to the total field $\boldsymbol {B}$, not just $\boldsymbol {\nabla }\chi$. When the bootstrap current is high enough, the $\eta _\parallel \boldsymbol {j}_\parallel$ term is dominant, and, since $\boldsymbol {\nabla }\chi$ is the dominant component of $\boldsymbol {B}$, we have $\boldsymbol {j}_\parallel \approx \boldsymbol {j}^\parallel$ in the lowest order. As we 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.15) with $\eta (\kern1.5pt j^\parallel )^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_{\parallel } = 0$ and dropping equation (2.10).

Finally, we point out that in the tokamak limit, when $\chi = F_0\phi$, the new reduced models derived in this section match the currently implemented JOREK tokamak models (Hoelzl *et al.* Reference Hoelzl, Huijsmans, Pamela, Becoulet, Nardon, Artola, Nkonga, Atanasiu, Bandaru and Bhole2021). The $\varPsi$ evolution equation, once we set $\zeta = \varOmega = 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 $\chi = F_0\phi$ (Hoelzl *et al.* Reference Hoelzl, Huijsmans, Pamela, Becoulet, Nardon, Artola, Nkonga, Atanasiu, Bandaru and Bhole2021). 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 $\chi = F_0\phi$.

## 3. Energy conservation

We begin by showing that the new models conserve energy. If we apply the first projection operator in (2.6) to some vector field $\boldsymbol {Q}$, multiply the result by a test function $\varPhi ^*$ and integrate over the plasma volume, we can get the following expression using the identity $\boldsymbol {\nabla } f\boldsymbol {\cdot }\boldsymbol {\nabla }\times \boldsymbol {U} = -\boldsymbol {\nabla }\boldsymbol {\cdot }(\boldsymbol {\nabla } f\times \boldsymbol {U})$ and integration by parts:

where we let $\varPhi ^* = 0$ on $\partial V$, so the surface integral term is zero. Doing the same with the third projection operator and test function $\zeta ^*$, we have

We can also apply the second projection operator (2.6) to $\boldsymbol {Q}$, multiply by a test function $v_{\parallel }^*$ and integrate over the plasma volume:

Now, if we let $\boldsymbol {Q}$ be the vector momentum equation, $\boldsymbol {Q} \equiv \rho \partial \boldsymbol {v}/\partial t + \rho (\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v} + \boldsymbol {v} P - \boldsymbol {j}\times \boldsymbol {B} + \boldsymbol {\nabla } p = 0$, as well as $\varPhi ^* = \varPhi$, $v_{\parallel }^* = v_{\parallel }$ and $\zeta ^* = -\zeta$, and sum up the above three equations, we get

where the left-hand side is zero due to equations (2.9), (2.10) and (2.11). More importantly, if we set $\zeta = 0$ and drop equation (2.11), (3.4) still holds, since with $\zeta ^* = -\zeta = 0$ equation (3.2) becomes $0 = 0$, and (2.11) no longer has to be satisfied in order for the left-hand side of (3.4) to remain zero. Similarly, if we set $v_{\parallel } = 0$ and drop equation (2.10), (3.4) continues to hold due to equation (3.3) becoming $0 = 0$ and no longer requiring equation (2.10) 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 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.12) are satisfied. These two equations are simply the vector components of Faraday's law in the $\boldsymbol {\nabla }\psi _v$ and $\boldsymbol {\nabla }\chi$ 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 $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B} = 0$ constraint. In the reduced MHD case, Faraday's law is also satisfied, as one can easily see from multiplying equation (2.18) by $\boldsymbol {\nabla }\chi /B_v$ and taking the curl. Thus, we can take the dot product of Faraday's law with $\boldsymbol {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^\parallel$ in the reduced MHD case. To complete the derivation of the energy conservation law, we rewrite (3.4) as

where we used the continuity equation. We then divide (2.14) by $\gamma -1$ and integrate it, together with (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 $\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {n} = 0$ and $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {n} = 0$, where $\boldsymbol {n}$ is a unit normal vector to the boundary. These conditions, along with $\varPhi = \zeta = 0$ on $\partial V$, which we had assumed before to make the boundary integral terms in (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 $\partial \rho /\partial t + \boldsymbol {\nabla }\boldsymbol {\cdot }(\rho \boldsymbol {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 $\varPhi ^*$, $\zeta ^*$ and $v_{\parallel }^*$ are finite-element basis functions. Since the finite-element solutions $\varPhi$, $\zeta$ and $v_{\parallel }$ are linear combinations of the basis functions, (3.6) will continue to hold. Similarly, since the $\varPsi$ and $\varOmega$ solutions will be linear combinations of the basis functions and the Bubnov–Galerkin method ensures that the weak form of (2.12) is satisfied whenever the test function is a basis function, the integral of (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 (2.14) 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 § 4. The practical accuracy of energy conservation in the standard JOREK tokamak model, which is what the new model presented in § 2 reduces to in the tokamak limit, is also considered in Hoelzl *et al.* (Reference Hoelzl, Huijsmans, Pamela, Becoulet, Nardon, Artola, Nkonga, Atanasiu, Bandaru and Bhole2021).

The necessity of neglecting perpendicular current in the resistive term of the first equation (2.12) when using the reduced MHD model can be clarified by deriving a magnetic field equation for the case when it is not neglected. Let $\varPsi _1$ and $\varPsi _2$ satisfy the equations

then $\varPsi = \varPsi _1 + \varPsi _2$ satisfies the original first equation in (2.12), without the neglect of perpendicular current. The first of the two equations above is equivalent to ${\partial \boldsymbol {B}_1/\partial t = -\boldsymbol {\nabla }\times \boldsymbol {E}_1}$, where $\boldsymbol {B}_1 = \boldsymbol {\nabla }\chi + \boldsymbol {\nabla }\varPsi _1\times \boldsymbol {\nabla }\chi$ and $\boldsymbol {E}_1 = -\boldsymbol {v}\times \boldsymbol {B} + \eta \boldsymbol {j}^\parallel$. Dotting with $\boldsymbol {B} = \boldsymbol {B}_1 + \boldsymbol {B}_2$ and adding $\boldsymbol {B}\boldsymbol {\cdot }\partial \boldsymbol {B}_2/\partial t + \boldsymbol {\nabla }\boldsymbol {\cdot }(\boldsymbol {E}_2\times \boldsymbol {B})$ to both sides, where $\boldsymbol {B}_2 = \boldsymbol {\nabla }\varPsi _2\times \boldsymbol {\nabla }\chi$ and $\boldsymbol {E}_2 = \eta \boldsymbol {j}^\perp$, we get

The evolution equation for $\varPsi _2$ can be written as $(\partial \boldsymbol {B}_2/\partial t)^{\psi _v} = -(\boldsymbol {\nabla }\times \boldsymbol {E}_2)^{\psi _v}$, and so the last two terms on the right-hand side of the above equation will only cancel partially, leaving behind the non-conservative terms $B_{\beta _v}\partial B_2^{\beta _v}/\partial t + B_{\beta _v}(\boldsymbol {\nabla }\times \boldsymbol {E}_2)^{\beta _v} + (\boldsymbol {\nabla }\times \boldsymbol {E}_2)^\chi$.

Finally, we show why the projection operators used in our original model allowed for non-physical 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.* (Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019), except using the projection operators (2.6) instead of the operators (2.4). Note that the projection operators (2.6) are orthogonal to all but one term in the covariant velocity ansatz (3.11):

Now, applying the first projection operator in (2.4) to a vector field $\boldsymbol {Q}$ and then multiplying by a test function $\varPhi ^*$ and integrating, we get

where we again used integration by parts to get the right-hand side. Similarly, for the third projection operator in (2.4), we can write

and for the second projection operator, we simply multiply by $v_{\parallel }^*$ and integrate, without doing any transformations:

In order to get $\int _V \rho \boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {Q}\,\textrm {d}V = 0$ when we set $\boldsymbol {Q} \equiv \partial \boldsymbol {v}/\partial t + (\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v} + \boldsymbol {v} P/\rho - \boldsymbol {j}\times \boldsymbol {B}/\rho + \boldsymbol {\nabla } p/\rho = 0$, we need $\varPhi ^* = -\tilde {\varPhi }$, $v_{\parallel }^* = \widetilde {v_{\parallel }}$, $\zeta ^* = \tilde {\zeta }$ and $\rho = \mathrm {const}$. The first three conditions pose no problem in the full MHD case, but if we try to set $\zeta = 0$, it does not imply $\tilde {\zeta } = 0$, and thus when we drop the evolution equation for $\zeta$, (3.14) 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 $\zeta \neq 0$. Thus, both the reduction and spatially varying density can lead to non-physical generation of kinetic energy.

## 4. 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 stellarator-capable models, as presented in Nikulsin *et al.* (Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019), with the standard JOREK tokamak model without parallel flow (Hoelzl *et al.* Reference Hoelzl, Huijsmans, Pamela, Becoulet, Nardon, Artola, Nkonga, Atanasiu, Bandaru and Bhole2021). Note that in the tokamak limit $\chi = F_0\phi$, the latter is equivalent to the the model without parallel flow from the set of new stellarator-capable models, as introduced in § 2. In addition, we point out that the standard tokamak reduced MHD model in JOREK is known to accurately reproduce tearing modes when compared with full MHD (Haverkort *et al.* Reference Haverkort, de Blank, Huysmans, Pratt and Koren2016; Pamela *et al.* Reference Pamela, Bhole, Huijsmans, Nkonga, Hoelzl, Krebs and Strumberger2020), thus comparing with the standard tokamak model is sufficient.

In the test cases considered below, $F_0 = 10\ \mathrm {T}\ \textrm {m}$. The initial conditions were set up by solving the Grad–Shafranov equation with $FF'(\psi _n) = 1.173(1-\psi _n)$ in units of $\mathrm {T}^2\ \textrm {m}$, where $\psi _n = (\psi - \psi _{\textrm {axis}})/(\psi _{\textrm {edge}}-\psi _{\textrm {axis}})$, and $p(\psi _n)=0$. In the tokamak limit, the $\psi$ in the Grad–Shafranov equation is related to the $\varPsi$ in the magnetic field ansatz (2.2) by $\psi = F_0\varPsi$. The density was initialized to be constant at $3.346\times 10^{-7}\ \mathrm {kg}\ \mathrm {m}^{-3}$, corresponding to $10^{20}$ deuterium ions per cubic metre, and $\varPhi$ was initialized to zero. A viscosity of $\mu = 5.159\times 10^{-8}\ \mathrm {kg}\ (\mathrm {m}\ \mathrm {s})^{-1}$ was used in all simulations, which corresponds to $10^{-5}$ in JOREK units for $\mu _\perp ^t$ in the standard tokamak model.Footnote ^{2} When using the original stellarator model, the kinematic viscosity was set to be $\nu = \mu /\rho _0 = 0.1542\ \mathrm {m}^2\ \mathrm {s}^{-1}$ ($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 $\boldsymbol {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:

where $\delta \boldsymbol {U}^n = \boldsymbol {U}^{n+1} - \boldsymbol {U}^n$. For more details on the numerical implementation, see Czarny & Huysmans (Reference Czarny and Huysmans2008), Huysmans & Czarny (Reference Huysmans and Czarny2007) and Hoelzl *et al.* (Reference Hoelzl, Huijsmans, Pamela, Becoulet, Nardon, Artola, Nkonga, Atanasiu, Bandaru and Bhole2021).

To begin with, we tested the original stellarator model (Nikulsin *et al.* Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019) in the linear regime by calculating the tearing mode growth rates at different resistivities and comparing them with the growth rates obtained from the standard tokamak model (Hoelzl *et al.* Reference Hoelzl, Huijsmans, Pamela, Becoulet, Nardon, Artola, Nkonga, Atanasiu, Bandaru and Bhole2021). 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,\ldots ,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 $\beta =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.* Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019) performs decently in the linear regime and at $\beta =0$, two problems can arise after nonlinear saturation. First, as discussed in the previous section, non-physical 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 simulation with a resistivity of $10^{-5}$ JOREK units ($1.9382\times 10^{-5}\ \Omega \ \textrm {m}$). In figure 3, $-\textrm {d}E/\textrm {d}t$ is plotted and compared with the physical energy loss rate. We define $E = \int _V\mathcal {E}\,\textrm {d}V$, where $\mathcal {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 $-\textrm {d}E/\textrm {d}t$ 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 $\sim$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 $-\textrm {d}E/\textrm {d}t$ and physical energy loss rate for a simulation with artificial dissipation in the form of a hyperviscosity term (with $\nu _h = 1.08\times 10^{-3}\ \mathrm {m}^4\ \mathrm {s}^{-1}$, corresponding to $7\times 10^{-10}$ in JOREK units) in the evolution equation for $\varPhi$. 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 behaviour 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 $\varPsi$ evolution equation (discussed below).

One can introduce a finite pressure into the original stellarator model by using (2.15) 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 $\beta$ becomes non-negligible, likely because the pressure terms in the original stellarator model (Nikulsin *et al.* Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019) (not shown in table 1) contain either $\boldsymbol {\nabla }\rho$ or $\partial ^{\parallel } 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 $-\textrm {d}E/\textrm {d}t$ 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 § 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 $\varPsi$ evolution equation can often produce ill-conditioned time stepping matrices, resulting in numerical instabilities. In fact, the $\varPsi$ equation is what destabilizes the kinetic energy, causing it to explode in figure 4(*a*). If we replace the $\varPsi$ evolution equation in the original stellarator model (see table 1) by (2.18), 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 $\psi$ evolution equation in the standard tokamak model by the corresponding equation from the original stellarator model (table 1), thus testing the $\varPsi$ 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 slightly 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 $\varPsi$ evolution equation in the original stellarator model in table 1 is equivalent to equation (2.18) 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 non-physical gain of kinetic energy, or depositing energy into it if there is a non-physical 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-$\beta$ simulations. It is much better to use equation (2.15) 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 $\beta$ is not negligible.

## 5. 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, $\varPhi$ and $v_{\parallel }$, 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 first compare momentum conservation properties of the original model presented in Nikulsin *et al.* (Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019) with those of the new model derived in § 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 non-physical generation of momentum within the plasma due to approximations made in reduced MHD.

### 5.1. 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.* (Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019). The action of the first projection operator (2.4) on a vector $\boldsymbol {Q}$ can be written as $\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla }\times [\boldsymbol {\nabla }\chi \times ( \boldsymbol {e}_\chi \times \boldsymbol {Q})] = \boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla }\times (\boldsymbol {e}_\chi Q^\chi - \boldsymbol {Q})$. Thus, if $\boldsymbol {Q}$ is the full MHD momentum equation (2.1), the action of the projection operator corresponds to dropping the two vector components of a vorticity-like equation, $\boldsymbol {\nabla }\times (\boldsymbol {e}_\chi Q^\chi - \boldsymbol {Q})$, that are perpendicular to $\boldsymbol {\nabla }\chi$. 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 $\boldsymbol {\nabla }\chi$. This equation can be written as

where the reduced velocity $\hbox{$\boldsymbol{\rlap{-} v}$} = -\boldsymbol {\nabla }\chi \times (\boldsymbol {e}_\chi \times \boldsymbol {v}) = \boldsymbol {\nabla }\varPhi \times \boldsymbol {\nabla }\chi /B_v^2$ and the reduced vorticity $\hbox{$\boldsymbol{\rlap{-} \omega}$} = \boldsymbol {\nabla }\times \hbox{$\boldsymbol{\rlap{-} v}$}$ were introduced, and the viscosity term is not considered. If the components of this equation perpendicular to $\boldsymbol {\nabla }\chi$ 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:

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 $\hbox{$\boldsymbol{\rlap{-} \omega}$}$ and $\boldsymbol {j}$ will be directed strictly along $\boldsymbol {\nabla }\chi$. If we allow $g^{ik}$, $\varPhi$ or $\varPsi$ to vary along $\boldsymbol {\nabla }\chi$, the same calculation will show that $\hbox{$\boldsymbol{\rlap{-} \omega}$}$ has non-zero perpendicular components. This will cause $\partial \hbox{$\boldsymbol{\rlap{-} \omega}$}/\partial t$ to have non-zero components perpendicular to $\boldsymbol {\nabla }\chi$, which cannot be cancelled by any other terms since there are no more time derivatives involving $\varPhi$ in the equation. Similarly, if we let any of the other quantities vary along $\boldsymbol {\nabla }\chi$, the last term and the seventh term on the right-hand side (pressure), the first term and the seventh term on the right-hand side (density), the last term on the left-hand side ($P$) and the third term on the left-hand side ($v_{\parallel }$) will have non-zero perpendicular components, which will not be cancelled by any other terms. If the conditions (5.2) are met, then only the sixth and eighth terms on the left-hand side have non-zero 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.6), the vorticity-like equation analogous to (5.1) can be written as

Equation (2.9) is just the contravariant $\chi$ component of this equation. In the above equation, the third and eleventh terms on the left-hand side will have perpendicular components even when condition (5.2) is met; namely $\boldsymbol {\nabla }^{\perp }(\rho B_v^{-2}\partial v _{\parallel }/\partial t)\times \boldsymbol {\nabla }\chi$ from the third term and $\boldsymbol {\nabla }^{\perp }(Pv_{\parallel } B_v^{-2})\times \boldsymbol {\nabla }\chi$ from the eleventh term. In addition, while the eighth term on the left-hand side of (5.4) is analogous to the sixth term on the left-hand side of (5.1), there is no term in (5.4) analogous to the eighth term of (5.1) to cancel the non-zero perpendicular components of its eighth term. Thus, a more stringent condition is required to ensure that the perpendicular components of equation (5.4) 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_{\parallel } = 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.* Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019) has better momentum conservation properties than the new reduced model (§ 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.

### 5.2. Global momentum conservation error (numerical examples)

As a test case, we use a simple ballooning mode in an X-point plasma in a tokamak with an aspect ratio of 3 (figure 6). In the simulations shown here only the Fourier modes with $0 \leqslant n \leqslant 6$ were included. The simulation parameters were as follows: $F_0 = 3\ \mathrm {T}\ \mathrm {m}$, a time step of 3 Alfvén times, a resistivity of $\eta = 3.8764\times 10^{-6}\ \Omega \ \mathrm {m}$ ($2\times 10^{-6}$ JOREK units) and a viscosity of $\mu = 2.293\times 10^{-7}\ \mathrm {kg}\ (\mathrm {m}\ \mathrm {s})^{-1}$, corresponding to $\mu _\perp ^t = 4\times 10^{-6}$ JOREK units. The initial conditions were set up by solving the Grad–Shafranov equation with

$T(\psi _n) = 0.015(1 - 0.66\psi _n)[1 - \tanh ((\psi _n - 0.94)/0.08)]/2 + 3\times 10^{-4}$ and $\rho (\psi _n)=[1-\tanh ((\psi _n - 0.94)/0.08)]/2 + 0.01$, where $FF'$ is in units of $\mathrm {T}^2\ \mathrm {m}$, $T$ is in $(10^{20}k_B\mu _0)$ K and $\rho$ is in $3.346\times 10^{-7}\ \mathrm {kg}\ \mathrm {m}^{-3}$. In addition, a hyperresistivity of $\eta _h = 5.8146\times 10^{-10}\ \Omega \ \mathrm {m}^{-2}$ ($3\times 10^{-10}$ JOREK units) was used in the simulations. We compare the new reduced model with $v_{\parallel } = 0$ with the new reduced model with $v_{\parallel }\neq 0$. Since in the tokamak limit these models match the standard JOREK tokamak models, we used the tokamak models for the actual 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 non-zero momentum tends to appear. Here, we choose the $x$ direction as the $R$ direction when $\phi =0$ and the $y$ direction as the $R$ direction when $\phi ={\rm \pi} /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 $\varPhi ^* = \ln R$ in (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_{\parallel }\neq 0$ is worse by more than an order of magnitude than that in the model with $v_{\parallel } = 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.

## 6. Conclusion

In continuation of our previous work (Nikulsin *et al.* Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019), 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 non-physical gain or loss of kinetic energy. Our original plan, as presented in Nikulsin *et al.* (Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019), 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-$\beta$ 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 non-physical 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.* Reference Franck, Hölzl, Lessig and Sonnendrücker2015; Hoelzl *et al.* Reference Hoelzl, Huijsmans, Pamela, Becoulet, Nardon, Artola, Nkonga, Atanasiu, Bandaru and Bhole2021) in the tokamak limit, i.e. when $\chi = F_0\phi$, where $\boldsymbol {\nabla }\chi$ is the vacuum magnetic field and $\phi$ 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_{\parallel }\neq 0$.

## Acknowledgements

The authors thank Rohan Ramasamy, Florian Hindenlang, Boniface Nkonga, Guido Huijsmans, Vinodh Bandaru and Erika Strumberger for fruitful discussions. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

*Editor Paolo Ricci thanks the referees for their advice in evaluating this article.*

## Declaration of interests

The authors report no conflict of interest.

## Appendix A

In most reduced MHD models used for tokamaks, including the tokamak reduced MHD model implemented in JOREK (Franck *et al.* Reference Franck, Hölzl, Lessig and Sonnendrücker2015; Hoelzl *et al.* Reference Hoelzl, Huijsmans, Pamela, Becoulet, Nardon, Artola, Nkonga, Atanasiu, Bandaru and Bhole2021), field compression is neglected altogether via the ansatz $\boldsymbol {B} = F_0\boldsymbol {\nabla }\phi + \boldsymbol {\nabla }\varPsi \times \boldsymbol {\nabla }\phi$, 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, ${\boldsymbol {\nabla }\phi \boldsymbol {\cdot }\boldsymbol {\nabla }\times (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 § 2.6 of Franck *et al.* (Reference Franck, Hölzl, Lessig and Sonnendrücker2015) with $u = v_\parallel = 0$:Footnote ^{3}

where $R$ is the distance from the axis of symmetry (major radius coordinate), $p$ is pressure, $\varPsi$ is the poloidal flux and $\hbox{${\rlap{-} j}$} = \varDelta ^*\varPsi$. One can easily see that this equation can be obtained by applying the projection operator $\boldsymbol {\nabla }\phi \boldsymbol {\cdot }\boldsymbol {\nabla }\times (R^2$ to the static equilibrium condition $\boldsymbol {j}\times \boldsymbol {B} = \boldsymbol {\nabla } p$. Under the assumption of axisymmetry, the derivative in the last term in (A1) is zero. In addition, since pressure is a flux function, we have $\boldsymbol {\nabla } p = p'\boldsymbol {\nabla }\varPsi$. Thus, (A1) becomes

Meanwhile, the Grad–Shafranov equation reads

Substituting this into the second term in (A2), we have

where the last two terms are identically zero since $p'$ and $FF'$ are both functions of only $\varPsi$, so their Poisson bracket with $\varPsi$ must be zero. Finally, the last (parallel momentum) equation in § 2.6 of Franck *et al.* (Reference Franck, Hölzl, Lessig and Sonnendrücker2015) with $u = v_\parallel = 0$ reads

and can be obtained by projecting $\boldsymbol {j}\times \boldsymbol {B} = \boldsymbol {\nabla } p$ onto $\boldsymbol {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 $\boldsymbol {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 $\boldsymbol {v} = 0$ even if we were to find a $(\varPsi ,p)$ satisfying $\boldsymbol {j}\times \boldsymbol {B} = \boldsymbol {\nabla } 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 § 4 with our original model, we had to introduce a static force balancing term, which was not considered in Nikulsin *et al.* (Reference Nikulsin, Hoelzl, Zocco, Lackner and Günter2019). 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 $\boldsymbol {B}$ is the magnetic field given by the ansatz, $\mu _0\boldsymbol {j} = \boldsymbol {\nabla }\times \boldsymbol {B}$ and $\boldsymbol {B}_f$ and $\boldsymbol {j}_f$ are the full MHD magnetic field and current, i.e. the actual field and current, not their ansatz-based 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 $-FF'\,\boldsymbol {\nabla }\varPsi |_{t=0}/R^2$.

We can show that this extra term does not violate global conservation of momentum. If we ignore the viscosity term, (A6) can be rewritten as

Integrating over the plasma volume and assuming that both the reduced and full magnetic fields are tangential to the plasma boundary, we obtain

The right-hand side is just the force exerted by the walls on the plasma. Note that the term $(B_f^2 - B^2)|_{t=0}/(2\mu _0)$ is approximately the energy density stored by compressing the background vacuum field in equilibrium. Since instabilities tend to not compress the magnetic field, the energy density stored in the compression should be roughly constant throughout the simulation, we should have $B^2 + (B_f^2 - B^2)|_{t=0} \approx B_f^2$, where $B_f^2/(2\mu _0)$ is the total magnetic energy density, which includes both field line bending and compression.

To conclude, we note that the momentum conservation properties discussed in § 5 remain unchanged with the addition of the force balancing term. To show this, we substitute $p = p|_{t=0} + \tilde {p}$ into (A6), where $\tilde {p}$ is defined as ${p - p|_{t=0}}$. Then, because $(\kern1.5pt\boldsymbol {j}_f\times \boldsymbol {B}_f)|_{t=0} = \boldsymbol {\nabla } p|_{t=0}$, those terms will cancel, and the only term remaining in (A6) that was not originally present in the vector momentum equation (2.1) is $-(\kern1.5pt\boldsymbol {j}\times \boldsymbol {B})|_{t=0}$. However, this is just the reduced Lorentz force term at $t=0$, which has the same structure as the time-varying reduced Lorentz force term that we have already dealt with. Most importantly, if $\partial ^{\parallel }\varPsi = \partial ^{\parallel } g^{ik} = 0$, then $\boldsymbol {j}\parallel \boldsymbol {\nabla }\chi$, 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

The only minor difference from condition (5.2) is the replacement of $p$ with $\tilde {p}$.