Ashish Bhole

Boniface Nkonga

José Costa

Guido Huijsmans

Stanislas Pamela

Matthias Hoelzl

Jorek Team

Stabilized bi-cubic Hermite Bézier finite element method with application to Gas-plasma interactions occurring during massive material injection in Tokamaks

Keywords: Stabilized FEM, Variational Multi-scale methods, Shock-capturing stabilization, bi-cubic Hermite Bézier FEM, Magnetohydrodynamics, Full MHD model, Massive material injection in tokamak plasma, Shattered pellet injection

come

Introduction

We are interested in the numerical modeling of strongly magnetized plasmas in tokamaks, particularly the MHD instabilities occurring in the tokamak plasmas. We assume that the compressible MHD equations describe the characteristic time scales of the MHD simulations under concern. Finite element methods (FEM) provide a flexible setup for developing a computational tool for partial differential equations (PDEs). Galerkin FEM codes developed for simulation of MHD instabilities in tokamak plasma include NIMROD [START_REF] Sovinec | Nonlinear magnetohydrodynamics simulation using high-order finite elements[END_REF], M3D-C1 [START_REF] Ferraro | Calculations of two-fluid magnetohydrodynamic axisymmetric steady-states[END_REF] and JOREK [START_REF] Czarny | Bézier surfaces and finite elements for MHD simulations[END_REF][START_REF] Huysmans | ELMs: MHD instabilities at the transport barrier[END_REF]. One of the motivations for choosing the numerical method for MHD is the divergence-free constraint (∇ • b b b = 0) on the magnetic field (b b b). In the physical model we consider, the magnetic vector potential (a a a) formulation is used such that the divergence-free magnetic field constraint is automatically satisfied: ∇ • (∇ ×a a a) = 0. However, the magnetic vector potential formulation leads to higher-order (more than 2) derivatives in the MHD equations (particularly reduced MHD equations). Therefore, a design of Galerkin FEM to provide C 1 -regularity across the finite elements is desired.

JOREK uses mixed Fourier spectral bi-cubic Hermite Bézier FEM that provides C 1 -regularity. It has been extensively used to simulate a range of tokamak plasma applications, including Edge localized modes (ELMs) and their control, disruptions and their control, and vertical disruption events (VDEs) [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF]. A tokamak is an axisymmetric device whose geometry is torus-shaped, and the 3D domain Ω is considered as the tensor product of a 2D poloidal domain Ω ξ with a 1D periodic domain in the toroidal direction ϕ = [0, 2π[. JOREK uses high-order isoparametric bicubic Hermite Bézier FEM in Ω ξ [START_REF] Czarny | Bézier surfaces and finite elements for MHD simulations[END_REF] while Fourier Spectral method in ϕ. Recently, higher-order extensions of Hermite Bézier FEM in Ω ξ have been implemented in JOREK [START_REF] Pamela | A generalised formulation of g-continuous bezier elements applied to non-linear mhd simulations[END_REF]. Historically, the numerical setup of JOREK has been used to solve reduced MHD equations, while the full MHD equations are a relatively recent addition. The full MHD model is expected to provide more complete modeling of plasma flows. It was implemented in JOREK and validated for simple geometries in [START_REF] Haverkort | Implementation of the full viscoresistive magnetohydrodynamic equations in a nonlinear finite element code[END_REF]. Later, it has been extended for the realistic tokamak geometries in [START_REF] Pamela | Extended full-mhd simulation of non-linear instabilities in tokamak plasmas[END_REF]. JOREK provides an excellent setup to perform comparative studies of the numerical solutions of full and reduced MHD models. Such investigations can improve the understanding of a region of validity of the full MHD models. However, full MHD models accommodate fast magnetosonic waves, as opposed to reduced MHD models, and hence they may demand higher resolution. Further, we face severe stability restrictions while using the full MHD model when shocks/discontinuities are involved in the numerical solution.

MHD instabilities play a critical role in magnetic confinement fusion power plants. Largescale instabilities, often called disruptions, are the most critical obstacle in the path toward smooth functioning of a magnetic confinement fusion power plant [START_REF] Hender | Chapter 3: MHD stability, operational limits and disruptions[END_REF][START_REF] Izzo | The role of MHD in 3d aspects of massive gas injection[END_REF]. Unmitigated disruptions can shut down a fusion power plant during which sudden loss of plasma energy is transferred to the tokamak wall damaging plasma-facing components. The disruption mitigation system (DMS) is currently under design for International Thermonuclear Experimental Reactor (ITER) intending to detect disruption and trigger a harmless artificial disruption. ITER DMS injects a massive amount of material into the tokamak plasma to trigger such harmless disruption. Several techniques for massive material injection (MMI) have been developed/are under development: Massive material injection (MGI) [START_REF] Hender | Chapter 3: MHD stability, operational limits and disruptions[END_REF], Pellet injection (PI), Shattered pellet injection (SPI) [START_REF] Combs | Alternative techniques for injecting massive quantities of gas for plasma-disruption mitigation[END_REF][START_REF] Baylor | Shattered pellet injection technology design and characterization for disruption mitigation experiments[END_REF] etc. At present, SPI is considered to be a favorable candidate for The ITER DMS [START_REF] Lehnen | Disruptions in ITER and strategies for their control and mitigation[END_REF]. A typical disruption has distinct phases: thermal quench (TQ), current quench (CQ), and formation of runaway electrons (RE). Each of these phases involves complex physics, and its lack of deep understanding motivates the disruption-related investigations. The simulation codes such as JOREK are under development to predict such complex physics and validate against theory and experiments. The development of a realistic impurity model is also one of the ongoing development to study disruptions numerically [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF]. On the other hand, with the existing models, MMI simulations for the realistic physical parameters are computationally expensive to perform [START_REF] Hoelzl | First predictive simulations for deuterium shattered pellet injection in asdex upgrade[END_REF] due to the lack of a proper numerical stabilization mechanism. The use of the full MHD model further increases the cost of computations. MMI simulations for large devices such as Joint European Torus (JET), ITER, and DEMO are expected to face further severe numerical challenges. The goal of this work is to develop a numerical stabilized mechanism for the underlying Galerkin FEM to reduce the computational cost of MHD simulations.

Numerical schemes for convection-dominated flows must consider the effects of unresolved scales to ensure the stability of the numerical approach. In the context of compressible hydrodynamics, the pioneering work of von Neumann and Richtmyer [START_REF] Neumann | A method for the numerical calculations of hydrodynamical shocks[END_REF], and its 2D extension by Wilkins [START_REF] Wilkins | Calculation of elastic plastic flow[END_REF][START_REF] Wilkins | Use of artificial viscosity in multidimensionnal fluid dynamic calculations[END_REF] the effects of unresolved scales are formulated as artificial viscosity. Godunov [START_REF] Godunov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF] was the first to introduce an explicit evaluation of the sub-scale effects on the resolved scales via Riemann problems. However, these popular formulations are mainly associated with finite volume and discontinuous Galerkin (DG) methods. Continuous Galerkin finite element methods (FEM) provide a flexible numerical setup for computational applications of partial differential equations (PDE)s. These methods give rise to centered approximations of differential operators and can lead to dispersion errors while solving for the convection-dominated flows. Indeed, Galerkin FEMs do not provide mechanisms to control the sub-scales effects on the resolved scales: stabilization. The vast literature exist for stabilized Galerkin FEM such as Taylor-Galerkin methods [START_REF] Donea | Time-accurate solution of advectiondiffusion problems by finite elements[END_REF][START_REF] Ambrosi | A taylor-galerkin method for simulating nonlinear dispersive water waves[END_REF], bubble functions [START_REF] Brezzi | A relationship between stabilized finite element methods and the galerkin method with bubble functions[END_REF][START_REF] Brezzi | Further considerations on residual-free bubbles for advective-diffusive equations[END_REF], Streamline Upwind Petrov-Galerkin (SUPG) methods [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF], Variational Multiscale (VMS) [START_REF] Hughes | Multiscale and Stabilized Methods[END_REF][START_REF] Hughes | Stabilized Methods for Compressible Flows[END_REF] methods etc. The VMS approach provides attractive guidelines for developing stabilized FEMs where numerical stabilization is achieved by an additional contribution to Galerkin's weak formulation. This contribution mimics the effect of the unresolved scales over the resolved scales. In this way, an upwinding process is introduced that leads to numerical diffusion essentially in the flow direction and avoids crosswind diffusion effects in the multidimensional flows. The critical point of developing VMS-based stabilization is the design of a stabilization (intrinsic time-scale) matrix that controls the numerical stabilization added. The order of the accuracy of the underlying Galerkin FEM is preserved.

VMS stabilized FEM gives an excellent approximation to smooth solutions. However, it does not provide a mechanism for BV stability. High-order numerical methods often lack robustness in the vicinity of shocks where the numerical solution may develop Gibbs oscillations and severely affect the numerical method's stability. Shock-capturing stabilization strategies are often used to enforce the total variation stability where the numerical solution develops sharp gradients. The idea is to introduce the artificial viscosity adaptively only in the regions surrounding shocks. Several strategies have been developed for shock capturing stabilization [START_REF] Johnson | On the convergence of shock-capturing streamline diffusion finite element methods for hyperbolic conservation laws[END_REF][START_REF] Tezduyar | Stabilization and shock-capturing parameters in supg formulation of compressible flows[END_REF][START_REF] Bassi | High-order discontinuous galerkin computation of axisymmetric transonic flows in safety relief valves[END_REF][START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF][START_REF] Dzanic | Utilizing time-reversibility for shock capturing in nonlinear hyperbolic conservation laws[END_REF] which are often dependent on the physical problem and mesh topology. In MMI applications, strong and localized mass sources can give rise to hydrodynamic shocks. As MHD equations are nonlinear and can admit discontinuous solutions, the shock-capturing technique is needed to stabilize the underlying numerical method.

This work aims to develop a numerical stabilization strategy for mixed Fourier spectral bicubic Hermite Bézier FEM applied to plasma fluid models based on VMS formulation and shockcapturing technique. Similar approaches have been successfully used in [START_REF] Speleers | Isogeometric analysis with powellsabin splines for advection-diffusion-reaction problems[END_REF][START_REF] Giorgiani | A stabilized powell-sabin finite-element method for the 2d euler equations in supersonic regime[END_REF] where high order Powell-Sabin FEM is stabilized using SUPG and shock-capturing stabilization to solve Euler's equations and in [START_REF] Trelles | Variational multiscale method for nonequilibrium plasma flows[END_REF] where VMS formulation and discontinuity capturing technique is used to stabilize the underlying FEM to solve non-equilibrium plasma with the application to industrial problems. The numerical stabilization for the reduced MHD models in JOREK has been achieved by Taylor-Galerkin stabilization. Sometimes, fourth-order diffusive terms are also used to damp very high wave-number oscillations in the numerical solution. In this work, we develop numerical stabilization techniques for the underlying FEM and demonstrate their use to simulate gas-plasma interactions in the tokamak plasma. The manuscript is organized as follows. Section 2 presents the fluid models for the magnetized plasma extended for MMI. The description of the FEM used in JOREK and the formulation of the stabilized FEM is the subject of Section 3. An application of the stabilized FEM to MMI in tokamak plasma is presented in Section 4. Finally, Section 5 outlines the conclusions and opens future perspectives.

Governing equations

A plasma consists of charged and uncharged particles (atoms, molecules, ions, and electrons) and exhibits a collective behavior due to the interactions among these particles. 'Kinetic plasma theory' describes the collective behavior of constituent charged particles of the plasma. In the kinetic models, the information of a collection of charged particles is expressed in terms of the distribution functions f α of a species α and their evolution is governed by the Boltzmann or Vlasov equations. For plasma applications of interest here, characteristic length and time scales of the kinetic models, Debye length and the electron plasma frequency take the values of order ≈ 10 -5 m and ≈ 10 12 Hz respectively. These scales are very restrictive from the computational application point of view. In this context, the fluid models offer an attractive choice for computational applications with large domains such as ITER. The fluid models are obtained from the Vlasov/Boltzmann equations by taking their first three velocity moments and bringing closure to the fluid equations under the assumption that the distribution function f α is close to the Maxwellian distribution [START_REF] Goedbloed | Principles of Magnetohydrodynamics: With Applications to Laboratory and Astrophysical Plasmas[END_REF]. The resulting fluid models describe a plasma reasonably well when the plasma is dominated by Coulomb collisions [START_REF] Braginskii | Transport processes in a plasma[END_REF]. Unfortunately, the prohibitive cost of the kinetic approach often leads to the use of the fluid approach even for weakly collisional plasmas. Even then the fluid models have been proven to be a powerful tool to study MHD instabilities.

Tokamak plasma mainly consists of ions of hydrogen isotopes (deuterium) and electrons. In MMI experiments neutral gases (such as deuterium, neon, argon, etc.) are injected into tokamak plasma and the atomic processes occurring between gas and plasma give rise to the charged ions at different ionization levels. To study such applications, a physical model needs to be accounted for gas-plasma interactions among the multiple fluid species. The visco-resistive MHD model used in this work treats the 'impurities' as an additional fluid species which is assumed to be convected with the main plasma velocity, independently of its charge state. Furthermore, the impurities are assumed to be at coronal equilibrium, the consequences of which are discussed along with the impurities modeling in [START_REF] Hu | Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection[END_REF] in the context of reduced MHD. The reduced MHD model with impurities transport has been implemented in JOREK and used to study MMI in tokamak plasma [START_REF] Nardon | Progress in understanding disruptions triggered by massive gas injection via 3d non-linear mhd modelling with jorek[END_REF][START_REF] Nardon | Fast plasma dilution in ITER with pure deuterium shattered pellet injection[END_REF][START_REF] Hoelzl | First predictive simulations for deuterium shattered pellet injection in asdex upgrade[END_REF][START_REF] Hu | Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection[END_REF][START_REF] Nardon | Thermal quench and current profile relaxation dynamics in massive-material-injection-triggered tokamak disruptions[END_REF][START_REF] Hu | Collisionalradiative non-equilibrium impurity treatment for JOREK simulations[END_REF]. In this work, we extend the visco-resistive full MHD model implemented in [START_REF] Pamela | Extended full-mhd simulation of non-linear instabilities in tokamak plasmas[END_REF] to include impurity transport.

The full MHD model extended to include transport of the impurity species is written as:

∂a a a ∂t + v × b b b = -η(j j j -j j j ⋆ ) + ∇Φ (1)
ρ ∂v ∂t + ρv • ∇v + ∇p -j j j × b b b = µ∇ 2 v -S S S v (2)
∂ρ ∂t + ∇ • (ρv) = ∇ • (D∇ρ) + S ρ + S ρ f (3) ∂p ∂t + v • ∇p + γp∇ • v + L (∂, n f , n e , t e ) = (γ -1)∇ • (K + n f D κ ) ∇t + S p (4)
∂ρ f ∂t + ∇ • (ρ f v) = ∇ • (D f ∇ρ f ) + S ρ f (5)
These equations govern the evolution of model variables: the magnetic vector potential (a a a), the plasma velocity (v), the total density (ρ), the internal energy (but written in terms of the total pressure (p)) and the partial density (ρ f ) of impurity species at all the levels of ionization.

The above system of equations is a single fluid, single temperature model where the assumption is made that main plasma ions and electrons are at the same temperature. Without this assumption, one would get the single fluid, two-temperature model where Eq. (4) in the above system of PDEs is replaced by the equations for the internal energies of ions and electrons.

∂p e ∂t + v • ∇p e + γp e ∇ • v + L (∂, n f , n e , t e ) = (γ e -1)∇ • (K e + n f D κ ) ∇t e + S pe (6)
∂p i ∂t + v • ∇p i + γp i ∇ • v = (γ i -1)∇ • (K i ∇T i ) + S pi (7)
The ratio of specific heats γ is assumed to be the same for both electron and ion species and is chosen as 5/3 which corresponds to a monoatomic gas. The absence of the permeability of the free space µ 0 from the MHD equations is due to the choice of the normalization [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF]. The total density ρ is the sum of the partial densities for the main tokamak plasma and impurities species.

The potential vector a a a and the velocity v are written into the cylindrical (or toroidal) coordinate system:

a a a = A r e r + A z e z + ψ r e ϕ and v = v r e r + v z e z + v ϕ e ϕ (8)
where r and z are the poloidal coordinates (that define the poloidal plane) and ϕ is the toroidal angle (that denotes the direction perpendicular to the poloidal plane). The vectors (e r , e z , e ϕ ) form an orthonormal basis which from now on will be designated by the cylindrical basis. The toroidal component of the magnetic potential (A ϕ ) is written in terms of the magnetic flux (ψ = rA ϕ ) [START_REF] Haverkort | Implementation of the full viscoresistive magnetohydrodynamic equations in a nonlinear finite element code[END_REF]. The magnetic field (b b b) and the current (j j j) are given by the following relations

b b b = F r e ϕ + ∇ × a a a and j j j = ∇ × b b b (9)
where F is a given function associated with the axisymmetric ideal MHD equilibrium such that the magnetic field remains divergence-free. The induction Eq. ( 1) is written for the evolution of a a a rather than b b b which requires the inversion of the curl operator and the specification of a gauge (Φ). A convenient choice of the Weyl's gauge is made in the induction equation [START_REF] Haverkort | Implementation of the full viscoresistive magnetohydrodynamic equations in a nonlinear finite element code[END_REF] such that:

∇Φ ≡ 0. ( 10
)
The total particle density is denoted by n 0 . The impurities species contains the neutral atoms and the ions at different ionization levels. The total impurity particle density (n f ) is the sum of the particle densities (n z ) with different ionization levels (z):

n f ≡ ρ f m f = nz z=0 n z
where n z denotes the maximum possible levels and m f denotes the effective mass of the impurities species. The quasi-neutrality assumption (the mass of electrons is negligible compared to the mass of ions) allows us to write the mean charge (z f ) for the impurities species as:

z f = nz z=1 z z n z n f (11)
where z z = ze is the charge associated with the ionization level z. The equation of state formulates a link between the plasma temperature (t), pressure, and partial and total densities. The total pressure is the sum of partial pressures. In addition to the quasi-neutrality of the plasma and the assumption that the electrons, ions, and impurities are at the thermal equilibrium, the total pressure in the normalized units is written as:

p = (ρ + α f ρ f ) t with α f = m i m f z f + 1 z i + 1 -1, (12)
where m i and z i denote the mass and the electric charge of the main ion species (ions in the tokamak plasma). To compute the total pressure, the ratio nz n f is required which is given by a probability function of the local electron temperature (t e ) and the particle density of electrons (n e ) under the assumption of coronal equilibrium [START_REF] Hu | Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection[END_REF] as:

n z n f ≃ P (z, n e , t e )
where t e = t/2 and n e ≃ n 0

The function P is obtained from Atomic Data and Analysis Structure (ADAS) database [41]. The variation of P (z, n e , t e ) according to the electron density n e is negligible, such that we can use a constant density n e ≃ n 0 in our applications. Note that, if the injected gas is the hydrogen isotope from which the tokamak plasma is generated, then we have m i = m f and z f = z i . In this case, α f = 0 and the definition of the total pressure reduces to p = ρ t.

We assume that the atomic reactions weakly affect the partial pressures of each species. The coronal equilibrium assumption results in an instantaneous change in the ionization states resulting in difficulties in treating the ionization energy and the corresponding recombination radiation which would not be present in a self-consistently evolving non-equilibrium model. To avoid such artificial recombination radiation, the ionization energy is treated as a potential energy [START_REF] Hu | Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection[END_REF] and the associated contribution L (∂, ρ f , t e ) is added to the equation of the internal energy.

L (∂, n f , n e , t e ) = (γ -1) ∂ ∂t (n f κ f ) + ∇ • (n f κ f v) with κ f ≡ κ f (n e , t e ) (13)
In the case of the two-temperature model, using conservation of the total energy, this potential energy contribution is added to the electron energy equation as can be seen from Eq. ( 6). The function κ f (n e , t e ) denotes the potential energy of atomic ionization and recombination and is evaluated as:

κ f (n e , t e ) = nz z=0 κ z n z n f = nz z=0 κ z P (z, n e , t e ) where κ z = z k=0 κ k+1 k
The potential energy κ f is the energy required to ionize an atom from the neutral state to the ionized state with level z, while κ k+1 k denotes the energy needed to move an electron from the level k up to level k + 1. The database ADAS [41] provides the potential energies κ k+1 k for the species under consideration. The tensor D κ is a symmetric positive definite and is obtained by applying Fick's law to the fluctuations of the convection velocity, about the average plasma velocity (v), associated with n f κ f .

The resistivity (η) and viscosity (µ) can be specified as anisotropic tensors [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF]. However, in this work, a simplistic choice of assuming these parameters as scalars is made. The resistivity is specified with Spitzer-like temperature dependence as well as the effective charge (z eff ) dependence [START_REF] Wesson | Tokamaks. International series of monographs on physics[END_REF]:

η = η 0 f (t, z eff ) with z eff = nz z=1 z 2 z P (z, n e , t e ) nz z=1 z z P (z, n e , t e ) (14)
The scalar viscosity is also implemented with Spitzer-like temperature dependence µ = µ 0 t -3/2 where η 0 and µ 0 are the respective values at the magnetic axis. The tensors D, D f and K (in the case of the two-temperature model K e and K i ) are all symmetric positive definite. Fick's law applied to the fluctuations of velocities about the average velocity (v) introduces the tensors D and D f . Fourier's law of thermal conduction introduces the thermal conductivity tensors K, K e , and K i . Each of these tensors (denoted by X) is always decomposed into the parallel and perpendicular components to the magnetic field as: since in a strongly magnetized plasma, the dynamics parallel and perpendicular to the magnetic field are generally very different. The scalar coefficients X ∥ and X ⊥ control the diffusions parallel and perpendicular to the magnetic field and can be modeled to have temperature dependence. The perpendicular coefficients can also be modeled to mimic the kinetic turbulent transport [START_REF] Pamela | Extended full-mhd simulation of non-linear instabilities in tokamak plasmas[END_REF]. The term j j j ⋆ denotes a current source term. It can be used to preserve the equilibrium current profile approximately throughout the simulation. The current source term can also be used to model a consistently evolving bootstrap current [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF]. The source terms S ρ and S ρ f denote the mass sources of the tokamak plasma and externally injected impurities, respectively. The source term in the momentum equation contains the contribution of the diffusion terms, density sources as well as the momentum source term S m :

S S S v = S m -(S ρ + S ρ f + ∇ • (D∇ρ) ) v
The source terms in the equation of the internal energy contain the Ohmic heating term, thermal energy source (S e ), particle source effects, and radiation losses.

S p = (γ -1) η(j j j -j j j ⋆ ) • j j j + S e -v • S m + v • v 2 S ρ + S ρ f + ∇ • (D∇ρ) -n e n f L rad (t e ) ,
where L rad (t e ) denotes the radiation power function obtained under coronal equilibrium assumption [START_REF] Summers | Radiative power loss from laboratory and astrophysical plasmas. i. power loss from plasmas in steady-state ionisation balance[END_REF][START_REF] Hu | Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection[END_REF].

The system of equations ( 1)-( 5) with closure relation [START_REF] Baylor | Shattered pellet injection technology design and characterization for disruption mitigation experiments[END_REF] can be written in a compact form as: ∂w w w ∂t + L w w w (∂, w w w)s w w w (w w w) = 0 with w w w = ( a a a, ρv, ρ, p

+ (γ -1)ρ f κ f , ρ f ) t (15)
where w w w denotes the vector of physical variables, L w w w denotes the differential operator while s w w w is the vector of the terms at the right-hand side in the equations including sources. For the above system of equations to be well-posed, it must be supplemented by the initial and boundary conditions. In the numerical simulations of tokamak modeling, usually, the aim is to investigate the MHD instabilities around the plasma equilibrium configuration. Therefore, the initial condition is given by the equilibrium configuration of the ideal MHD equations, to which small perturbations are added, and the MHD equations evolved using the numerical approximation. Physics-based boundary conditions are supplied to the numerical solver.

Initial conditions:

The equilibrium for the steady and stationary ideal MHD equations and the axis-symmetric configuration is governed by the famous Grad-Shafranov equation (GSE) [START_REF] Grad | Hydromagnetic equilibria and force-free fields[END_REF][START_REF] Shafranov | On magnetohydrondynamical equilibrium configurations[END_REF] and is written as:

r 2 ∇ • 1 r 2 ∇ψ = -F dF dψ - dp dψ
where F is an axisymmetric equilibrium function of ψ and does not evolve in time. The magnetic field is defined in terms of F as b b b = ∇ × a a a + (F/r) e ϕ such that Gauss's law is satisfied exactly.

The GSE is a nonlinear elliptic PDE and its solution requires the specification of the boundary data and the profiles F (ψ) and p(ψ). Usually, this data is specified analytically for simple cases or can be extracted from 'geqdsk' files [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF]. For the numerical tests considered here we impose Dirichlet boundary conditions on ψ at the boundaries of the computational domain.

Boundary conditions. The boundary conditions can be set in a flexible way such that the decomposition of the boundary for Dirichlet and Neumann conditions can be different for each variable. In most situations, boundaries are identified whether they are aligned to the magnetic flux surfaces (constant ψ surfaces) or not. By default, all the variables are kept constant in time at the boundaries that are aligned to the flux surfaces (Dirichlet boundary conditions). At the boundaries intersecting with the flux surfaces (e.g., in the divertor region or for the multi-block grids extended to the true physical wall [START_REF] Pamela | A wallaligned grid generator for non-linear simulations of mhd instabilities in tokamak plasmas[END_REF]) the sheath boundary conditions are applied as commonly done in divertor physics codes [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF]. The velocity is forced to be equal to the ion sound speed: v • n = ± √ γ t while the temperatures are constrained by the following condition on the heat flux:

ρ v • v 2 + γ ρ t γ -1 v • n -(K + n f D κ ) ∇t • n = γ sh ρ t v • n
where γ sh is the total sheath transmission factor that has typical values of 7-8.

Although the full MHD equations (15) written above have been described from the tokamak modeling point of view, in general, they govern gas-plasma interactions that occur in a wide range of fields such as Astrophysics, industrial plasma application, etc. The numerical approach used to solve these MHD equations is the topic of the next section.

Stabilized Finite Element Method

The full MHD model described in Section 2 is the system of nonlinear PDEs given by Eq. [START_REF] Neumann | A method for the numerical calculations of hydrodynamical shocks[END_REF] in the unknowns w w w(t, x) for t ⊂ R + defined over the spatial domain x ∈ Ω ⊂ R 3 . This system is compactly written in the residual form as: R R R (w w w) = 0 where, R R R (w w w) := ∂w w w ∂t + L w w w (∂, w w w)s w w w (w w w) = 0 (16)

In the numerical approach considered here, the Galerkin FEM is applied to the above system which is re-written in terms of the interpolated variables u u u(t, x):

u u u = (A r , A z , ψ, v r , v z , v ϕ , ρ, t, ρ f ) t
and the corresponding residual form becomes:

R R R (w w w) = M (u u u) R R R u (u u u) = 0 where, R R R u (u u u) := ∂u u u ∂t + L (∂, u u u) -s (u u u) = 0 ( 17
)
where L (∂, u u u) denotes the differential operator of the system given by R R R u (u u u) while s (u u u) denotes the terms at the right-hand side along with the sources. The mass matrix M is given by:

M =       I 0 0 0 0 0 ρI 0 0 0 0 0 1 0 0 0 0 t ρ t f 0 0 0 0 1       (18)
where ρ and t f comes from the definition of the extended pressure which is given by: P = (ρ + α f ρ f ) t + (γ -1)ρ f κ f as follows:

ρ = ρ + α f ρ f + ρ f t ∂α f ∂t + (γ -1)ρ f ∂κ f ∂t t f = α f t + (γ -1)κ f
The initial condition is given by u u u(t = 0, x) = u u u 0 (x) ∀x ∈ Ω. The boundary conditions are expressed as:

u u u(t, x) = u u u D ∀x ∈ Γ D and ∂u u u ∂ n (t, x) = u u u N ∀x ∈ Γ N
where the boundary of the domain is denoted by ∂Ω ≡ Γ D ∪ Γ N with Γ D ∪ Γ N = ∅. The symbols Γ D and Γ N denote the parts of the boundary where Dirichlet and Neumann boundary conditions are to be applied while n denotes a unit vector normal to the boundary ∂Ω. The boundary conditions can be specified differently for different variables and hence the decomposition of ∂Ω can also be different for each variable. The initial and boundary conditions to be supplemented to this system are described in Section 2.

In the context of the Galerkin FEM, the variational or weak formulation of the problem ( 17) is written as: Find the trial functions u u u ∈

V := {u u u|u u u(t, x) ∈ H(Ω) nv } such that for all test functions u u u * ∈ V := {u u u * |u u u * (x) ∈ H(Ω) nv } following holds: w w w * • R R R (w w w) dΩ = u u u * • R R R u (u u u) dΩ = 0 ( 19
)
where R R R (w w w) = M (u u u) R R R u (u u u), the modified test functions are u u u * = M t w w w * and nv denotes the number of variables of the system of PDEs. In Galerkin FEM the test functions (w w w * ) and the trial functions (u u u) are sought from the same space of functions. Such a formulation gives rise to the central approximation of the differential operators and can lead to dispersion errors in the numerical solution for convection-dominated problems. The Galerkin FEM does not provide any inherent mechanism for numerical stabilization that can stabilize the spurious waves due to dispersion errors. Furthermore, the high-order methods can develop high wave-number Gibbs oscillations in the vicinity of shocks and can lead to severe stability problems.

The fluid plasma models of interest in this work contain a weakly hyperbolic operator that governs the convection-driven nonlinear dynamics and any central approximation of such fluid models needs additional treatments to mitigate numerical instabilities. Therefore we need a numerical model to stabilize the underlying Galerkin FEM that can remove the spurious waves from the numerical solution without compromising the accuracy of the numerical solution. We present the numerical stabilization model based on the VMS formulation in subsection 3.1. The applications of interest in this work can also have a presence of shocks in the numerical solution of the fluid plasma models. Therefore, we need a numerical technique that can detect the shocks and add artificial viscosity in their vicinity to stabilize the underlying high-order FEM. We present such a shock-capturing stabilization method in subsection 3.2. The two stabilization models are applied to the bi-cubic Hermite Bézier FEM whose description is presented in subsection 3.3.

Variational multiscale formulation : VMS

In the VMS framework, the weak formulation ( 19) is decomposed into coarse and fine scales. The coarse scales (u u u) are assumed to lie in the finite-dimensional finite element sub-space V h while the fine scales (u u u ′ ) in an infinite dimensional sub-space Z with V = V h ⊕ Z. The numerical solution is assumed to be decomposed as into resolved scales (u u u ∈ V h ) and unresolved scales (u u u ′ ∈ Z) such that u u u = u u u + u u u ′ . The similar decomposition is applied to the tests functions: u u u * = u u u * + u u u ′ * . Using the scale decomposition, the variational weak form given by Eq. ( 19) becomes: Find u u u ∈ V h and u u u ′ ∈ Z such that following holds:

(u u u * + u u u ′ * ) • R R R u (u u u ′ + u u u ′ ) dΩ = 0 ∀u u u * ∈ V h , ∀u u u ′ * ∈ Z (20)
The primary purpose of the VMS formulation is to derive an estimation of the unresolved scales u u u ′ which will be used to model the numerical stabilization for the resolved scales. This process is similar to an up-winding mechanism that stabilizes finite volumes schemes when applied to convection-dominated flows [START_REF] Maire | Multi-scale godunov-type method for cell-centered discrete lagrangian hydrodynamics[END_REF][START_REF] Balsara | Multidimensional hlle riemann solver: Application to euler and magnetohydrodynamic flows[END_REF][START_REF] Balsara | A two-dimensional riemann solver with self-similar sub-structure -alternative formulation based on least squares projection[END_REF][START_REF] Balsara | Multidimensional riemann problem with self-similar internal structure-part iii -a multidimensional analogue of the hlli riemann solver for conservative hyperbolic systems[END_REF][START_REF] Dumbser | A new efficient formulation of the hllem riemann solver for general conservative and non-conservative hyperbolic systems[END_REF]. The variational weak from ( 20) is then split into the problems for coarse and fine scales as:

u u u * • R R R u (u u u + u u u ′ ) dΩ = 0 ∀u u u * ∈ V h (21) u u u ′ * • R R R u (u u u + u u u ′ ) dΩ = 0 ∀u u u ′ * ∈ Z ( 22
)
The residual is approximated using Taylor's series expansion about u u u as:

R R R u (u u u + u u u ′ ) ≈ R R R u (u u u) + L (u u u ′ ) where L ≃ ∂R R R u ∂u u u (u u u) . ( 23
)
where the operator L is the full Jacobian of R R R u (u u u) with respect to u u u. Inverting the above equation gives u u u ′ ≃ -τ τ τ R R R u (u u u) where τ τ τ is an approximation of L -1 . Substituting u u u ′ in Eq. ( 22)

and using the adjoint duality, Eq. ( 21) can be written as:

u u u * • R R R u (u u u) dΩ - L t (u u u * ) • (τ τ τ R R R u (u u u) ) dΩ = 0 ∀u u u * ∈ V h (24)
The first term in the above equation is the classical Galerkin finite element approximation that solves for the resolved scales u u u. The second term is the effect of the unresolved scales (sub-scales) on the resolved scales.

The evaluation of the exact stabilization matrix (L -1 ) is the main ingredient of the VMS stabilized FEM [START_REF] Hughes | Multiscale and Stabilized Methods[END_REF][START_REF] Hughes | Stabilized Methods for Compressible Flows[END_REF]. For linear operators, it is possible to express the stabilization matrix (L -1 ) without any approximation [START_REF] Hughes | Multiscale and Stabilized Methods[END_REF]. However, in general, evaluating the exact form of L -1 poses a more difficult problem than the original one. For nonlinear problems, an optimal way to obtain L -1 is still an open question. Often, the derivation of an approximated stabilization matrix τ τ τ is based upon heuristic arguments and some physical properties of the physical model [START_REF] Billaud | A simple stabilized finite element method for solving two phase compressible-incompressible interface flows[END_REF][START_REF] Trelles | Variational multiscale method for nonequilibrium plasma flows[END_REF][START_REF] Bayona | Solution of low mach number aeroacoustic flows using a variational multi-scale finite element formulation of the compressible navier-stokes equations written in primitive variables[END_REF]. In this work, we use a diagonal stabilization matrix given by:

τ τ τ = h e λ e max diag C a a ar , C a a az , C a a a ϕ , C vr , C vz , C ρ , C t , C ρ f ( 25
)
where h e is a local mesh size and λ e max is the local maximum wave speed of the physical model over a finite element. To estimate h e approximately, we assume that each element e is formed by linear edges and h e is taken as the minimum of lengths of two line segments formed by joining midpoints of the linear edges of e. The local wave speed is estimated as the maximum speed of the magneto-sonic wave over each element. The free parameters C i ≥ 0, ∀i ∈ 1, ..., nv in the above expression are the constant coefficients that control the amount of the numerical stabilization added in each equation of a system of PDEs. The ratio (h e /λ e max ) denotes a timescale and goes to zero as the mesh size h e → 0 since λ e max is independent of the mesh size. We note that if L is symmetric and when R R R u (u u u) is approximated as L(u u u), the resulting stabilization term is dissipative. Then the numerical stabilization model is written as the sum of dissipative contributions:

L t (u u u * ) • (τ τ τ R R R u (u u u) ) dΩ ≃ ℓ L t ℓ (u u u * ) • (τ τ τ ℓ L ℓ (u u u)) dΩ (26)
It is a common practice to use approximations of the residual by excluding the terms of the higher-order derivatives [START_REF] Shakib | A new finite element formulation for computational fluid dynamics: X. the compressible euler and navier-stokes equations[END_REF][START_REF] Bazilevs | Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows[END_REF][START_REF] Trelles | Variational multiscale method for nonequilibrium plasma flows[END_REF] to reduce the computational cost and complexity. We follow here a similar strategy to write the VMS-based numerical stabilization model:

D vms u u u, u u u ⋆ = L t a (u u u * ) • (τ τ τ L a (u u u) ) (27)
where L a (u u u) denotes approximation of the residual such that the second-order derivative terms are excluded. The operator L a for the full MHD equations considered here is written as:

L a =               -v × ∇× 0 0 0 0 0 v • ∇ T ρ ∇ ∇ 0 0 ρ ∇ T v • ∇ 0 0 0 (γ -1) T ∇ T 0 v • ∇ 0 0 0 0 0 v • ∇               (28)
In the further simplified option, the VMS-based numerical stabilization model can written as:

D vms u u u, u u u ⋆ = L t d (u u u * ) • (τ τ τ L d (u u u) ) (29)
where L d denotes the diagonal part of the operator L:

L d =               -v × ∇× 0 0 0 0 0 v • ∇ 0 0 0 0 0 v • ∇ 0 0 0 0 0 v • ∇ 0 0 0 0 0 v • ∇               (30)
Furthermore, it is assumed that the small scales in the above equation vanish at the element boundaries meaning that the operator L (or any of its simplified version) is piece-wise defined on each finite element of a mesh.

Shock-capturing stabilization

Shock-capturing stabilization introduces artificial viscosity only in the vicinity of shocks such that high wave-number Gibbs oscillations near shocks are removed from the numerical solutions and at the same time discontinuities in the variables are captured. Such an adaptive numerical stabilization to be added in the Galerkin weak formulation as an additional contribution is written as follows:

D sc u u u, u u u ⋆ = ∇u u u ⋆ : τ sc τ sc τ sc ∇u u u (31)
where τ sc τ sc τ sc is solution-dependent shock-capturing symmetric positive definite stabilization matrix. In the application of interest here, we primarily deal with the hydrodynamic shock arising due to the MMI in tokamak plasma and hence we write the shock-capturing strategy for hydrodynamic variables only.

In strongly magnetized plasma, the magnetic field introduces a strong anisotropy in plasma flows and for generality, the anisotropic shock-capturing stabilization matrix (τ sc τ sc τ sc ) formulated here is written as:

τ sc τ sc τ sc = τ sc c ∥ b b b ⊗ b b b + c ⊥ (I -b b b ⊗ b b b) where b b b = b b b ∥b b b∥
where τ sc is a solution-dependent scalar function that detects shocks. Such an approach is similar to that used in [START_REF] Tezduyar | Stabilization and shock-capturing parameters in supg formulation of compressible flows[END_REF][START_REF] Bassi | High-order discontinuous galerkin computation of axisymmetric transonic flows in safety relief valves[END_REF][START_REF] Giorgiani | A stabilized powell-sabin finite-element method for the 2d euler equations in supersonic regime[END_REF] where anisotropic directions are mostly pressured gradientbased. Using the anisotropic stabilization matrix in Eq. ( 31) the shock-capturing stabilization method takes the form:

D sc u u u, u u u ⋆ = τ sc c ⊥ ∇u u u ⋆ : ∇u u u + τ sc c ∥ -c ⊥ b b b • ∇u u u ⋆ • (b b b • ∇u u u)
where c ∥ and c ⊥ are the user-defined parameters to control the amount of the shock-capturing stabilization to be added in the parallel and perpendicular direction to the magnetic field. Fixing optimum values of these parameters requires some tuning and in general, their values depend upon the physical parameters and local mesh size. The isotropic shock-capturing stabilization method is recovered if one sets c ∥ = c ⊥ . The function τ sc used here is motivated by the shockcapturing stabilization technique used in [START_REF] Bassi | High-order discontinuous galerkin computation of axisymmetric transonic flows in safety relief valves[END_REF] for the discontinuous Galerkin method applied to 2D Navier-Stokes equations. A similar technique has also been used in [START_REF] Giorgiani | A stabilized powell-sabin finite-element method for the 2d euler equations in supersonic regime[END_REF] for C 1 Galerkin FEM applied to 2D Euler equations. The form of τ sc used in this work is written as:

τ sc = h 2 e (d p + d s ) p f p with f p = h e ||∇p|| p ( 32
)
where p denotes the total pressure of the plasma given by [START_REF] Baylor | Shattered pellet injection technology design and characterization for disruption mitigation experiments[END_REF]. The dimensionless quantity f p acts as a shock detector. The numerical solution dependent multiplier term h e (d p + d s ) denotes the shock-strength and modulates the value of τ sc . Another mesh-dependent multiplier h e implies that stabilization goes to zero as a mesh is refined. The shock-strength terms d p and d s are defined as:

d p = ∂p ∂ρ L(ρ) + ∂p ∂ρ f L(ρ f ) + ∂p ∂t L(t) and d s = ∂p ∂ρ S(ρ) + ∂p ∂ρ f S(ρ f ) + ∂p ∂t S(t) (33)
where the quantities S(ρ), S(ρ f ) and S(t) denote the source terms in the density, impurity density and temperature equations respectively. To this point, we have developed the numerical strategies to stabilize a Galerkin FEM but we have not specified the finite element sub-spaces yet. Hence the formulation of the numerical stabilization remains valid for any Galerkin FEM. We aim to stabilize the mixed Fourier-spectral bi-cubic Hermite Bézier FEM [START_REF] Czarny | Bézier surfaces and finite elements for MHD simulations[END_REF] that has been implemented in the computational framework of JOREK. This FEM has been extensively used to simulate many plasma flow applications to tokamaks using reduced MHD models and most of the simulations use Taylor-Galerkin stabilization [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF][START_REF] Nkonga | VMS Finite Element for MHD and Reduced-MHD in Tokamak Plasmas[END_REF]. The objective of this work is to implement the developed numerical stabilization strategies in JOREK for simulating plasma flows using the reduced and full MHD models. Next, we briefly describe the high-order Galerkin FEM that has been implemented in JOREK, the solution strategy for the time integration, and the weak form of the full MHD model.

Finite element method

High-order continuity is particularly needed when dealing with physical models that include highorder spatial derivatives, for example, the reduced MHD models implemented in JOREK [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF] have third-order derivative terms. In addition, it is desired that unperturbed ideal MHD equilibrium should be accurately preserved for a long time. Moreover, numerical precision is vital in solving nonlinear problems where high-order FEM provides obvious advantages. Isoparametric bi-cubic Hermite Bézier (C 1 ) FEM has been implemented in the poloidal plane while Fourier spectral method in the periodic toroidal direction.

In the Galerkin finite element framework, the vector of the variables (u u u) lies in the finitedimensional finite element sub-space V V V h such that:

w w w * •R R R (u u u) dΩ = 0 ∀w w w * ∈ V V V h where R R R (u u u) ≡ R R R (w w w(u u u)) = R R R t ∂u u u ∂t , u u u +R R R x (u u u) (34)
where R R R t (.) denotes the time derivative terms while R R R x (.) denotes the remaining terms in the strong form of the system [START_REF] Wilkins | Use of artificial viscosity in multidimensionnal fluid dynamic calculations[END_REF]. The weak form of the equations in this system is obtained by applying the following projections:

w w w ⋆ a a ar = ( n i e e e r , 0, 0, 0, 0)

t w w w ⋆ ρ f = (0, 0, 0, 0, n i ) t ( 35
)
where n i ∈ V h denote a scalar test functions. Using these projections the weak form is compactly written in the residual form as:

A t A t A t ∂u u u ∂t , u u u, n i + A x A x A x (u u u, n i ) = 0 ∀n i ∈ V h (36)
with

A t A t A t ∂u u u ∂t , u u u, n i = R R R t • w w w ⋆ a a ar , R R R t • w w w ⋆ a a az , R R R t • w w w ⋆ ψ , R R R t • w w w ⋆ vr , R R R t • w w w ⋆ vz , b b b r R R R t • w w w ⋆ vr + b b b z R R R t • w w w ⋆ vz + b b b ϕ R R R t • w w w ⋆ v ϕ , R R R t • w w w ⋆ ρ , R R R t • w w w ⋆ t , R R R t • w w w ⋆ ρ f t A x A x A x (u u u, n i ) = R R R x • w w w ⋆ a a ar , R R R x • w w w ⋆ a a az , R R R x • w w w ⋆ ψ , R R R x • w w w ⋆ vr , R R R x • w w w ⋆ vz , b b b r R R R x • w w w ⋆ vr + b b b z R R R x • w w w ⋆ vz + b b b ϕ R R R x • w w w ⋆ v ϕ , R R R x • w w w ⋆ ρ , R R R x • w w w ⋆ t , R R R x • w w w ⋆ ρ f t
The projection for the toroidal component of the momentum equation uses a linear combination with the components of the magnetic field. This approach, proposed in [START_REF] Haverkort | Implementation of the full viscoresistive magnetohydrodynamic equations in a nonlinear finite element code[END_REF], removes noisy Lorentz force terms (j j j × b b b) by eliminating fast waves from the toroidal component of the momentum equation. The interpolation of the variable u u u in the finite element space takes the following form:

u u u (t, X, ϕ) = j u u u j (t) n j (X, ϕ)
The index j runs over the number of basis functions N j , u u u j (t) denotes the degrees of freedom (coefficients of the interpolation) and X denotes the poloidal coordinates (r, z) T . For any test function n i ∈ V h , the VMS based numerical stabilization model ( 27) applied to the full MHD equations ( 17) (with u u u * = M t w w w * ) is written as:

D D D vms (u u u, n i ) = D vms u u u, w w w ⋆ a a ar , D vms u u u, w w w ⋆ a a az , D vms u u u, w w w ⋆ ψ , D vms u u u, ρw w w ⋆ vr , D vms u u u, ρw w w ⋆ vz , b b b r D vms u u u, ρw w w ⋆ vr + b b b z D vms u u u, ρw w w ⋆ vz + b b b ϕ D vms u u u, ρw w w ⋆ v ϕ , (37)
D vms u u u, tw w w ⋆ t + D vms u u u, w w w ⋆ ρ , D vms u u u, ρw w w ⋆ t , D vms u u u, t f w w w ⋆ t + D vms u u u, w w w ⋆ ρ f t dΩ
The test functions n i are described in detail below. The shock-capturing stabilization scheme applied to the full MHD equation ( 17) is written as:

D D D sc (u u u, n i ) = 0, 0, 0, D sc u u u, ρw w w ⋆ vr , D sc u u u, ρw w w ⋆ vz , b b b r D sc u u u, ρw w w ⋆ vr + b b b z D sc u u u, ρw w w ⋆ vz + b b b ϕ D sc u u u, ρw w w ⋆ v ϕ , (38)
D sc u u u, tw w w ⋆ t + D sc u u u, w w w ⋆ ρ , D sc u u u, ρw w w ⋆ t , D sc u u u, t f w w w ⋆ t + D sc u u u, w w w ⋆ ρ f t dΩ
The weak form of the full MHD model with the VMS and shock-capturing stabilization terms are compactly written as:

A t A t A t ∂u u u ∂t , u u u, n i + A x A x A x (u u u, n i ) = D D D vms (u u u, n i ) + D D D sc (u u u, n i ) ∀n i ∈ V h (39)
All the operators in the equation above are linear in u u u ⋆ and nonlinear in u u u. The first contribution on the left-hand side comes from the time derivative terms while the second contribution comes from the rest of the terms in the physical model. Contributions on the right-hand side denote the VMS-based and shock-capturing terms respectively. Hence, the implementation of the stabilized FEM consists of the addition of two new contributions to the existing numerical algorithm. It involves computations of the additional terms however all computations are local and no additional communications associated with the parallel programming are involved. Similar strategies have been used for 3D non-equilibrium plasma simulations in [START_REF] Trelles | Variational multiscale method for nonequilibrium plasma flows[END_REF] of some industrially-relevant plasma flows where L 2 stabilization is achieved using VMS formulation and BV stabilization is achieved using the discontinuity capturing method. In [START_REF] Giorgiani | A stabilized powell-sabin finite-element method for the 2d euler equations in supersonic regime[END_REF] stabilized Powell-Sabin FEM is validated for 2D Euler's equations where upwinding is achieved using SUPG and the shock-capturing stabilization uses a similar strategy as described above.

Time integration. The semi-discrete form of the system given by Eq. ( 36) is discretized in time by a general three-level implicit second-order backward difference (also known as Gears) [START_REF] Beam | Alternating direction implicit methods for parabolic equations with a mixed derivative[END_REF] method. The method has been adapted for variable time stepping and is written as:

A t A t A t (1 + r n ζ)u u u n+1 -(1 + ζ(1 + r n ))u u u n + ζu u u n-1 ∆t n , u u u n , n i + A x A x A x θu u u n+1 + (1 -θ)u u u n , n i = D D D vms θu u u n+1 + (1 -θ)u u u n , n i + D D D sc θu u u n+1 + (1 -θ)u u u n , n i (40)
where r n is the ratio of successive time steps ∆t n-1 /∆t n . The choice of parameters ζ = 0 and θ = 1/2 gives Crank-Nicolson method. Gears scheme is given by ζ = 1/2 and θ = 1 while the first order implicit Euler method corresponds to ζ = 0 and θ = 1. When θ > 0, the numerical scheme is implicit, and in each step, we solve a nonlinear system using a Newton procedure. In practice, one Newton iteration with an approximated Jacobian is found to be enough, and the resulting linear system is solved by a preconditioned GMRES method implemented via external libraries (see [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF] for details).

Poloidal and Toroidal decomposition of

V h .
Since tokamak geometry is axisymmetric, the computational domain considered is a tensor product of a 2D poloidal domain and 1D toroidal periodic domain. The finite element space V h is therefore written as the product of the function spaces in the poloidal plane and toroidal direction.

Ω = Ω X ⊗ [0, 2π] =⇒ V h = V X h ⊗ V ϕ h =⇒ n i (X, ϕ) = n X i (X) n ϕ i (ϕ)
Due to the periodicity of the domain in the toroidal direction, a natural choice of the Fourier Spectral method is made where the space V ϕ h is spanned by Fourier modes.

V ϕ h = {1} ∪ span cos (kϕ) , sin (kϕ) , k ∈ ϑ ⊂ N ⋆ and n ϕ i (ϕ) ∈ V ϕ h
where ϑ is a set of positive integers and the tests functions n ϕ i (ϕ) have high-order regularity. With the Fourier method chosen in the toroidal direction, it is up to the design of the space of poloidal basis functions that will give C 1 continuous FEM.

bi-cubic Hermite Bézier FEM

To achieve C 1 -regularity in a 2D poloidal plane, the isoparametric bi-cubic Hermite Bézier FEM with the curved quadrangular elements in the physical space is designed [START_REF] Czarny | Bézier surfaces and finite elements for MHD simulations[END_REF]. The mapping from the reference element τ = [0, 1] × [0, 1] to the physical element (e) is expressed using the third-order Bernstein polynomials basis as:

X = X e (s, t) = 3 i=0 3 j=0 P e ij B i,j (s, t) , X = r z ∈ e, (41)
where P e ij denote the control points for an element e and B i,j are bi-cubic Bézier basis functions. We introduce the following nomenclatures: the reference element τ (s, t), the parametric space ζ ζ ζ ≡ (ζ 1 , ζ 2 ) t and the arc-length space ξ ξ ξ ≡ (ξ 1 , ξ 2 ) t . The mapping among these spaces is defined as:

X(ξ 1 , ξ 2 ) = X (ξ 1 (ζ 1 (s)), ξ 2 (ζ 2 (t))) (42)
In the physical space, a finite element e has curved edges along which ξ 1 and ξ 2 take constant values. The curves defined by the constant values of ξ 1 and ξ 2 are the straight lines in the parametric space ζ ζ ζ which is decomposed into a set of structured trapezoidal elements. A trapezoidal finite element Q e in the parametric space is defined by the boundaries:

ζ e+ 1 (s) = ζ 1,kl + s δζ e k , ζ e- 1 (s) = ζ 1,k+1,l+1 -s δζ e k+1 ζ e+ 2 (t) = ζ 2,kl + t δζ e l , ζ e- 2 (t) = ζ 2,k+1,l+1l -t δζ e l+1
such that the four vertices of a trapezoidal element

Q e are ζ ζ ζ k,l = ζ e+ 1 (0) ζ e+ 2 (0) , ζ ζ ζ k+1,l = ζ e+ 1 (1) ζ e- 2 (1)
, ζ ζ ζ k+1,l+1 = ζ e- 1 (0) ζ e- 2 (0) , ζ ζ ζ k,l+1 = ζ e+ 1 (1) ζ e- 2 (1)
, There exist linear relations between the controls points P e ij , the coordinates of the vertices X v and the derivatives of ∂ ζ1 X, ∂ ζ2 X and ∂ 2 ζ1ζ2 X at the vertices. However, these derivatives in the parametric space do not have any physical meaning. Nevertheless, we can derive linear relations between these derivatives and the derivatives with respect to arc-length coordinates ∂ ξ1 X, ∂ ξ2 X and ∂ 2 ξ1ξ2 X at the vertices. To enforce the C 1 -regularity in the physical space, arclength derivatives are shared at any vertex sharing four elements. These linear relations can be used to re-write the interpolation in Eq. (41) as [START_REF] Czarny | Bézier surfaces and finite elements for MHD simulations[END_REF]:

X = X e (s, t) = 4 v=1 3 o=0 X o v β o,e v H o v (s, t) (43)
where index v and o denote the vertex number of an element and the degrees of freedom associated with that vertex v and H o v (s, t) are bi-cubic Hermite-Bézier basis functions. These functions associated with the first vertex are written as:

X o v = X v β o,e v = 1 H 0 1 (s, t) = H o (s)H o (t) X 1 v = ∂X ∂ξ 1 v β 1,e v = ± dξ 1 dζ 1 dζ e 1 ds v H 1 1 (s, t) = H 1 (s)H o (t) X 2 v = ∂X ∂ξ 2 v β 2,e v = ± dξ 2 dζ 2 dζ e 2 dt v H 2 1 (s, t) = H o (s)H 1 (t) X 3 v = ∂ 2 X ∂ξ 1 ∂ξ 2 v β 3,e v = β 1,e v β 2,e v H 3 1 (s, t) = H 1 (s)H 1 (t) with H o (z) = (1 + 2z)(1 -z) 2 and H 1 (z) = z(1 -z) 2 .
The basis functions for other vertices can be obtained from H o 1 (s, t) using following replacements:

H o 2 (s, t) = H o 1 (1 -s, t) , H o 3 (s, t) = H o 1 (1 -s, 1 -t) , H o 4 (s, t) = H o 1 (s, 1 -t) , ∀o
The scale factors β o,e v associated to a given vertex are different for each element e sharing the vertex v. The signs of β o,e v depend upon the orientation of each element with respect to the reference element. In the isoparametric framework, the interpolation of any scalar function f (X) in terms of bi-cubic Hermite Bézier basis functions is written as:

f (X) = f (X e (s, t) ) = 4 v=1 3 o=0 f o v β o,e v H o v (s, t) (44)

Meshing for bi-cubic Hermite Bézier FEM

As a demonstration, let us consider the example of an annular domain Ω A with the center given by (r 0 , z 0 ) t :

Ω A = X = r 0 + r cos θ z 0 + r sin θ , r 0 ≤ r ≤ r ⋆ , θ 0 ≤ θ ≤ θ ⋆
where the subscripts 0 and ⋆ denote minimum and maximum values respectively. Then the definitions of arc-length and parametric space become:

ξ 1 = r, ξ 2 = rθ, ζ 1 = r, ζ 2 = θ
The vertices of a trapezoidal element in the parametric space are given by: where indices k and l denote the discretization in r and θ with possibly unequal spacing δr k and δθ l respectively. The arc-length derivatives are given by:

ζ e 1 = r 0 + s δr k , ζ e 2 = θ 0 + t δθ l
∂X ∂ξ 1 = cos θ sin θ , ∂X ∂ξ 2 = -sin θ cos θ , ∂ 2 X ∂ξ 1 ∂ξ 2 = 1 r -sin θ cos θ
The scale factors associated with an element e are written as:

β 1,e k,l = δr k , β 1,e k+1,l = -δr k , β 2,e k,l = r k δθ l β 2,e k,l+1 = -r k δθ l
The application of the above definitions to construct a finite element grid in the parametric and physical spaces is shown in Fig. [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF] for the chosen domain in the parametric space (shown in the left frame in Fig. (1)). The bi-cubic Hermite Bézier FEM leads to the curved elements in the physical space as shown in the right frame and is used to generate grids that are isomorphic to polar grids in the poloidal plane. Polar grids aligned to poloidal magnetic flux surfaces are commonly used for MHD simulations using JOREK. These surfaces define the ideal MHD equilibrium configuration and can be described as Morse functions. Topological features of a Morse function are used to decompose the physical domain into several sub-domains connected by the singular points, out of which at least one is polar [START_REF] Guillard | Tokamesh : A software for mesh generation in Tokamaks[END_REF]. This polar sub-domain is discretized using bi-cubic Hermite Bézier formulation to construct a polar grid. Figure (3c) shows an example of such an advanced multi-block grid with two singular points: one magnetic axis and one saddle point. Sophisticated meshing in JOREK uses the bi-cubic Hermite-Bézier formulation for an accurate representation of the vacuum chamber of tokamaks devices [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF] as well. We can see that the polar grids constructed using the bi-cubic Hermite Bézier formulation have a geometric singularity as r 0 → 0 which is a point of concern numerically. The numerical treatment to be applied for such a singularity has been developed in [START_REF] Ashish | Treatment of polar grid singularities in the bi-cubic hermite-bézier approximations: isoparametric finite element framework[END_REF].

Weak formulation of the full MHD model

The discrete weak form of the system of the full MHD equations is written as: Find the trial functions u u u ∈ V nv h such that for all test functions n i ∈ V h following holds:

w w w * • R R R (u u u) dΩ h = 0 ( 45
)
The projections of the scalar equations are straightforward. For the induction equations, the projections are obtained using the cylindrical basis. For the momentum equation, the projections are obtained using e R , e Z , and b b b. The discrete weak form of each component equation in the full MHD system is written below:

∂a a a ∂t -v × b b b • e r n i = - b b b • ∇ × (ηe r n i ) -e r • j j j ⋆ n i + b a a a R ∂a a a ∂t -v × b b b • e z n i = - b b b • ∇ × (ηe z n i ) -e z • j j j ⋆ n i + b a a a Z ∂a a a ∂t -v × b b b • e ϕ rn i = - b b b • ∇ × (ηre ϕ n i ) -re ϕ • j j j ⋆ n i + b ψ (ρD t v) -∇ (p + π) + b b b • ∇b b b • e r n i = - ∇v : ∇ (µe r n i ) + S S S v • e r n i + b v R (ρD t v) -∇ (p + π) + b b b • ∇b b b • e z n i = - ∇v : ∇ (µe z n i ) + S S S v • e z n i + b v Z ρ (D t v) -∇p • b b bn i = - ∇v : ∇ (µb b bn i ) + S S S v • b b bn i + b v ϕ ∂ρ ∂t + ∇ • (ρv) n i = -((Dρ) • ∇n i -S ρ n i + -S ρ f n i ) + b ρ ∂p ∂t + ∇ • (vp) + (γ -1) p∇ • v n i = - (K f ∇t) • ∇n i -S p n i + b p ∂ρ f ∂t + ∇ • (ρ f v) n i = -((D f ∇ρ f ) • ∇n i -S ρ f n i ) + b ρ f
For conciseness, we have used the notation for the material derivative D t v = ∂v ∂t +v •∇v, omitted the volume integrals dΩ, denoted by b i boundary integral terms and used the following notations in the above equations:

π = b b b • b b b 2 , p = (ρ + α f ρ f ) t + (γ -1)ρ f κ f and K f = (γ -1) (K + n f D κ )
4 Application to SPI

The developed stabilized bi-cubic Hermite Bézier FEM is applied to the simulation of MMI in the JET-like plasma. In context of the reduced MHD models, the neutrals modeling have been developed in JOREK and applied in [START_REF] Fil | Modeling of massive gas injection triggered disruptions in tokamak plasmas[END_REF][START_REF] Nardon | Progress in understanding disruptions triggered by massive gas injection via 3d non-linear mhd modelling with jorek[END_REF][START_REF] Nardon | Fast plasma dilution in ITER with pure deuterium shattered pellet injection[END_REF][START_REF] Hoelzl | First predictive simulations for deuterium shattered pellet injection in asdex upgrade[END_REF][START_REF] Nardon | Thermal quench and current profile relaxation dynamics in massive-material-injection-triggered tokamak disruptions[END_REF] and the impurities modeling in [START_REF] Hu | Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection[END_REF]. Similar studies for MMI in tokamak plasma have also been performed using M3D-C1 in [START_REF] Ferraro | 3d two-temperature magnetohydrodynamic modeling of fast thermal quenches due to injected impurities in tokamaks[END_REF] and using NIMROD in [START_REF] Kim | Shattered pellet injection simulations with nimrod[END_REF]. The physics involved in such applications are highly complex and include distinct phases such as TQ, CQ, production of REs, etc. Numerical simulation of such physics for realistic values of the physical parameters is a challenging and computationally expensive task due to severe stability restrictions on the numerical methods. In [START_REF] Hoelzl | First predictive simulations for deuterium shattered pellet injection in asdex upgrade[END_REF], to save computational time, most numerical simulations presented for SPI in ASDEX Upgrade use parameters higher than their realistic values. Furthermore, the plasma volume for ASDEX Upgrade, JET, and ITER respectively are about 13, 100, and 840 m 3 . Given the size of computational domains of present-day and future tokamaks, it is desirable to develop stable, computationally efficient, and robust simulation tools.

In subsections below, we present numerical simulations for SPI into JET-like plasma using stabilized bi-cubic Hermite Bézier FEM where the plasma is modeled with the reduced or full MHD model (see Section (3)). The physical setup of the problem is similar to the JET-like simulation presented in [START_REF] Hu | Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection[END_REF] and is described briefly below.

Problem Description : JET pulse No. 85943

We consider a JET-like plasma that resembles JET pulse No. 85943 at time 62.4 s with the toroidal magnetic field B ϕ ≈ 3 T, total plasma current I p ≈ 2 MA, core electron temperature t e (ψ N = 0) ≈ 3. The initial condition is obtained by numerically solving the GSE using bi-cubic Hermite Bézier FEM on a polar grid. The data required to solve GSE: F F ′ , density and temperature is shown in Figure (2a) as the profiles with respect to the normalized magnetic flux (ψ) where ψ = ψ-ψ0 ψe-ψ0 with ψ 0 and ψ e as the values of ψ at the magnetic axis and the plasma edge respectively. The boundary conditions on ψ to be specified on the polar grid are shown in Figure (2b). Contours of the numerical solution of GSE (ψ eq ) are shown in Figure (3a) on the top of the polar grid used. This ψ-field has one minimum which is the location of the magnetic axis and one saddle point, also known as X-point. For computational efficiency purposes, JOREK uses a strategy to construct multi-block grids whose curves are aligned to constant ψ-surfaces on which the MHD equations are evolved. Such a multi-block grid used for simulations presented here is shown in Figure (3b). The figure also shows ψ contours at which the safety factor q takes rational values q = 1 to 7 (marked radially outwards). q = 1 surface lies close to the magnetic axis, while q = 6, 7 surfaces are close to each other and the plasma edge. In MHD instability theory, the rational q surfaces are probable locations at which MHD instabilities occur. The solution (ψ eq ) of GSE also forms the initial conditions for the MHD problem:

a a a(X, 0) = ψ eq (X) R + ϵ êϕ , v(X, 0) = 0, ρ(X, 0) = ρ(ψ eq (X)), t(X, 0) = t(ψ eq (X)) (46)
We consider Dirichlet boundary conditions except at parts of the boundaries that are not aligned to ψ eq , where physics-based Mach-1 and Sheath boundary conditions are applied to the velocity and temperature respectively (see Section 2). The aim is to simulate plasma response to neon or argon SPI. The sketch in Fig. (3c) shows a cone (projected in the poloidal plane) in which SPI fragments exist in the computational domain. The cone's apex (O) denotes the location where a pellet is shattered to form SPI fragments. The point P denotes the location around which the fragments are assumed to appear in the computational domain at time t = t spi . Simulations presented here assume that a solid pellet of neon or argon gas with the total number of atoms N a is shattered into N f number of spherical fragments. The initial velocities V V V ⋆ i and positions X X X ⋆ i of each fragment are chosen randomly about a fixed velocity V V V spi and position vector --→ OP = X X X spi , such that:

∥V V V ⋆ i ∥ = ∥V V V spi ∥ + ν i (δv) spi 2 and X X X ⋆ i = X X X spi + ℓ spi V V V ⋆ i ∥V V V spi ∥ where ν i denotes a random number in the interval [-1, 1]. Magnitudes of fragment velocities V V V * i are distributed randomly around the magnitude of V V V spi in the interval [-(δv) spi /2, (δv) spi /2].
The directions of V V V * i are spread randomly around that of V V V spi such that the angle between V V V * i and V V V spi does not exceed (δα) spi radians. The vector X X X spi is given in the cylindrical coordinates and V V V spi is defined in the Cartesian coordinates as:

X X X spi = r x spi e e e r (ϕ x spi ) + z x spi e e e z + ϕ x spi e e e ϕ (ϕ x spi ) and V V V spi = V x spi e e e x + V y spi e e e y + V z spi e e e z For all the simulations presented here, we use N f = 100 and specify that the fragments are injected into the plasma at the time t spi = 0.03564 ms. Furthermore, we specify X X X spi and V V V spi using following parameters:

r x spi = 3.05, z x spi = 1.7, ϕ x spi = 4.51
V x spi = -89.44, V y spi = -178.885, V z spi = 0 ℓ spi = 0.15, (δv) spi = 200 and (δα) spi = 0.349/2

The initial size of each fragment follows the statistical fragmentation model (see [START_REF] Hu | Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection[END_REF]). The evolution of the position X X X i (t) and size r i ≡ r i (t) of each fragments is governed by the following set of equations [START_REF] Parks | Effect of transonic flow in the ablation cloud on the lifetime of a solid hydrogen pellet in a plasma[END_REF][START_REF] Sergeev | Studies of the impurity pellet ablation in the high-temperature plasma of magnetic confinement devices[END_REF][START_REF] Parks | The ablation rate of light-element pellets with a kinetic treatment for penetration of plasma electrons through the ablation cloud[END_REF]:

dX X X i dt = V V V i and d dt r 5 3 i = -R i n 1 3 e t 5 3 e with X X X i (t ⋆ i ) = X X X ⋆ i , r i (t ⋆ i ) = r ⋆ i ( 47
)
where R i denotes the ablation coefficient and t e ≡ t e (t, X X X i ) and n e ≡ n e (t, X X X i ) denotes the electron temperature and particle density respectively. The injection velocity is assumed to be constant in time and hence the trajectories of SPI fragments are straight lines. The ablation process begins once the injected fragments interact with the main tokamak plasma. Since we assume spherical fragments, the infinitesimal variation of the ablated volume of each fragment is given by -4πr 2 i dr i . Therefore, 'the ablation rate' is given by:

ρf,i (t) = -4πr 2 i dr i dt ρ ⋆ i = 12πR i ρ ⋆ i 5 r 4 3 i n 1 3 e t 5 3
e where ρ ⋆ i denotes the density of the solid pellet. The source term in the impurity density equation is specified as a Gaussian function around the position of each fragment (X X X i ). As a result, the impurities deposited into the main plasma are given by the sum of all fragments:

S ρ f (t, X X X ) = N f i=1 ρf,i (t) ω i exp - (r -r i ) 2 + (z -z 2 i ) δℓ 2 i exp - (ϕ -ϕ i ) 2 δϕ 2 i ( 48
)
where δℓ i and δϕ i denote the parameters to control the shape of a Gaussian cloud. The spreading volume ω i is integral in the entire domain of the spatial Gaussian profiles of the particle. The Runge-Kutta approximation of the ODE [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF] gives the position and the radius of each particle at a given time during the simulation.

Numerical results for SPI in JET-like plasma

In this subsection, we present the application of stabilized FEM to simulate SPI in JET-like plasma. First, we take a simplified situation where the SPI source is assumed axisymmetric and demonstrate the use of the shock-capturing method for reduced as well as full MHD modeling of the plasma. Later, we mainly focus on the full MHD model to demonstrate the use of numerical stabilization to simulate 3D SPI sources in the plasma described by somewhat realistic values of the parameters. The current source j j j ⋆ is removed from these simulations. The numerical results obtained are discussed with the help of some visualizations and integrated quantities.

Computations of the integrated quantities such as thermal energies, plasma current, and magnetic energies are shown in Appendix A.

Axisymmetric sources

The axisymmetric source term is modeled by specifying a high value δϕ i = 1000 in Eq. ( 48) that flattens the shape of the source term in its toroidal extent. While poloidal shape of the source term is specified by δℓ i = 0.08. The resistivity is modeled with the dependence t e and z eff (Eq. ( 14)) and is specified by the value at the magnetic axis η 0 = 4.22 × 10 -8 Ohm.m. The current source j j j ⋆ is removed from these simulations. The simulations are performed with the different values of the isotropic shock-capturing stabilization coefficients to discuss the effect of the numerical stabilization. Specifying coefficients c ∥ = c ⊥ = c implies the use of isotropic shockcapturing stabilization while their values are used to control the amount of the stabilization.

Reduced MHD

Axisymmetric SPI source is introduced in the JET-like plasma modeled by the single temperature reduced MHD model [START_REF] Hu | Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection[END_REF]. The source term is assumed to be formed by shattering a argon pellet, with total number of N a = 1.5 × 10 22 atoms, into N f = 100. The temperature-dependent coefficient of viscosity is specified with the value at the magnetic axis µ = 1.18 × 10 -8 kg.m -1 .s -1 . The perpendicular and parallel particle diffusivities are set to 13.45 × 10 -5 and 6.7264 m 2 .s respectively. The perpendicular and parallel heat conductivity is set to 4.72 × 10 -7 and 946.44 kg.m -1 .s -1 respectively. The simulations use the multi-block grid similar to that shown in Fig. (3b) with 4130 elements in the poloidal plane and toroidal harmonics 0, 1, 2, 3. BDF2 time integration method is used with ∆t = 0.00148 ms and the simulations are performed using isotropic shock-capturing stabilization terms with the coefficient c ∥ = c ⊥ = c in all equations except the induction equation. Figure (4a) and (4b) shows the evolution of the ablation rates and thermal energies respectively in the reduced MHD simulations. Without shock-capturing terms (denoted by c ∥ = c ⊥ = c = 0) the linear solver for finite element matrix fails to converge around 0.59 ms. Indeed, the development of discontinuities brings high wave-numbers into the numerical solution that can not be resolved with the grid used. Stabilized simulations however go further in time as discontinuities are smeared. The difference between the results obtained with non-stabilized and stabilized simulation with c ∥ = c ⊥ = c = 0.1 is insignificant as can be seen in figures (4a) and (4b). With the use of higher numerical stabilization (c ∥ = c ⊥ = c = 1) simulation runs for the longer times and the numerical solution starts to deviate from that in the non-stabilized simulations.

Figure (5a) shows the snapshot of the total particle density (n 0 ) in [10 20 ] plotted on ϕ = 0 plane and at t ≈ 0.562 ms. The straight line shown in black color joins two points (2.02, -0.56) and (3.02, 1.67) in the poloidal plane and passes through the location of the impurity source at t ≈ 0.562 ms. Trajectories of SPI fragments are spread about this line. Figure (5b) shows the comparison of the pressure plotted along this straight line in ϕ = 4.51 plane obtained with the simulations using different values of shock-capturing stabilization coefficients. Peaks in pressure lineout plots are due to the impurity sources at that location. As the impurity source travels inside the tokamak plasma a discontinuity in the pressure is developed (see Fig ( 5b)). The use of the shock-capturing terms adds the numerical stabilization locally near the shock to smear the discontinuity while everywhere else does not significantly affect the numerical solution. Figure (5b) shows that an increase in the amount of stabilization further smears the discontinuity.

Note that the simulations presented in this subsection use a low-resolution grid. The purpose of this test was to demonstrate that the shock-capturing stabilization terms act locally to improve the stability of the underlying FEM. At the same time, the test also points out to difficulties in choosing optimal values of the stabilization coefficients c ∥ = c ⊥ = c. Analyzing such non-linear numerical methods is often very difficult and fixing optimum values of the stabilization coefficients is based on experience. Furthermore, in the applications of interest here, these coefficients depend upon the physical parameters such as η, µ, D, K as well as the grid resolution.

Full MHD

Here, an axisymmetric SPI source is introduced in the JET-like plasma modeled by the single temperature full MHD model (Section 2). The source term is assumed to be formed by shattering a neon pellet, with total number of atoms N a = 1.5 × 10 21 atoms, into N f = 100. The temperature-dependent coefficient of viscosity is specified with the value at the magnetic axis µ = 2.36 × 10 -7 kg.m -1 .s -1 . The perpendicular particle diffusions are set to 67.264×10 -5 m 2 .s while parallel particle diffusions are set to zero. The simulations use the multiblock grid shown in Fig ( 3b) with 7645 elements in the poloidal plane and toroidal harmonics 0, 1, 2, 3 while BDF2 method for time integration with ∆t = 0.000297 ms. The simulations are run using VMS stabilization terms based on the diagonal part of the operator L (Eq. ( 29)) with the coefficient C i = 10 -2 in all the equations. Additionally, the shock-capturing stabilization terms with included with the coefficients c ∥ = c ⊥ = c in all equations except the induction equation. Figure (6a) and (6b) shows comparison of the evolution of ablation rates and thermal energies respectively in the non-stabilized and stabilized full MHD simulations under discussion here. Without shock-capturing terms (denoted by c ∥ = c ⊥ = 0) the linear solver fails to converge around 1.4 ms. Use of the shock-capturing stabilization with c ∥ = c ⊥ = c = 1 allows simulation to run beyond 1.4 ms. In the stabilized simulations the ablation rates are underestimated slightly and show variations after 1.4 ms. Such variations are due to increased MHD dynamics that nonstabilized simulation fails to capture with the prescribed resolution. The evolution of the thermal energies in the two simulations looks almost identical. A decrease in thermal energies coincides with sudden variations in the ablation rates. The stabilized simulation continues to run beyond 2.5 ms.

Figure (7a) shows the snapshot of the total particle density (ρ 0 ) in [10 20 ] plotted on ϕ = 0 plane and at t ≈ 0.095 ms. The straight line shown in black color joins two points (2.02, -0.56) and (3.02, 1.67) in the poloidal plane and passes through the location of the impurity source. Trajectories of SPI fragments are spread about this line. Figure (7b) shows the comparison of the pressure plotted along this straight line in ϕ = 4.51 plane obtained with the simulations using different values of shock-capturing stabilization coefficients. Peaks in pressure lineout plots are due to the impurity sources at that location. As the impurity source travels inside the tokamak plasma a discontinuity is developed in the pressure field (see Fig ( 7b). Once again, the shock-capturing terms add the numerical stabilization locally near the discontinuity to smear it while everywhere else does not significantly affect the numerical solution. The purpose of this test is to demonstrate that the shock-capturing stabilization terms act locally to improve the stability of the underlying FEM. With an appropriate choice of the stabilization coefficient, it may be possible to run well-resolved simulations with improved stability where the ablation rates do not differ much from the non-stabilized simulation.

3D sources

Simulations presented in this subsection assume that a pellet of neon is shattered into N f = 100 fragments which appear at t SPI = 0.03564 ms in the computational domain. Localized 3D impurity source is used by specifying δℓ i = 0.08 and δϕ i = 2.0. The constant coefficient of viscosity is set to µ = 2.36 × 10 -7 kg m -1 s -1 . The perpendicular particle diffusions are set to 67.264 × 10 -5 m 2 .s while parallel particle diffusions are set to zero. The perpendicular and parallel heat conductivity is set to 1.18 × 10 -6 and 946.44 kg.m -1 .s -1 respectively. Simulations use the multi-block grid as shown in Fig ( 3b) with 7645 elements in the poloidal plane and the toroidal harmonics 0, ..., 7 while BDF2 method for time integration.

The stabilized FEM is used to simulate 3D SPI sources in tokamak plasma modeled with the single temperature full MHD model presented in Section 2 and a comparison with the nonstabilized simulation (without any numerical stabilization) is discussed. The stabilized simulations use the VMS-based stabilization based on the diagonal part of the operator L (Eq. ( 29)) with C i = 10 -2 in all the equations and additionally the shock-capturing stabilization in all the equations except in the induction equation with coefficients c ∥ = c ⊥ = c. The numerical results are discussed below with a focus on shock-capturing stabilization as it is found to be effective in improving the stability of the underlying FEM.

Effect of numerical stabilization

For the full MHD simulations performed with and without numerical stabilization, the resistivity is specified with the dependence on t e and z eff (Eq. ( 14)) with the value at the magnetic axis η 0 = 4.22 × 10 -8 Ohm.m. Three stabilized simulations differ by the values of the isotropic shock-capturing stabilization coefficient equal to 1, 10, and 50. Figure [START_REF] Pamela | Extended full-mhd simulation of non-linear instabilities in tokamak plasmas[END_REF] shows the plasma equilibrium configuration at the time when the impurity source just appears in the computational domain. Figure (8a) shows the initial temperature field plotted in the poloidal plane and the neon SPI fragments (shown by black dots) near the boundary of the computational domain. The fragments are distributed such that initially some of them are even outside the computational domain. As the fragments are far away from the tokamak plasma core region, their effect is not yet felt by the plasma core and the temperature field remains axisymmetric (see Figure (8b)). During the simulations, SPI fragments travel inside the tokamak plasma core with the prescribed constant velocity and neutral atoms of neon from SPI fragments begin to ablate. Modeling of the atomic reactions (ionization, recombination, and radiation) between ablated atoms and the tokamak plasma is included in the MHD equations. Such gas-plasma interactions lead to complex physical phenomena that are challenging to simulate numerically. Moreover, they demand restrictive time steps due to the development of shocks as discussed in the subsections above.

Figure (9a) shows the number of time steps taken by the different simulations vs t in ms where changes in the slopes denote the reduction in time steps. For the sake of presentation, the time steps are marked in the normalized units in Fig. (9a). During simulations, the development of discontinuities and increased MHD dynamics introduce high wavenumbers in the solutions which restrict the stability and leads to failure of convergence of the linear solver. Typically, a simulation is restarted from the time instants before such failures by reducing ∆t. Such reduction of ∆t brings the numerical method to a stable region and allows longer simulations. However, without any numerical stabilization mechanism, ∆t required is practically very small which makes a simulation computationally very expensive. To reach the physical time of almost 1.7 ms the non-stabilized simulation (c = 0) takes more than 9000-time steps and demands further smaller time steps to advance. This simulation was not pursued further to save computational time. The use of numerical stabilization, in particular shock-capturing stabilization, is seen to improve the stability of the method and allows the use of larger ∆t. To reach the physical time of almost 1.7 ms, stabilized simulations take less than 5000-time steps, depending upon the coefficient (c). The simulation with c = 1 takes more than 20000-time steps to reach physical time t ≈ 3 ms. With c = 10 however, the simulation reaches almost the same physical time in less than 6000time steps. A further advantage in terms of time steps is seen with c = 50. The ability to use larger time steps in the stabilized simulations allows us to capture complex dynamics beyond 2 ms in less computational time. Figure (9b) and (9c) shows the evolution of the ablation rates and total number of neon atoms ablated into the plasma. It can be seen that the shock-capturing terms do not significantly affect the ablation process up to 2 ms. After 2 ms large variations in the ablation rates, which are associated with the increased nonlinear MHD dynamics, are seen. These variations are associated with the rapid deposition of impurities (neon atoms) in plasma (Fig. (9c)). Figure [START_REF] Izzo | The role of MHD in 3d aspects of massive gas injection[END_REF] shows the evolution of the magnetic energies in the different toroidal modes in all simulations. The non-stabilized simulation does not proceed beyond t ≈ 1.75 ms with ∆t * = 0.5 normalized units and shows the sign of the failure at the end where the magnetic energies increase suddenly. In this simulation, at initial time instants, the energy cascading shows unphysical behavior (Fig. (10a)) where magnetic energies in mode n = 7 dominate the energies in the lower modes. Such cascading denotes a lack of resolution where higher modes cannot transfer their energies to lower modes. Up to the same resolution, one needs to take into account the effect of the unresolved scales on the resolved scales via a numerical stabilization mechanism. In the stabilized simulations the energy cascading becomes acceptable where the lower modes are seen to dominate the higher modes (Fig. (10b),(10c) and (10d)). Beyond t = 2 ms a sharp increase in magnetic energies is seen which is associated with the increase in the MHD activities. During this highly nonlinear phase, further reductions in the time steps are needed for simulations to be within the numerical stability limits. All stabilized simulations show a sharp rise of the magnetic energies in all modes which is associated with SPI fragments crossing rational q = 2 surface (marked in Fig. Top row in Fig. [START_REF] Combs | Alternative techniques for injecting massive quantities of gas for plasma-disruption mitigation[END_REF] shows snapshots of n f [10 20 m -3 ] plotted in the poloidal plane ϕ = 0 at t ≈ 1.6 ms and obtained using different amount of the numerical stabilization. The difference between the numerical solution obtained without numerical stabilization and with c ∥ = c ⊥ = c = 1 is not significant. Increase in the amount of the shock-capturing stabilization however increases the difference from the non-stabilized results. The particle densities are underestimated in the simulations with c = 10 and 50. Bottom row in Fig. [START_REF] Combs | Alternative techniques for injecting massive quantities of gas for plasma-disruption mitigation[END_REF] shows n f snapshots at t ≈ 2.73 ms in a highly non-linear phase of MHD dynamics. The non-stabilized simulation has not been run until this physical time. Visualizations of the numerical solution obtained with c = 1 and 10 look alike while the dynamics with c = 50 look different. The shock-capturing stabilization is the numerical solution dependent and hence a non-linear technique to achieve stability. When used more than required, it can be sensitive to nonlinear dynamics. Visualization for the simulation with c = 50, in comparison to those with c = 1 and 10, shows that the stabilization has altered the dynamics and is more than required. Below we discuss the numerical results obtained with c = 10 to highlight the complex dynamics of thermal quench occurring in disruptions. Figure [START_REF] Baylor | Shattered pellet injection technology design and characterization for disruption mitigation experiments[END_REF] shows 3D visualizations at some time instants from the simulation with the shockcapturing stabilization coefficient c ∥ = c ⊥ = c = 10. The color map denotes the magnitude of t in keV while iso-surfaces marked by green and red colors denote the impurity particle density (n f ) 35 % and 70 % of its maximum value respectively. In Fig. (12a) the impurity density cloud at t ≈ 0.17 ms is seen to be highly localized in the poloidal plane and does not affect the tokamak plasma core significantly. At this time t-field is almost axisymmetric. At t ≈ 2.31 ms the impurity density cloud has moved inside the tokamak plasma core (Fig. (12b)) and has expanded as a result of its transport modeling. The plasma temperature is seen to drop up to locations of SPI fragments and the temperature field is no more axisymmetric. At t ≈ 2.61 ms the impurity density is seen to have further spread in the plasma core and the temperature field has developed small-scale structures (Fig. (12c)). This rapid spread is due to the crossing of fragments on the q = 2 surface which triggers MHD activities and strong convection carries the impurities into the plasma core. Figure (12d) shows that the impurities have spread almost in the entire plasma core further cooling the plasma. These 3D visualizations show an overview of the complex processes that occur in the gas-plasma interactions.

Figure [START_REF] Lehnen | Disruptions in ITER and strategies for their control and mitigation[END_REF] shows the snapshots of the impurity particle density (n f ) in 10 20 m -3 (in first row) and Poincaré plots of the magnetic field (in second row) plotted at the poloidal plane ϕ = 0 and some time instants. These visualizations are from the simulation with the shock-capturing stabilization coefficient c ∥ = c ⊥ = c = 10. Positions of the SPI fragments are marked using black dots in the Poincaré plots. At t = 2.18 ms the n f -field shows that the impurities are being convected in the domain. The corresponding Poincaré plot of the magnetic field suggests that the plasma temperature has dropped down in the outer region of the tokamak plasma where ergodization of the plasma has occurred. Such ergodization denotes the loss of magnetic confinement. In the outer layer of the tokamak plasma, mode structures associated with m = 3, n = 1 and m = 2, n = 1 are seen which denotes the growth in the corresponding modes. At t ≈ 2.25 ms n f -field begins to show small scales structures where m = 2, n = 1 mode is seen to grow. Further, at t ≈ 2.89 ms the impurities are seen to spread in most of the domain showing many small-scale structures. This is due to the strong convection triggered by the growth of m = 2, n = 1 mode as fragments reach q = 2 surface. During the fast dynamics over which neon atoms are spread in the plasma, the plasma temperature drops rapidly in the core region. The sharp rise seen in the magnetic energies (Fig. [START_REF] Izzo | The role of MHD in 3d aspects of massive gas injection[END_REF]) is also associated with this phenomenon. Such rapid cooling of plasma leads to the loss of its confinement. At t ≈ 2.93 ms m = 1, n = 1 mode is seen to be growing and plasma remains ergodized. The third row in Fig. [START_REF] Lehnen | Disruptions in ITER and strategies for their control and mitigation[END_REF] shows the shock-capturing term τ sc plotted at the poloidal plane ϕ = 0 and some time instants. The term is seen to detect the location of discontinuities in pressure, takes high values only in the vicinity, and therefore highlights the location at which artificial viscosity is being added. In all figures, the discontinuities in pressure are detected near boundaries and inside the domain near the high values of the impurities. Detection near the boundaries is associated with the initial plasma profile in which the pressure has a jump. Detections inside the domain are associated with the moving shocks arising due to gas-plasma interactions during the SPI simulations. The artificial viscosity added with such an adaptive coefficient improves the stability of the numerical method and allows the simulation to go beyond 2 ms with reasonable time steps. Such a stabilization mechanism is not provided by the underlying FEM or VMS-based stabilization.

Figure (14a) shows the evolution of the plasma temperature at the magnetic axis (t 0 ) in keV which decreases suddenly between 2 and 3 ms. The decrease in t 0 is associated with the rapid cooling of plasma, a distinct phase of a typical disruption which is TQ. During this phase large variations in the ablation rate (Fig. (9b)), sharp rise in the total atoms ablated (Fig. (9c)) and increase in the magnetic energies (Fig. (10c)) is seen. The similar drop in t 0 has also been noted in [START_REF] Hoelzl | First predictive simulations for deuterium shattered pellet injection in asdex upgrade[END_REF] for deuterium SPI in ASDEX Upgrade and in [START_REF] Hu | Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection[END_REF] for neon SPI in ITER plasma. Figure (14b) shows the comparison of the evolution of the thermal energies in the simulations. In the stabilized simulation with c = 10, around t ≈ 2.25 ms, the plasma thermal energy drops rapidly. This drop corresponds to the rapid decrease in T 0 as shown in Figure (14a). It is important to note here that the thermal energies in stabilized and non-stabilized simulations are comparable which denotes that the shock-capturing stabilization does not alter the underlying physics. Figure (14c) shows the evolution of the plasma current I p where I p drops slowly during the simulations and show a sign of a spike in I p , a phenomenon also seen in experiments, at the end of TQ. Such behavior of I p is also seen in the SPI simulations with the reduced MHD model presented in [START_REF] Hu | Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection[END_REF]. Once again it can be seen that I p evolution in all simulations is almost identical denoting that the numerical stabilization does not alter the underlying physics. In experiments, CQ phase follows after I p spike whose prediction is not an aim of this study, mainly because the physical model used does not include the Ohmic heating term, which can not be ignored especially after TQ when plasma temperature drops down and as a result the resistivity increases. Note that the plasma thermal energy and I p are the integrated quantities over the whole plasma domain.

In Table [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF] an estimate of the computational time taken for each simulation to reach t ≈ 1.5 ms is shown. During each time step, the linear system is solved using the GMRES method, and the average time taken by the GMRES solver at each time step is noted in Table [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF]. This time is multiplied by the number of time steps taken by each simulation to reach t ≈ 1.5 ms. It can be seen that the simulation with c = 10 takes less than twice the time taken by the non-stabilized simulation to reach t ≈ 1.5 ms to obtain the estimate for the computational time taken. This estimate is written in hours in the last column of Table [START_REF] Hoelzl | The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas[END_REF]. As the amount of the numerical stabilization increases, the average time taken by the GMRES solver also increases however, due to the ability to use larger time steps in the stabilized simulations, we gain in the overall computational time required. The stabilized simulation with c = 10 shows more than twice the speed up as compared to the non-stabilized simulation. Whereas the simulation with c = 50 takes higher time than with c = 10 as higher stabilization takes a longer time for GMRES to

Another shock-capturing method

We demonstrated the advantage of the stabilized FEM on the challenging problem above and now we try to elevate the difficulty of the problem by performing the simulation with further lower resistivity. By keeping all other parameters and the grid the same as above, we perform simulations by specifying the resistivity by an order lower i.e. η 0 = 4.22×10 -9 Ohm.m. Although grid resolution for this simulation may not be sufficient, we hope that the numerical stabilization will include the effect of unresolved scales on the resolved scales. From the numerical experiments, we found that the shock capturing term given by Eq. ( 32) does not provide the required stability. Therefore, based on experience, we use a different strategy for the modulation of the shock sensor that gives stronger estimates of τ sc . The form of τ sc used in this subsection is written as:

0 1 2 3 t [ms] 1 2 3 T 0 [keV] C=0 C=1 C=10 C=50 (a
τ sc = h 2 e β sc f p with β sc = 1 p ∂p ∂ρ L(ρ) + ∂p ∂ρ f L(ρ f ) + 1 ρ ∂p ∂t L(p) (49)
The choice of β sc is motivated by the fact that it gives higher estimates of τ sc and hence provides higher estimates of stabilization near shocks/discontinuities. It helps stabilize the spurious waves generated by interpolation of Atomic data in the impurities modeling as well. To compare the effect of different shock-capturing modulations, below we discuss the simulations listed in Table (2). In Fig. [START_REF] Neumann | A method for the numerical calculations of hydrodynamical shocks[END_REF] the evolution of the ablation rates and total impurities deposited are plotted for simulations A, B, and C. Behavior of the ablation rate, including large variations after 2 ms, is similar in simulations A and B. In both simulations, these large variations are accompanied by the rapid deposition of the impurities (see Fig. (15b)). In simulation C, the smooth evolution of the ablation rate and total ablation rate is seen. This is because the shock-detection criterion given by Eq. ( 33) with c ∥ = c ⊥ = c = 10 gives a very diffusive numerical scheme. This simulation runs with ∆t * = 5 normalized units throughout, which is significantly faster than any other simulation performed in this work. Figure (17a) and (17b) shows evolution of t 0 [keV] at the magnetic axis and thermal energies [MJ] respectively in the three simulations. In simulation C, t 0 does not drop as rapidly as in the other two simulations. A similar trend is seen in the evolution of thermal energies. Comparable evolution of t 0 and thermal energies are obtained in simulations A and B. Figure (17c) shows the evolution of I p in the three simulations where simulation A and B shows a sign of I p -spike at the end of the thermal quench. Since high numerical stabilization smooths the dynamics, the thermal quench and I p spike are absent in simulation C.

Simulation η 0 [Ohm.m] Shock-strength criterion c ∥ = c ⊥ = c A 4.22 × 10 -8 Eq. ( 33
Rows in Fig. [START_REF] Godunov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF] shows the visualization of the plasma temperature [keV] in the simulations A, B, and C respectively, at some time instants. Again the visualizations from simulations A and B look comparable. The plasma temperature starts to drop from the outer region as neon SPI fragments travel in the tokamak plasma core. With the rapid spread of the neon atoms in the tokamak plasma core after 2 ms the MHD dynamics increase which is seen from the small-scale structures in 3 rd and 4 th columns of Fig. [START_REF] Godunov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF]. During this phase, the plasma temperature is seen to drop rapidly. Such a sequence of dynamics was also expected in simulation C but it is smoothed out because of the high numerical stabilization used via the new shock-capturing method (Eq. ( 49)).

Figure [START_REF] Donea | Time-accurate solution of advectiondiffusion problems by finite elements[END_REF] shows the evolution of the magnetic energies in simulations B and C. Comparing the magnetic energies in simulation B (Fig. can be seen that the sharp rise in the magnetic energies appears nearly at the same time. Effect due to the lack of grid resolution in simulation C is visible in Fig. (19b) where the magnetic energies of higher modes are seen to be crossing those of lower modes in Fig. [START_REF] Donea | Time-accurate solution of advectiondiffusion problems by finite elements[END_REF]. Nevertheless, a sharp rise in the magnetic energies of about 2 ms is also seen in simulation C. We hope that lowering the value of c will give the numerical solution that is well resolved for the grid used.

The numerical stabilization methods in this work are designed such that as the grid is refined (h e → 0) the stabilization will go to 0.

Summary and Perspectives

Stabilized bi-cubic Hermite Bézier FEM is presented with an application to gas-plasma interaction occurring during MMI in a tokamak plasma. The stabilization for the Galerkin FEM is formulated using two methodologies: first, the VMS formulation that considers the effect of unresolved scales on the resolved scales to add upwinding in the underlying FEM, and second, the shock-capturing method that adaptively adds artificial viscosity only in the vicinity of 0 1 the shocks/discontinuities. The stabilization methodologies are applied to the bi-cubic Hermite Bézier FEM and implemented in the computational set-up of JOREK to solve MHD equations numerically. These MHD equations model complex, nonlinear gas-plasma interactions arising in the tokamak plasma. In applications considered here, nonlinear modeling of resistivities has strong local variations and gas-plasma interactions at different densities produce shocks in the solution. Such phenomena restrict the stability of the numerical methods severely making the simulations extremely expensive. The use of stabilized bi-cubic Hermite Bézier FEM is demonstrated by simulating SPI in JET-like tokamak plasma. It is seen that the shock-capturing strategy used can improve the stability of the bi-cubic Hermite Bézier FEM and reduce the computational cost required in simulating the complex physics of disruptions. Simulations presented here show physical features, including TQ, are similar to those in the previous studies performed with a reduced MHD model [START_REF] Hu | Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection[END_REF].

The numerical strategies presented here address the issues of numerical challenges faced in the gas-plasma interactions in tokamak plasma and can be applied to other tokamak applications involving shocks/discontinuities as well. Finding the optimal estimates of the stabilization matrix/coefficients remains a challenge because of the highly nonlinear plasma flows and the large number of parameters involved. Furthermore, the full MHD model used here does not include the Ohmic heating term in the energy equation. At low resistivities, its contribution is generally negligible. However, during TQ the plasma temperature decreases, and as a result the resistivity increases. The effect of the Ohmic heating term can be significant in disruptions simulations to capture the dynamics of CQ, the phase that follows TQ. Finding better estimates of the stabilization coefficients for low-resolution meshes and implementation of missing terms are reserved for future work.

The stabilized FEM presented can be easily extended for the two-temperature reduced and full MHD models. MMI in tokamak plasma applications involves many parameters such as the amount of the material injected; the number of fragments; their positions and velocities; configuration of the main plasma equilibrium; physical diffusivities to name a few. A scan in these parameters can be performed to gain insight into the MHD instabilities and complex physics involved. Further, the developed stabilized method can be extended straightforwardly to recently implemented higher-order Hermite Bézier FEMs [START_REF] Pamela | A generalised formulation of g-continuous bezier elements applied to non-linear mhd simulations[END_REF] in JOREK which are shown to be more computationally efficient than the original C 1 variant. Equipped with the numerical stabilization strategies developed here, stabilized higher-order FEM can lead to an effective and robust numerical tool for tokamak applications with large computational domains such as JET and ITER. While adaptive numerical stabilization is shown to gain computational efficiency, for larger applications such as ITER or DEMO, this alone may not be sufficient. Other techniques such as adaptive mesh refinement and multi-grid methods may be considered for further improvement of the computational efficiency.

X

  = X ∥ (b b b ⊗ b b b) + X ⊥ (I -b b b ⊗ b b b) where b b b = b b b ∥b b b∥ and I is the identity matrix.

  = ( n i e e e z , 0, 0, 0, 0) t w w w ⋆ ψ = (rn i e e e ϕ , 0, 0, 0, 0) t w w w ⋆ ρ = ( 0, 0, n i , 0, 0) t w w w ⋆ t = ( 0, 0, 0, n i , 0) t w w w ⋆ vr = (0, n i e e e r , 0, 0, 0) t w w w ⋆ vz = (0, n i e e e z , 0, 0, 0) t w w w ⋆ v ϕ = (0, n i e e e ϕ , 0, 0, 0)

Figure 1 :

 1 Figure 1: An example of annular mesh in the parametric (left) and in the physical (right) spaces. The color map in both cases denotes the values of the small radius (r).

Figure 2 :

 2 Figure 2: JET-like plasma equilibrium : (a) The profiles characterising plasma equilibrium b) boundary conditions on ψ to be specified for a polar grid plotted vs azimuthal direction θ.

Figure 3 :

 3 Figure 3: a) The contours for the numerical solution of GSE plotted on the top of the polar grid used. The solution field shows one minimum and one saddle point where ψ = 0. b) Multi-block grid with grid lines of the each block is aligned to the numerical solution of GSE. Red contours denote the surfaces where safety factor q takes rational values equal to 1 to 7 (radially outwards). (c) A schematic of SPI into JET-like plasma: JOREK grid denoting JET-like computational domain. The SPI fragments are shattered into a curved tube before releasing into the tokamak plasma. The red lines show the cone in which SPI fragments travel in the tokamak plasma.

Figure 4 :

 4 Figure 4: Numerical results axisymmetric source term modeling argon SPI in JET-like plasma computed with different shock-capturing stabilization coefficient c ∥ = c ⊥ = c. Evolution of (a) ablation rates (b) thermal energies for argon SPI simulations with the axisymmetric source.

Figure 5 :

 5 Figure 5: Numerical results axisymmetric source term modeling argon SPI in JET-like plasma computed with different shock-capturing stabilization coefficient c ∥ = c ⊥ = c. (a) The total particle density n 0 [10 20 ] plotted at ϕ = 0 plane and t ≈ 0.562 ms. The line is traced approximately along the direction of SPI velocity and passing through the location of the source term. (b) The pressure plotted along this line at the same time instant.

Figure 6 :

 6 Figure 6: Numerical results for the axisymmetric source term modeling of argon SPI in JET-like plasma computed with different shock-capturing stabilization coefficient c = c ∥ = c ⊥ . Evolution of (a) ablation rates (b) thermal energies for argon SPI simulations with the axisymmetric source.

Figure 7 :

 7 Figure 7: Numerical results axisymmetric source term modeling argon SPI in JET-like plasma computed with different shock-capturing stabilization coefficient c = c ∥ = c ⊥ . (a) The particle density n f [10 20 ] plotted at ϕ = 0 plane and t ≈ 0.562 ms. The line is traced approximately along the direction of SPI velocity and passing through the location of the source term. (b) The pressure plotted along this line at the same time instant.

Figure 8 :

 8 Figure 8: (a) At the beginning of SPI simulations some fragments (denoted by black dots) are spatially distributed (they may even be outside the plasma). The colors of the contours denote the magnitude of t in keV. (b) 3D temperature field at the beginning of the simulation showing axisymmetry.

Figure 9 :

 9 Figure 9: Comparison of the simulations of Neon SPI in JET-like plasma without and with numerical stabilization. The following quantities are plotted vs t [ms] (a) The number of time steps (N t ) taken by the simulations. The labels to the curves are the values of the time step ∆t ⋆ written in the normalized units for the sake of the presentation. The scale to convert ∆t ⋆ in ms is = 2.973 × 10 -4 and hence ∆t = ∆t * × 2.97 × 10 -4 ms. (b) The ablation rate. (c) The total number of impurity atoms ablated in the plasma.

= 50 Figure 10 :

 5010 Figure 10: Evolution of the magnetic energies in each toroidal harmonic in the SPI simulation with η 0 = 4.22 × 10 -8 [Ohm.m] and c ∥ = c ⊥ = c.

Figure 11 :

 11 Figure 11: Snapshots of the impurity particle density n f [10 20 m -3 ] plotted in the poloidal plane ϕ = 0 with the shock-capturing stabilization coefficient c ∥ = c ⊥ = c. Top and bottom rows show visualization at t ≈ 1.6 and 2.73 ms respectively.

Figure 12 :

 12 Figure 12: 3D visualizations of SPI in JET-like plasma at some time instants. The pseudo-color plotted at two planes 180 degrees apart from each other denotes the values of T in keV. The green and red iso-surfaces denote the impurity particle density n f 35 % and 70 % of its maximum value respectively.

Figure 13 :

 13 Figure 13: Snapshots of the impurity particle density n f [10 20 ] (first line), Poincaré plots of the magnetic field (second line), and the shock-capturing parameter τ sc (last line) at the mentioned time instants. The figures in the first and second lines are plotted in the poloidal plane ϕ = 0. The black dots in the figures with Poincaré plots denote the positions (projected onto the poloidal plane) of the SPI fragments while the color denotes values of the plasma temperature in keV. These visualizations are from the simulation with the shock-capturing stabilization coefficient c ∥ = c ⊥ = c = 10.

Figure 15 :

 15 Figure 15: Evolution of a) the ablation rate and b) the total number of impurity atoms ablated in the plasma during Neon SPI simulations in JET-like plasma with different resistivities and shock-capturing coefficient c ∥ = c ⊥ = c.

Figure ( 16 )

 16 Figure[START_REF] Wilkins | Calculation of elastic plastic flow[END_REF] shows snapshots of the impurity particle density (n f ) [10 20 m -3 ] plotted in ϕ = 0 plane for the three simulations at approximately same time instants. Visualizations from different physical phases are comparable for simulations A and B which use different shock-capturing methods. Both simulations show fast dynamics over which neon atoms are spread rapidly to cool the plasma and lead to thermal quench. While the use of very high numerical stabilization smooths the dynamics in simulation C and no distinct phase of thermal quench is seen.Figure (17a) and (17b) shows evolution of t 0 [keV] at the magnetic axis and thermal energies [MJ] respectively in the three simulations. In simulation C, t 0 does not drop as rapidly as in the other two simulations. A similar trend is seen in the evolution of thermal energies. Comparable evolution of t 0 and thermal energies are obtained in simulations A and B. Figure(17c) shows the evolution of I p in the three simulations where simulation A and B shows a sign of I p -spike at the end of the thermal quench. Since high numerical stabilization smooths the dynamics, the thermal quench and I p spike are absent in simulation C.Rows in Fig.[START_REF] Godunov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF] shows the visualization of the plasma temperature [keV] in the simulations A, B, and C respectively, at some time instants. Again the visualizations from simulations A and B look comparable. The plasma temperature starts to drop from the outer region as neon SPI fragments travel in the tokamak plasma core. With the rapid spread of the neon atoms in the tokamak plasma core after 2 ms the MHD dynamics increase which is seen from the small-scale structures in 3 rd and 4 th columns of Fig.[START_REF] Godunov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF]. During this phase, the plasma temperature is seen to drop rapidly. Such a sequence of dynamics was also expected in simulation C but it is smoothed out because of the high numerical stabilization used via the new shock-capturing method (Eq. (49)).Figure[START_REF] Donea | Time-accurate solution of advectiondiffusion problems by finite elements[END_REF] shows the evolution of the magnetic energies in simulations B and C. Comparing the magnetic energies in simulation B (Fig. (19a)) with those in simulation A (Fig. (10b)) it

  Figure[START_REF] Wilkins | Calculation of elastic plastic flow[END_REF] shows snapshots of the impurity particle density (n f ) [10 20 m -3 ] plotted in ϕ = 0 plane for the three simulations at approximately same time instants. Visualizations from different physical phases are comparable for simulations A and B which use different shock-capturing methods. Both simulations show fast dynamics over which neon atoms are spread rapidly to cool the plasma and lead to thermal quench. While the use of very high numerical stabilization smooths the dynamics in simulation C and no distinct phase of thermal quench is seen.Figure (17a) and (17b) shows evolution of t 0 [keV] at the magnetic axis and thermal energies [MJ] respectively in the three simulations. In simulation C, t 0 does not drop as rapidly as in the other two simulations. A similar trend is seen in the evolution of thermal energies. Comparable evolution of t 0 and thermal energies are obtained in simulations A and B. Figure(17c) shows the evolution of I p in the three simulations where simulation A and B shows a sign of I p -spike at the end of the thermal quench. Since high numerical stabilization smooths the dynamics, the thermal quench and I p spike are absent in simulation C.Rows in Fig.[START_REF] Godunov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF] shows the visualization of the plasma temperature [keV] in the simulations A, B, and C respectively, at some time instants. Again the visualizations from simulations A and B look comparable. The plasma temperature starts to drop from the outer region as neon SPI fragments travel in the tokamak plasma core. With the rapid spread of the neon atoms in the tokamak plasma core after 2 ms the MHD dynamics increase which is seen from the small-scale structures in 3 rd and 4 th columns of Fig.[START_REF] Godunov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF]. During this phase, the plasma temperature is seen to drop rapidly. Such a sequence of dynamics was also expected in simulation C but it is smoothed out because of the high numerical stabilization used via the new shock-capturing method (Eq. (49)).Figure[START_REF] Donea | Time-accurate solution of advectiondiffusion problems by finite elements[END_REF] shows the evolution of the magnetic energies in simulations B and C. Comparing the magnetic energies in simulation B (Fig. (19a)) with those in simulation A (Fig. (10b)) it

  (a) t ≈ 2.18 ms (b) t

Figure 16 :Figure 17 :

 1617 Figure 16: Snapshots of the impurity particle density (n f ) [10 20 ] m -3 at some time instants. Plots in the first, second and third rows are from the simulations A, B and C respectively listed in Table (2). All the plots show visualizations in the poloidal plane ϕ = 0.

(

Figure 18 :

 18 Figure 18: Snapshots of the temperature T [keV] at some time instants. Plots in the first, second and third rows are from the simulations A, B and C respectively listed in Table (2). All the plots show visualizations in the poloidal plane ϕ = 0.

Figure 19 :

 19 Figure 19: Evolution of the magnetic energies in each toroidal harmonic during neon SPI simulations in JET-like plasma with different resistivities and shock-capturing coefficient c ∥ = c ⊥ = c.

Table 1 :

 1 Time taken by GMRES solver in different simulations with c ∥ = c ⊥ = c to reach the physical time t ≈ 1.5 ms using 64 MPI processes.

		0.1 0.2 0.3 0.4 0.5 Thermal Energy [MJ]	C=0 C=1 C=10 C=50					2.2 2.4 1.6 2.0 Ip [MA] 1.8	C=0 C=1 C=10 C=50
			0	1	t [ms]	2	3	0	1	t [ms]	2	3
	)				(b)					(c)
	Figure 14: Evolution of the (a) plasma temperature at the magnetic axis [keV] (b) plasma
	thermal energy [MJ] and (c) plasma current I p [MA].
	c no of time steps Average time in [s] re-	Total GMRES
			quired for one GMERS	time in hours ≈
			solve
	0	8096	6.0461					13.6
	1	3988	10.6423				11.8
	10	1182	15.2078				5
	50	1182	17.6222				5.78

Table 2 :

 2 Parameters to identify the simulations for the discussion the effect of different shockcapturing methods differing in shock modulations. All other parameters and the grid used in these simulations are identical.

	)	1

Acknowledgment

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 -EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. The simulations presented here have been performed using the ITER HPC cluster, Marconi-Fusion supercomputer operated by Cineca, Italy, and JFRS-1 operated by IFERC-CSC in Japan.

A Computation of integrated quantities

The integrated quantities plotted in the Figures are computed as follows. The thermal energy (TE) is computed as p dΩ where p is the total pressure given by Eq. [START_REF] Baylor | Shattered pellet injection technology design and characterization for disruption mitigation experiments[END_REF]. Magnetic energies (ME) in each toroidal harmonic n are computed as 1 J ϕ dΩ where J ϕ is ϕ-component of the current vector j j j.