Extended full-MHD simulation of non-linear instabilities in tokamak plasmas

Non-linear magnetohydrodynamic (MHD) simulations play an essential role in active research and understanding of tokamak plasmas for the realization of a fusion power plant. The development of MHD codes such as JOREK is a key aspect of this research effort. In this paper, we present an operational version of the full-MHD model implemented in JOREK, a significant advancement from the reduced-MHD model used for previous studies, where assumptions were made on the perpendicular dynamics and the toroidal magnetic field. The final model is presented in detail, and benchmarks are performed using both linear and non-linear simulations, including comparisons between the new full-MHD model of JOREK and the previously extensively studied reduced-MHD model, as well as results from the linear full-MHD code CASTOR3D. For the cases presented, this new JOREK full-MHD model is numerically and physically reliable, even without the use of numerical stabilization methods. Non-linear modeling results of typical tokamak instabilities are presented, including disruption and edge-localized-mode physics, most relevant to current open issues concerning future tokamaks such as ITER and DEMO. VC 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http:// creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0018208


I. INTRODUCTION
Industrial electricity production using nuclear fusion power would greatly contribute to the reduction of greenhouse gas emissions and of long-lived radioactive nuclear waste, while providing electricity to society without the limit of an exhaustible natural resource. A favorable candidate for industrial fusion reactors is the tokamak device. Tokamaks use a helical magnetic field that winds itself around a toroidal vacuum chamber. The periodic nature of the torus ensures that charged particles, which approximately follow magnetic field lines, are not lost at the end of open field lines like in linear plasma devices. However, this periodicity can be subject to resonance and instabilities. Large-scale instabilities typically involve both the plasma and the magnetic field and are often studied in the fluid picture using magnetohydrodynamics (MHD).
There is a wide variety of MHD instabilities in tokamak plasmas, some of which can reduce or limit the operational capabilities of the machine. Edge-localized-modes (ELMs) are instabilities that eject plasma from the confined region onto the material surfaces of the first wall of the machine; these instabilities can lead to large heat-fluxes that may reduce the lifetime of the material surfaces. [1][2][3][4] Toroidal Alfv en eigenmodes (TAEs) can be excited by the 3.5 MeV alpha-particles born from fusion reactions; these can significantly limit the plasma pressure and are a concern for future reactors where the burning plasma will produce a large amount of alpha-particles. [5][6][7][8] Global MHD instabilities, during which the entire plasma is affected, can lead to disruptions; during such events, the kinetic and magnetic energies of the plasma can be transferred to the wall, leading to material heatfluxes and/or wall-currents that can damage the machine and its structure. [9][10][11][12][13][14] In order to study, understand, and predict these instabilities, non-linear MHD simulations are performed using numerical codes such as JOREK, 15,16 M3D-C1, 17,18 NIMROD, 19,20 XTOR, 21 BOUTþþ, 22,23 MEGA, 24,25 and HALO. 7 Previous studies of MHD instabilities with the JOREK code relied on a reduction of the full-MHD system. This reduction, known as reduced MHD, assumes that the toroidal magnetic field is constant in time and that the perpendicular velocity (perpendicular with respect to the magnetic field) is approximately poloidal. 26,27 The latter assumption is mostly kept for simplicity in formulating the equations. It does not enforce incompressibility of the fluid (i.e., r Áṽ 6 ¼ 0), and parallel dynamics is retained (i.e.,ṽ jj 6 ¼ 0), which is essential for a reliable description of tokamak instabilities. 28,29 The former assumption, however, leads to a physical simplification of the representation of tokamak dynamics. In particular, fixing the toroidal magnetic field in time removes some of the waves from the model, commonly referred to as fast magneto-sonic waves. These fast-waves propagate by compressing the magnetic field and the fluid in the direction perpendicular to the magnetic field. They are inherently stable with a wavelength and amplitude much smaller than the macro-instabilities of interest in tokamak devices. However, these waves can pollute and restrict the solution of numerical codes, and modeling codes often rely on distinct formulations or numerical stabilization to avoid numerical issues with these fast waves. Therefore, reduced-MHD is a powerful model in the sense that fast-waves are entirely removed from the system, and it is clearly sufficient to describe the non-linear dynamics of most tokamak instabilities.
However, it is already widely acknowledged that some tokamak instabilities, particularly internal kink instabilities, cannot properly be captured by reduced-MHD, 30 and so, including a full-MHD model (as opposed to reduced-MHD) in JOREK is necessary. In this paper, we present the full-MHD model implemented in JOREK and applied to linear and non-linear instabilities. A first implementation of the full-MHD model had already been presented in Ref. 31. Unfortunately, this initial study suffered implementation issues that restricted simulations to low-beta instabilities only. All numerical issues have now been resolved, such that modeling of all instabilities of interest is now possible. Section II of this paper presents the key numerical methods and the full-MHD physics model. Section III addresses the linear benchmark of the model for several types of MHD instabilities and how it compares to the reduced-MHD model. Section IV demonstrates the capability of the code to deal with large non-linear instabilities such as ELMs and disruptions. Finally, Sec. V summarizes the work and lays out the further improvements required for future studies of tokamak instabilities.

II. PHYSICS MODEL AND NUMERICAL METHODS
In this section, we describe the essential physical and numerical ingredients of the full-MHD model implemented in JOREK, including visco-resistive and diffusive effects, sources, diamagnetic rotation and neoclassical friction, boundary conditions, and normalization.

A. Numerical methods
Before introducing the physics model, some of the key numerical methods of the JOREK code are presented. These methods are the same as in previous publications, 16,32 but they are repeated here to point out their role in the numerical stability of the full-MHD model.
JOREK uses a finite-element grid in the poloidal plane, with finite Fourier series in the toroidal direction. The finite elements are quadrangular bi-cubic Bezier elements using a nodal formulation that ensures continuity of variables and their derivatives across elements (commonly defined as C1-or G1 continuity). 32 This is also valid for the poloidal space coordinates (R,Z) such that the elements are isoparametric. The Fourier series used in the toroidal direction can be chosen with arbitrary toroidal periodicity, meaning that with three Fourier modes and a toroidal periodicity of 4, the physical toroidal mode numbers simulated would be n ¼ 4, 8, and 12. Typically, when simulating linear stability scans, a single Fourier mode is used, while scanning the toroidal periodicity. For non-linear simulations, multiple Fourier modes are used, with a periodicity of 1 or 2.
The time discretization used in JOREK is fully implicit, with a choice between the Crank-Nicolson or the Gear scheme. In the simulations presented here, the Gear scheme 33 was used. These methods result in a linearized system of equations, for which a sparse matrix needs to be inverted. This is done either with a direct solve or using a pre-conditioned GMRES (Generalized Minimal Residual Solver) iterative solver. In practice, the direct solve approach is only used for the toroidally axisymmetric kinetic equilibrium n ¼ 0, while for non-linear cases, the GMRES method is employed. The GMRES pre-conditioner is equivalent to a Block-Jacobi pre-conditioner, which is obtained by solving each toroidal harmonic block of the matrix individually (i.e., without the harmonic coupling), which is done in parallel with a direct solve. The direct solve (either for the pre-conditioner harmonic-blocks or for the whole system) is done using open-source solvers such as MUMPS, 34 PASTIX, 35 or STRUMPACK. 36 The weak formulation method 37 is used for all equations, which are multiplied by a test-function and integrated over the element surfaces. One of the main advantages of the weak formulation is that it allows integration by parts of divergence terms, which allows the removal of all second order derivatives from the system. This is a significant advantage when using C1-continuous finite-elements such as those implemented in JOREK 32 since all terms in the equations are guaranteed to be smoothly represented, thus improving numerical stability. The weak Form method, similar to a variational form approach, has stabilizing properties, which are thought to be important in this study.
The equations are normalized using two constants: vacuum permeability l 0 and the central density q 0 , as in previous JOREK publications. 16 This normalization is similar to the Alfv en time normalization such that for a deuterium plasma with central particle density n o ¼ 6 Â 10 19 m À3 , a normalized time unit corresponds to approximately 0.5 ls. Note that this is a pseudo-normalization, where not all variables are dimensionless in the final formulation. In particular, the magnetic field is not normalized.

B. Basis of the full-MHD formulation
In order to ensure that the magnetic field satisfies Gauss's law r ÁB ¼ 0, we defined the magnetic field as Physics of Plasmas ARTICLE scitation.org/journal/php whereÃ is the magnetic vector potential and F ¼ FðwÞ is a toroidally axisymmetric equilibrium function, defined to satisfy the initial Grad-Shafranov equilibrium (w is the poloidal magnetic flux). F is constant in time, and thus, the evolution of the toroidal magnetic field is determined byÃ alone. Note: this differs from the reduced MHD definition of B / ¼ F 0 =R, where F 0 is constant both in space and time. It should be noted that this formulation of the magnetic field ensures r ÁB ¼ 0 exactly, without any approximation, since r ÁB includes derivatives ofÃ only up to second order, which are well defined in bicubic finite-elements. At the boundary between two elements, although second-order derivatives lose their continuity, the condition r ÁB ¼ 0 is still satisfied, from the left element's variables and from the right element's variables, individually. Additionally, the weak form method allows the formulation of all equations avoiding the second-order derivative. The other variables of the system are the velocity vector fieldṽ, the mass density q, and the temperature T. The resistive MHD model is described by @Ã @t ¼ṽ ÂB À gJ À rU; (2) q @ṽ @t ¼ Àqṽ Á rṽ þJ ÂB À rp; where g is the resistivity. The current is defined asJ ¼ r ÂB and the total pressure as p ¼ qT. The ratio of specific heats c is taken to be that of a monatomic gas, 5 3 .

C. Choice of gauge
In Eq. (2), U is the electric potential obtained from uncurling Faraday's law r Â ½@ tÃ ¼ ÀẼ. As the curl operator annihilates the gradient operator that acts on the electric potential, this information is lost in the curl operation so that a choice of Gauge is necessary to determine U. Since the magnetic field and electric field are invariant with respect to the transformationÃ 0 ¼Ã þ rW and U 0 ¼ U À @ t W (for any scalar function W), a convenient choice is to use Weyl's gauge, 38 where U 0 ¼ 0 is chosen such that which simplifies the induction equation (2) to give @Ã @t ¼ṽ ÂB À gJ : In practice, this means that if an external electric field is applied to the plasma, the magnetic vector potential will shift in time (even in the stationary equilibrium state).

D. Diffusion coefficients and sources
In addition to the resistivity g, physical diffusion coefficients are also included for all other variables. The perpendicular diffusion of density and temperature needs to be balanced in order to retain the initial Grad-Shafranov equilibrium, and thus, sources are also introduced for density and temperature. The visco-resistive MHD equations are, thus, written as @Ã @t ¼ṽ ÂB À gJ ÀS j ; q @ṽ @t ¼ Àqṽ Á rṽ þJ ÂB À rp þ lr 2ṽ À C qṽ ; where the particle diffusion and the thermal conductivity have been split into perpendicular and parallel components, D ? and D jj and j ? and j jj , respectively. The parallel and perpendicular gradient operators are defined as The current sourceS j keeps the current profile steady for long simulations.S j also includes the bootstrap current source, which evolves as a function of the pressure gradient and is particularly essential for cyclic instabilities such as ELMs. Instead of a current sourceS j , it is also possible to use a loop voltage as boundary conditions. The particle source S q and the heating source S T , which are generally used only to balance against axisymmetric equilibrium diffusive terms, can also be used for other purposes. A good example is to use a time-evolving particle source S q to model pellet injections. In Sec. IV, S q is used to mimic a massive-gas-injection (MGI) disruption. Note that since the continuity equation (10) has been used in the derivation of the momentum equation (9), the particle source term must also be kept there. If ignored, large particle sources such as pellets would also lead to an artificial injection of momentum. The diffusive terms from the continuity equation, although generally much smaller, are also kept in the momentum equation. This is represented by the last term in Eq. (9), where C q ¼ r Á ðD ? r ? q þ D jj r jj qÞ þ S q . The physical diffusive coefficients are non-constant. Radial profiles for D ? and j ? can be used to mimic various levels of cross field kinetic turbulent transport, which cannot be described by MHD. This is important for situations such as the H-mode, where cross field turbulent transport is known to be strongly reduced in the pedestal region, at the plasma edge. In simulations of ELMs, for example, a radial drop in D ? and j ? is localized in the pedestal region to sustain the large pedestal pressure gradient of pre-ELM conditions. There are two possible choices to maintain the background equilibrium profiles. Considering the density profile, for example, the perpendicular diffusion coefficient can be chosen as D ? ðwÞ $ D 0 ð@q=@wÞ À1 , like in Refs. 39 and 40, such that the density diffusive flux D ? rq is constant radially, meaning that a spatially constant source will sustain the profile. In that case, if the particle source is too large (or too weak), the density profile will, then, increase (or decrease) rigidly. Alternatively, an arbitrary diffusion profile can be chosen, in which case the particle source needs to be spatially adapted to compensate for radially varying Physics of Plasmas ARTICLE scitation.org/journal/php density losses. In most linear tests (as in Sec. III), this does not matter much since the background equilibrium is typically fixed (i.e., only toroidal modes n > 0 are evolved). However, for non-linear tests in Sec. IV, the first method is used to maintain profiles in a steady state, until the 3D non-linear perturbations are introduced. Both the viscosity and resistivity have a Spitzer-like dependence on temperature, g ¼ g 0 T À3=2 and l ¼ l 0 T À3=2 , respectively, where g 0 and l 0 are the values of resistivity and viscosity on the magnetic axis. The parallel particle diffusion is typically kept zero, while the parallel thermal conductivity is formulated using the Braginskii model j jj ¼ j 0 T 5=2 , for which j 0 has well-defined physical amplitudes for ion and electron temperatures. Since the present model only includes a single temperature, the value of j 0 is typically chosen to be the average of ion and electron temperature coefficients.
At present, numerical (or hyper-) diffusion is not needed in the model for any equation, even with strongly non-linear cases, as will be shown in Sec. IV. In future works, it may be adequate to introduce hyper-resistivity for physical reasons. 13

E. Extended MHD
For many tokamak applications, additional non-ideal effects are necessary to address certain physical properties of non-linear MHD instabilities. In particular, rotation effects can have a significant impact on the linear stability and the non-linear dynamics of MHD instabilities. The toroidal momentum induced by the neutral-beam injection (NBI) can be significant, particularly in small and medium-size spherical tokamaks like COMPASS, NSTX-U, and MAST-U, although this rotation becomes marginal in larger devices such as ITER. 41 The diamagnetic rotation 42 has a damping effect on the stability of high toroidal mode numbers and plays a major role in the dynamics of quasi-periodic relaxations, such as ELMs and neoclassical tearing modes (NTMs). [43][44][45][46][47][48] Finally, neoclassical poloidal rotation also plays an important role in the rotation of precursor modes in ELMy H-mode simulations. 49 The full system of extended MHD equations, including diffusion, sources, and rotation effects, is described by where the ion diamagnetic velocity and the neoclassical poloidal friction tensor are defined, respectively, as with m i and e being the ion mass and charge, respectively, k i being the neoclassical heat diffusivity, andẽ h ¼B h =jB h j being the unit vector along the poloidal magnetic fieldB h ¼ B RẽR þ B ZẽZ . The poloidal velocity is, thus, defined as v h ¼ ðṽ þṽ iÃ Þ Áẽ h . Note that factor 2 in the definition of the diamagnetic and neoclassical velocities comes from the assumption that, with a single total temperature T, the ion pressure is assumed to be simply p i ¼ p=2.
Note that the diamagnetic effects have been implemented taking into account the gyro-viscous cancellation and the gyro-viscous heatflux cancellation, which eliminates several diamagnetic terms in the momentum and energy equations. 46,47 This diamagnetic cancellation is the common form, which assumes a constant magnetic field, not the full (more complex) form where the cancellation involves the magnetization velocity instead of the diamagnetic velocity. 50 The diamagnetic term in the induction equation (12) is obtained using the standard drift ordering, as in Ref. 46, where the total electric field is reduced to the form by neglecting the perpendicular force-balance mi eq ðJ ÂB À r ? pÞ, which is a manner of neglecting the Hall effects. This is the form used in Eq. (12), assuming that p e ¼ p=2.
It is not clear at present how the drift ordering approximations (particularly the absence of Hall terms) affect the non-linear dynamics of tokamak instabilities. It is reasonable to assume that linear stability is not affected, since with unperturbed equilibrium quantities, the assumption thatJ ÂB ¼ r ? p should hold. However, in non-linear cases, where the magnetic field (and current) and pressure are strongly perturbed away from their equilibrium balance, the Hall effects may become important. However, a detailed study of the limitation of this assumption is beyond the scope of this work and will be addressed in future work.
With the full-MHD model described above, fast-waves can pollute the solution and lead to numerical instabilities. The diamagnetic velocity was found to be particularly challenging, with short-wave oscillations on the scale of grid elements. It was found that the diamagnetic terms from the rR and rZ projections of the induction equation (i.e., the A R and A Z equations) are the source of these numerical oscillations, and neglecting those terms can strongly improve numerical stability, while retaining the stabilizing effects of the diamagnetic terms, which are clearly dominated by the toroidal component of the induction equation. Work is currently on-going to implement numerical stabilization methods, such as Taylor-Galerkin (TG) and variational-multi-scale (VMS) stabilization, 51,52 which could improve the numerical stability of diamagnetic terms in the A R and A Z equations.

F. Projection and coordinate system
The momentum and induction equations each need to be projected in order to obtain individual equations for each of their vector-field components. Although the projection could simply be made along the orthogonal cylindrical basis vectors (ẽ R ;ẽ Z ;ẽ / ), this Physics of Plasmas ARTICLE scitation.org/journal/php is numerically not the most stable choice, particularly for the momentum equation.
The magnetic vector potential, the velocity field, and the magnetic field are decomposed as which means that the toroidally axisymmetric component of A / is equivalent to the magnetic poloidal flux w, and the toroidal components V / and B / are the actual physical components of the velocity and magnetic field, respectively. Note that if an electric field is present (for example, if a loop voltage is used as boundary conditions), then A / will shift rigidly with time, but the normalized scalar w n will not change (assuming that the equilibrium is stationary). While the induction equation is simply projected along the vectors (ẽ R ;ẽ Z ;ẽ / ), the momentum equation is projected along the vectors (ẽ R ;ẽ Z ;B). As noted in Ref. 31, this choice of projection is essential for the numerical stability of the model. By removing thẽ J ÂB term in the equation for V / , theB projection is thought to remove unnecessary fast magneto-sonic wave components that pollute the solution otherwise. Note that this projection does not entirely remove the fast magneto-sonic waves from the system, and they are still supported by the two other momentum equation projections and the induction equations. However, removing unnecessary information about the toroidal propagation of fast magneto-sonic waves from this equations seems favorable to the numerical stability of simulations. Separation of different waves is also used in other MHD formulations. 53,54

G. Boundary conditions
There are two types of boundary conditions, depending on whether a boundary surface is aligned to a magnetic flux-surface or not. If a boundary coincides with a flux-surface, then all variables are fixed in time (Dirichlet), although this can be relaxed for the density and temperature. If a boundary intersects magnetic field lines, then a mixture of conditions is applied: the Dirichlet condition is enforced for the magnetic vector potentialÃ, free outflow of density is allowed (Neumann), while Mach-1 and sheath boundary conditions are applied to velocity and temperature, respectively, such that whereb is the unit vector along the magnetic fieldB; c ¼ 5=3 is the ratio of specific heats, and c sh is the ion sheath transmission factor, which is typically taken between 4.5 and 10.0, depending on the ion or electron temperature and on transient energy fluxes. 55

III. COMPARISON OF FULL-MHD AND REDUCED-MHD
In this section, we address a series of linear benchmarks for tokamak instabilities to compare the new full-MHD model against the previously implemented reduced-MHD. Linear benchmarks are conducted for core MHD instabilities and edge peeling-ballooning instabilities, both in circular and X-point plasmas. A comparison of peeling-ballooning modes is provided for the spherical tokamak MAST, for which it is widely, and wrongly, presumed that reduced-MHD cannot describe ELM and filament physics accurately, even though no geometrical approximation is required for the energyconservative derivation of reduced-MHD. 56,57 A. Linear benchmark for core n ¼ 1 modes The first two linear benchmarks are a low-b m ¼ n ¼ 1 internal kink mode and a low-b m ¼ 2, n ¼ 1 tearing-mode. Both instabilities were studied in previous publications, and they are included here for completeness and clarity. Both equilibria are similar but differ in the q-profile: the internal kink mode has a q-profile in the range [0.7, 1.6], crossing q ¼ 1 at w n ¼0.5, while the tearing mode has a q-profile in the range [1.7, 3.9], crossing q ¼ 2 at w n ¼0.3. The Grad-Shafranov equilibrium quantities and profiles for these two cases are described in more detail in Ref. 31. The tearing mode is in a low-temperature regime (T i ¼ 0.25 eV) with a magnetic field of 1 T and a major radius of 10 m, with ion gyroradius q i ¼ 0:7 Â 10 À4 m and ion collision time Both cases are run for a scan in resistivity. The kink mode is run with resistivity alone (without viscosity and without particle or thermal diffusion), while the tearing mode is run including all diffusions, with l 0 ¼ 10 À8 kg m À1 s À1 ; D ? ¼ 0:7 m 2 s À1 , and j ? ¼ 1:7 Â10 À8 kg m À1 s À1 . For simplicity, the resistivity and viscosity are taken to be spatially constant for both cases to ease comparison with other codes. Note that only the toroidal mode n ¼ 1 is simulated here such that coupling with higher toroidal modes is not present in these linear benchmarks. Figures 1 and 2 show the benchmark of the internal kink mode and the tearing mode, respectively, compared to the reduced-MHD model. Poloidal cross sections of n ¼ 1 perturbed quantities are shown for the toroidal magnetic potential A / and the temperature (for the full-MHD model), and the growth rates of the modes are plotted as a function of resistivity, compared to the reduced-MHD model. A comparison is also made against CASTOR3D, which is also a full-MHD code. 59,60 Both cases agree with linear theory in the ideal (low resistivity) regimes, with a g 1=3 scaling for the internal kink mode 61 and a g 3=5 scaling for the tearing mode. 62 Although the agreement between reduced-MHD and full-MHD is reasonable for both cases, the reduced-MHD model starts to deviate from the full-MHD solution at low resistivity for the internal kink mode. This is a typical case where reduced-MHD is expected to fail: for internal kink modes at finite-b. 30 Although this is a low b N ¼ 0.4% case, reduced-MHD already seems to be affected. At higher-b, the deviation becomes more pronounced.

B. Where reduced-MHD fails
The internal kink mode at finite-b is a good example to illustrate why full-MHD is absolutely necessary for the modeling of some tokamak instabilities, particularly core MHD, which is highly relevant to disruptions, one of the main focuses of current research and experiments. In this benchmark, the results of the linear calculation of an n ¼ m ¼ 1 internal kink instability obtained by the reduced-MHD model of JOREK, the full-MHD model of JOREK, and the linear full-MHD code CASTOR3D 59,60 are compared for different values of b. The equilibrium characteristics of this circular plasma, together with the diffusive parameters used for the simulation, are as follows: Note that the viscosity here was used only for the JOREK runs, while it is zero for CASTOR3D. As shown in Fig. 3, the full-MHD model of JOREK and CASTOR3D agrees well on the linear growth rates, while the reduced-MHD model fails to reproduce these results at finite values of b. This shortcoming of the reduced-MHD model is expected due to the neglect of parallel magnetic field fluctuations. 30 Qualitative differences between the results of the reduced-and full-MHD models are also seen in nonlinear simulations of low mode number core instabilities.

C. Peeling-ballooning modes in a circular plasma
The second linear benchmark presented here is the so-called CBM18 case. It is a circular plasma unstable to peeling-ballooning modes, with b N % 1:3%. This case is run for individual toroidal mode

ARTICLE
scitation.org/journal/php numbers between n ¼ 1 and n ¼ 20. The resistivity for this test-case is set to g 0 ¼ 6 Â 10 À6 X m, while all other diffusion coefficients are set to zero. Figure 4 shows the benchmark of the CBM18 case, when compared to the reduced-MHD model, and also results from MISHKA and CASTOR, which are linear full-MHD codes. 59,63,64 MISHKA solves the linear incompressible ideal full-MHD model, while CASTOR solves the resistive full-MHD model (also using g 0 ¼ 6 Â 10 À6 X m). Poloidal cross sections of n ¼ 20 perturbed quantities are shown for the toroidal magnetic potential A / and the temperature (for the JOREK full-MHD model), and the growth rates of the modes are plotted as a function of the toroidal mode number. The stability threshold of n ¼ 4 is well reproduced by the full-MHD model, and the same asymptotic increase in growth rates is observed with both models. The agreement with reduced-MHD and with the linear MHD calculations is reasonable.

D. Peeling-ballooning modes in an X-point plasma
Next, an X-point plasma is run for peeling-ballooning instabilities. This is an artificial equilibrium similar to a JET plasma. It is run using the kinetic stationary background equilibrium flows (n ¼ 0), together with a single toroidal harmonic, which is changed from n ¼ 1 up to n ¼ 20. The equilibrium characteristics, together with the diffusive parameters, are as follows: where the D ? and j ? values are given on the magnetic axis, with a radial profile that drops by a factor 10 at the top of pedestal. Note that, in this list, b N is defined as b t =ða:B t :I p Þ, and Dw ped is given as a normalized quantity, relative to w n . Figure 5 shows the benchmark of the X-point JET-like case, when compared to the reduced-MHD model. Poloidal cross sections of n ¼ 20 perturbed quantities are shown for the toroidal magnetic

ARTICLE
scitation.org/journal/php potential A / and the temperature (for the full-MHD model), together with the initial pedestal pressure gradient. The growth rates of the modes are plotted as a function of the toroidal mode number and compared with reduced-MHD. Note that this X-point case is run including the non-static background equilibrium, with stationary flows and Mach-1 boundary conditions, so that this level of agreement with reduced-MHD is excellent.

E. Peeling-ballooning modes in a spherical tokamak
The next benchmark of interest is the MAST plasma discharge #24763. This case was the subject of previous studies using reduced-MHD. 65 It is typically assumed that reduced-MHD is not adequate to study MHD instabilities in spherical tokamaks. The first derivations by Strauss 26 relied on a large inverse aspect-ratio assumption to demonstrate the validity of the reduced-MHD ansatz in tokamaks, but it was later demonstrated that the large inverse aspect-ratio is not needed to ensure energy conservation. 56,57 Nevertheless, this original large inverse aspect-ratio assumption has affected the reputation of the model even until now. Of course, there are instabilities where reduced-MHD fails, like the high-b internal kink mode, but for ELM physics, it seems not to be the case. As such, this benchmark is a validation of reduced-MHD as much as the new full-MHD model.
The MAST discharge #24763 is a nearly symmetrical double Xpoint equilibrium, with the following characteristics and diffusion coefficients: where the D ? and j ? values are given on the magnetic axis, with a radial profile that drops by a factor 10 at the top of pedestal. This MAST benchmark is also run including the complete nonstatic equilibrium with flows and Mach-1 boundary conditions. Figure 6 shows poloidal cross sections of n ¼ 20 perturbed quantities for the toroidal magnetic potential A / and the temperature (for the full-MHD model), together with the initial pedestal pressure gradient. The growth rates of the modes are plotted as a function of the toroidal mode number and compared with reduced-MHD. The agreement between both models shows that not only the new full-MHD model behaves as expected but also reduced-MHD is capable of capturing the linear stability of peeling-ballooning modes for spherical tokamaks. As will be shown later in Sec. IV, the non-linear dynamics of ELM filaments is also very similar between reduced-and full-MHD.

F. Diamagnetic effects and neoclassical friction
The diamagnetic terms were tested using the X-point JET-like benchmark presented above. Figure 7(a) shows the growth rate with and without diamagnetic effects, as a function of the toroidal mode number, for both reduced-and full-MHD. The agreement between the two models is reasonable, showing the expected stabilization of high-n modes.

ARTICLE scitation.org/journal/php
The neoclassical friction terms, as well as the toroidal momentum source, are tested using the same X-point case, to be compared against the reduced-MHD model. Figure 7(b) shows the comparison of velocity profiles in the toroidal and poloidal directions. This case is the same X-point JET-like plasma, which runs as a stationary axisymmetric kinetic equilibrium. For both the toroidal momentum and the neoclassical poloidal velocity, the full-MHD model is able to converge to the target prescribed profile. Note that both S NBI andṽ neo are prescribed only up to the separatrix. The target profile for the neoclassical velocity, as plotted in Fig. 7(b), is determined dynamically by the axisymmetric n ¼ 0 density and temperature profiles, as well as magnetic field, following Eq. (17). The toroidal momentum source S NBI , however, is a fixed profile here, which does not evolve in time.

G. Convergence and resolution
Convergence of growth rates as a function of spatial grid resolution is tested for the tearing and ballooning modes, using the tearing

ARTICLE
scitation.org/journal/php test-case described above and the CBM18 ballooning case. The grid resolution is scanned homogeneously in the radial and poloidal directions, scanning from ðn flux ; n tht Þ ¼ (27,180) to (90, 600), where n flux and n tht are the number of radial (flux surfaces) and poloidal (theta) grid-points, respectively, both of which are equidistant in real-space for this case. With finite elements, the local error is estimated as 37 where h is the element size and p is the polynomial order of the finite elements. In this case, ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi n flux n tht p Þ À1 is used as an approximation of the element size h, and since the elements are cubic, p ¼ 4. However, these tests evaluate the growth rates of toroidal modes, which are obtained by integrating the mode energies over the domain (i.e., over the elements), hence adding another factor ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi n flux n tht p Þ À1 to the error estimate. Thus, the error of the growth rates should scale with the fifth power of the spatial resolution, as ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi n flux n tht p Þ À5 . Although the tearing mode scales as expected, slower convergence is found for ballooning modes, such that the resolution required for full-MHD simulations is significantly higher than for reduced-MHD, with which the fifth power convergence is obtained for the CBM18 case and the tearing mode. Figure 8 shows the convergence of the growth rate error, as a function of spatial resolution. For the tearing mode, the fifth order scaling is found as expected, and beyond a high enough resolution, the error diminishes dramatically, suggesting that the growth rate is already fully converged. However, for ballooning modes using the CBM18 case, the growth rate convergence follows a scaling of the third order. Unfortunately, the reason for this poorer convergence has not been elucidated yet. This lower scaling of ballooning instabilities seems to be similar for all toroidal mode numbers.
In practice, this means that a higher resolution is required for full-MHD simulations than for reduced-MHD to obtain the same error level. It should be noted that this scan has been achieved by scaling the spatial resolution homogeneously everywhere in the plasma, meaning that the same amount of spatial resolution is set for the core and the pedestal. In typical ELM simulations, however, since the mode is localized in the pedestal region, spatial resolution can be diminished in the core region and increased in the pedestal. This enables a reduction of the cost in radial resolution, but the poloidal resolution, however, needs to be elevated everywhere in the plasma, since ballooning modes span both the high-and low-field sides. The benchmark study presented above for ballooning modes was done at the highest resolution possible within the current capabilities of the code. At this level, using the concentration of radial resolution at the pedestal, the error is negligible, below 1% of the growth-rate amplitude.
The converged growth rates, used to evaluate the errors in Fig. 8, are obtained using the highest resolution cases available. For the tearing mode, since saturation is obtained at a relatively low resolution, this is straightforward. For the CBM18 case, since the highest resolution in Fig. 8(b) is not yet saturated, another simulation was used, with non-equidistant resolution in the radial direction. By concentrating radial points around the pedestal, it is, thus, possible to use a lower n flux and therefore increase n tht such that the growth rates are well converged.
It should be noted that, unlike reduced-MHD, the convergence error for full-MHD is in the convenient direction, meaning that if the simulation is under-resolved, the instabilities are stabilized, whereas it is the opposite in reduced-MHD. Thus, artificial instabilities, which are possible in under-resolved reduced-MHD simulations, are avoided in full-MHD. Convergence tests should be performed for high-b cases to ensure that growth rates are not affected in full-MHD. Note that this aspect seems to differ from other well-established tokamak MHD codes such as NIMROD, 66 but there is unfortunately no clear explanation for this aspect of the simulations. The only insight that may be proposed is that the weak formulation method used in JOREK could have stabilizing properties in under-resolved cases, similar to a variational formulation. However, this contradicts the fact that the opposite behavior is found for reduced-MHD in JOREK, which also employs the weak form method, although the reduced-MHD model needs to deal with higher-order derivatives like toroidal vorticity. 16 FIG. 8. The error of the growth rate is plotted as a function of the spatial resolution ffiffiffiffiffiffiffiffiffiffiffiffiffiffi n flux n tht p for (a) the tearing mode and (b) the CBM18 ballooning mode. The error should converge as the 5th power of the spatial resolution, which is the case for the tearing mode but not for ballooning modes. The reason for the lower convergence of ballooning modes is thought to be due to the interaction of fast-waves, which require higher resolution at high b. [Associated dataset available at https://doi.org/10.5281/zenodo.3971940 (Ref. 58).]

IV. NON-LINEAR SIMULATIONS
In this section, we present three non-linear cases to demonstrate the ability of the model to address current open issues of tokamak MHD. These simulations are not meant to provide a validation against experiments but rather to push the model into highly non-linear and unstable regimes, in order to test the numerical stability of the code. The three cases are an MGI-triggered disruption in JET, ELM filament dynamics in MAST, and an ELM crash in JET.

A. Core MHD and disruptions
The first non-linear case addressed here concerns core-MHD, where a disruption is triggered using the massive-gas-injection method in a JET plasma. The JET experiment is chosen to be the pulse previously studied in Refs. 67 and 68, discharge #86887. The simulation is run with toroidal mode numbers n ¼ 1, 2, 3, 4, and 5, using the same diffusion coefficients as in Refs. 67 and 68, with a resistivity of g ¼ 1:5 Â 10 À6 X m and a viscosity of l ¼ 1:5 Â 10 À7 kg m À1 s À1 .
Since the full-MHD model does not yet include neutral physics as in Refs. 67 and 68, the MGI injection is mimicked by a density source S MGI localized at the outer midplane of the plasma. Where the source is localized, a temperature sink is also introduced to account for the ionization of the neutral gas, such that T sink ¼ ÀnS MGI =2, where n ¼ 13:6 eV is the ionization energy per neutral particle. The toroidally localized source is progressively displaced radially inward, starting from w n ¼ 0.9 until w n ¼ 0.6, spanning 4 ms.
Of course, this is an extremely simplified and crude way to address MGI physics, but it should be emphasized again that the purpose of this simulation is not to reproduce the experiment exactly but to trigger a non-linear crash and test the numerical reliability of the model. Future work will be dedicated to the inclusion of neutral and impurity physics as in Refs. 67 and 68 to address disruptions triggered by MGI and shattered-pellet-injection (SPI). However, such additional ingredients are not expected to affect the numerical stability of the model since neutral and impurity physics mainly enters the system of full-MHD as a density source and a temperature sink (i.e., in a similar manner as what was done here), although in more elaborate models, as in Refs. 67 and 68, supplementary terms enter the momentum equation. Figure 9 shows poloidal snapshots of the simulation as a function of time, for the variables T, q, and A / , together with a Poincare plot of the magnetic field. The evolution of the MHD triggered by the localized density source is found to be the same as observations made in Ref. 68. First, the n ¼ 1 mode becomes unstable in the core, later cascading to n ¼ 2 and higher mode numbers. The thermal quench occurs in two phases, first a loss of the outer layer of the plasma temperature through parallel conduction in a stochastic field region and then through convective mixing in the core during the second phase. Thus, the non-linear dynamics of the modes seems to be well reproduced, including the formation of magnetic islands and stochastic regions. Most importantly, the numerical stability of the model is well demonstrated, particularly considering that numerical hyperresistivity and hyper-viscosity are not necessary here, unlike in the reduced-MHD version, where they are critically essential.

B. MAST scrape-off layer filament dynamics
Another demonstrative example of the ability of the new model to handle non-linear dynamics is an ELM instability in the MAST tokamak, for which filaments are observed to travel far into the scrape-off layer (SOL). As in previous reduced-MHD MAST studies, 65,69 the same dynamics is found with full-MHD, with filaments traveling radially at speeds of about 0.5-3 km/s. These simulations were performed using the ELMy H-mode pulse #24763 with the single mode number n ¼ 20, as in Refs. 65 and 69. All diffusion coefficients were set to the same values as in Ref. 69 (except for viscosity, which is scanned from l ¼ 10 À7 to l ¼ 3:3 Â10 À6 kg m À1 s À1 ). One of the singular characteristics of filament dynamics in simulations is that the radial speed and propagation of the filaments strongly depend on the SOL visco-resistive parameters. Particularly, as observed in Ref. 69, at very low SOL viscosity, the filaments are sheared as they cross the separatrix, which stops their radial motion and leads to filaments traveling radially only up to a few cm in the SOL. More experimental behavior is found at high SOL viscosity, where filaments are observed to travel at constant speed until the domain boundary, where Dirichlet boundary conditions stop the filament evolution. Figure 10 shows poloidal snap-shots of the density in the MAST simulation, showing the time-evolution of filaments as they travel radially in the SOL. As filaments propagate radially, parallel transport convects/conducts plasma toward the target regions, resulting in filamentary lobes that are strongly sheared near the X-point due to the magnetic geometry of field lines. Two cases are plotted, using high viscosity (top row) and low viscosity (bottom row). Viscosity is observed to affect the evolution of filaments in the SOL. At lower viscosity, the initial growth rates of the ballooning modes are higher, resulting in a larger initial radial speed of filaments. However, this lower viscosity also results in a strong shearing of the filaments, which aborts the evolution of the filaments just outside the separatrix. At higher viscosity (top-row), the filaments travel at a lower speed but at constant velocity all the way into the far SOL. The same behavior was observed in Ref. 69, so it is concluded that reduced-and full-MHD are in good agreement, both for the linear stability and the non-linear dynamics of SOL filaments in spherical tokamaks like MAST.
Future works on MAST and MAST-U could contribute to the study of ELM burnthrough, 70 core-MHD and disruption studies, and the exploration of viscous models to describe accurately the evolution of SOL filaments. These will require further developments of the model, particularly the inclusion of neutral physics, which is not addressed in this work.

C. JET ELM simulation
Finally, the new model's capability to handle strongly non-linear MHD instabilities is illustrated with an ELM simulation in the JET tokamak, using multiple mode numbers and diamagnetic effects. The resistivity is taken to be a factor 20 above the experimental neoclassical amplitude. For ELM simulations at experimental resistivity, it is now understood that the temporal evolution of the pedestal buildup is needed, leading up to the ELM onset. 39,71 As demonstrated in Refs. 39 and 71, non-linear ELM onsets have an explosive behavior, which leads to larger ELM crashes than a simulation started from a linearly unstable pre-ELM state. Such simulations, however, require long computation times and a detailed setup, which will be addressed with full-MHD in future work.
The JET pulse used for this simulation is a well-studied ELMy H-mode experiment, #83334, which was studied in Refs. 39, 40, Physics of Plasmas ARTICLE scitation.org/journal/php and 72. The simulation is performed using the toroidal mode numbers n ¼ 1, 2, 3, …, 10. The pedestal losses obtained here are 21% for density and 23% for temperature (relative to the pedestal content, integrated in the range w n ¼ [0.85-1.0]). The mode n ¼ 5 begins the crash, until other modes start interacting in the later phase of the ELM. One intriguing observation made from ELM simulations with the full-MHD model is the clear formation of swirling eddies inside the pedestal as filaments leave the plasma, as shown in Fig. 11. Such eddies were never observed in reduced-MHD simulations, suggesting that the restriction of formulatingṽ ? as a function of a scalar electric potential may play a role in filament convection. However, it is not yet clear whether such differences play any role in global aspects of simulations, where quantitative comparisons of integrated measurements (like ELM size and divertor heat-fluxes) would need to be performed for multiple pulses, as in Refs. 39 and 72. Detailed studies of ELM physics using the full-MHD model will be required in the near future to provide extensive validation of the model against experiments. Such studies are now possible, and the non-linear simulations presented here clearly show that the model is numerically sound and ready for such extensive studies, at the same level as previous studies performed with the reduced-MHD model.  9. Simulation of an MGI-triggered disruption in JET. Columns from left to right: electron temperature T e , electron density n e , toroidal magnetic potential A / , and a Poincare plot of the magnetic field. Rows from top to bottom: time lapse of the simulation, at t ¼ 0 ms, t ¼ 2:8 ms, t ¼ 4:2 ms, and t ¼ 5:7 ms. Despite a very simplified MGI model lacking neutral physics, core-MHD is triggered and a thermalquench is obtained, demonstrating that the new full-MHD model is ready for simulations of disruption physics at the same level as the reduced-MHD model in previous studies. 68 [

V. SUMMARY AND FUTURE WORKS
This paper presents the significant advancement of MHD simulations using the JOREK code, with the inclusion of a full-MHD model. The model was described in detail, including extended MHD terms like momentum sources and diamagnetic effects. A linear benchmark of the new model was performed against the previously studied reduced-MHD model, showing satisfactory agreement. Linear tests were performed for core-MHD circular plasmas and ballooning modes in circular and X-point plasmas, including diamagnetic effects. Additionally, a spherical plasma was tested using MAST experimental data, showing that reduced-and full-MHD are in good agreement, even for small aspect-ratio machines like MAST. It is often asserted that reduced-MHD is unable to describe accurately the physics of MHD instabilities in spherical plasmas, but as was shown here, at least for ballooning modes, this is not the case. Finally, non-linear simulations were performed for a disruption in a JET plasma, for the dynamic evolution of SOL filaments during ELMs in the MAST tokamak, and for a non-linear ELM crash in a JET plasma, showing all the key characteristics expected from previous studies with the reduced-MHD model.
Beyond the demonstration that the new full-MHD model is reliable with respect to physics aspects of MHD instabilities, this paper also demonstrates the remarkable numerical stability of the model. Although numerical stabilization is often required to avoid issues resulting from fast-waves inherent to full-MHD, the simulations presented above demonstrate that no such stabilization is required for many of the experimentally relevant MHD instabilities in tokamak plasmas. It may be expected that in some extreme regimes, such stabilization could be required, and work is currently on-going to implement stabilization methods such as Taylor-Galerkin (TG) or variational-multi-scale (VMS) methods. 51,52 As such, the new full-MHD model is now ready for full production and quantitative validation against experiments, at least at the same level as previous studies performed with the reduced-MHD Filaments evolve radially at speeds in the range of 0.5-3 km/s, depending on viscosity and resistivity. As resistivity increases, the resistive ballooning regime leads to larger growth rates and thus higher radial speeds. Viscosity, however, has a more complex effect. Low pedestal viscosity leads to higher growth rates and higher initial filament speed, but low SOL viscosity leads to a strong shearing of filaments just after they cross the separatrix, as shown in the bottom case of (b). The top row of (b) shows a simulation with higher viscosity, where the filaments travel at a lower radial speed but all the way into the far SOL. Plot (a) shows a poloidal snap-shot of the whole plasma, with the zoom-box inside which the slices of plot (b) were obtained. [Associated dataset available at https://doi.org/10.5281/zenodo.3971940 (Ref. 58).]

ARTICLE
scitation.org/journal/php model. Of course, further developments of the model will now be required in the near future to address specific tokamak instabilities such as disruptions, detachment, and ELM burn-through, which requires neutral and impurity physics. However, it is not expected that these additional physics effects will deteriorate the numerical stability of the model since they enter the system as terms that were already included and tested here (mainly density sources and temperature sinks). Another important aspect of future developments will be to address more elaborate boundary conditions, with wall-extended grids as in Ref. 40 and coupling with the free-boundary model STARWALL, as in Refs. 73 and 74.