# A high-order residual-based viscosity finite element method for the ideal MHD equations

We present a high order, robust, and stable shock-capturing technique for finite element approximations of ideal MHD. The method uses continuous Lagrange polynomials in space and explicit Runge-Kutta schemes in time. The shock-capturing term is based on the residual of MHD which tracks the shock and discontinuity positions, and adds a sufficient amount of viscosity to stabilize them. The method is tested up to third order polynomial spaces and an expected fourth-order convergence rate is obtained for smooth problems. Several discontinuous benchmarks such as Orszag-Tang, MHD rotor, Brio-Wu problems are solved in one, two, and three spatial dimensions. Sharp shocks and discontinuity resolutions are obtained.

## Authors

• 3 publications
• 7 publications
07/22/2019

### A conforming DG method for the biharmonic equation on polytopal meshes

A conforming discontinuous Galerkin finite element method is introduced ...
12/06/2021

### A high order explicit time finite element method for the acoustic wave equation with discontinuous coefficients

In this paper, we propose a novel high order explicit time discretizatio...
12/09/2019

### An adaptive high-order piecewise polynomial based sparse grid collocation method with applications

This paper constructs adaptive sparse grid collocation method onto arbit...
11/01/2020

### A Riemann Difference Scheme for Shock Capturing in Discontinuous Finite Element Methods

We present a novel structure-preserving numerical scheme for discontinuo...
06/20/2022

### Monolithic parabolic regularization of the MHD equations and entropy principles

We show at the PDE level that the monolithic parabolic regularization of...
04/27/2021

### Simulation of non-linear structural elastodynamic and impact problems using minimum energy and simultaneous diagonalization high-order bases

We present the application of simultaneous diagonalization and minimum e...
10/17/2020

### Arbitrarily high-order exponential cut-off methods for preserving maximum principle of parabolic equations

A new class of high-order maximum principle preserving numerical methods...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Conservation laws play an important role in modeling and understanding physical processes. For example, conservation of mass, energy, and momentum are fundamental principles in gas and fluid dynamics. For ideal flow, these principles are modeled using the well-known system of Euler equations. Despite a long history and extensive research in understanding and solving Euler equations, many important mathematical and numerical questions remain unsolved. Especially for high-speed flows (higher than the speed of sound), solutions to the Euler system lose their regularity and many complex phenomena such as shocks, strong discontinuities, rarefaction, and contact waves develop. Besides, at very high temperature and energy gas transforms into plasma. Plasma is an ionized gas and the dynamic of plasmas is modeled by adding additional terms and equations involving the magnetic field into the Euler system. This resulting system is usually referred to as Magnetohydrodynamics or MHD in short.

Numerical approximation of MHD started as soon as researchers started to simulate the Euler equations, see (Lax_1954; Brackbill_1976). Existing numerical methods for compressible Euler equations are hard to extend to approximate the MHD system. If one reason is solving -additional equations, where is the space dimension, other two main reasons are developing high order robust shock-capturing techniques, satisfying the divergence-free constraint for the magnetic field, . State-of-the-art numerical methods to solve the MHD equations are finite differences, finite volumes, and discontinuous Galerkin (DG) schemes. Most of the finite difference and finite volume schemes are based on approximate solutions to the Riemann problem of MHD, see for instance (Brio_Wu_1988; Powell_et_al_1991; Dai_Woodward_1998; Balsara_et_al_1999; Bouchut_et_al_2007; Balsara_2010; Balsara_Dubser_Abgrall_2014) and references therein. Due to similarities of finite volume and DG schemes, an approximate Riemann solver approach was incorporated in the Galerkin variational formulation for the DG schemes, see (Warburton_Karniadakis_1999; Li_Shu_2005; Dumbser_2016). In general, the exact or approximate Riemann solvers are difficult to compute. Especially, the MHD system is not strictly hyperbolic and has a non-convex flux, and therefore it may admit non-regular waves which may result in non-unique Riemann solution (Torrilhon_2002). To overcome this difficulty, numerical methods that do not use Riemann solvers have been developed. For instance central DG schemes reported in (Tadmor_mhd1_2004; Li_Xu_Yakovlev_2011; Cheng_et_al_2013) approximate the MHD solutions without using any approximate or exact Riemann solvers.

The success of finite volume and DG schemes has not been observed in a continuous finite element (FE) framework. A reason for it could be the lack of high order, robust explicit stabilization techniques to stabilize continuous FE approximations of MHD. Traditional stabilization methods for FE approximations of conservation laws including Euler and MHD equations are based on the Galerkin-Least-Squares argument (GLS) augmented with residual-based shock-capturing terms, see for instance (Szepessy_1989; Johnson_et_al_1990) and references therein. GLS schemes are implicit by construction, difficult to make high order in time, and nontrivial to implement since the test functions are required to be time-dependent. Only a few works were done to simulate resistive MHD using GLS stabilization, see for instance (Badia_Codina_2013; Sitaraman_Raja_2013; Shadid_et_al_2016) and references therein.

It appears that a way of constructing a high-order explicit FE approximation of conservation laws could be by entirely disregarding the GLS terms from the stabilized formulation. Recently, Guermond_2008; Guermond_pasquetti_popov_JCP_2011; GuermondNaPo11; Nazarov_Larcher_2017 proposed the so-called entropy viscosity method, where the FE approximation is stabilized by adding an elliptic term, and the artificial viscosity coefficient is constructed to be proportional to an entropy residual of the system. For scalar conservation laws, the FE residual is one of the entropy residuals. Therefore, the residual-based artificial viscosity method can be obtained by disregarding the GLS terms in the stabilization method of (Szepessy_1989). In (Nazarov_2013)

, we proved that the FE approximation of residual-based artificial viscosity method applied to scalar conservation laws converges to the unique entropy solution, and we extended the method to solve more general systems in the framework of DG, spectral elements, finite differences, and radial basis functions (RBF), see

(Nazarov_Hoffman_2013; Marras_2015; Lu_2019; Striernstrom_2021; Tominec_2021). For other approaches to constructing the nonlinear artificial viscosity methods, we refer to the work of (Basting_2017; Kuzmin_2020; Mabuza_2020), where the stabilization is constructed based on a smoothness indicator of the solution.

In this work, we propose a new residual-based shock-capturing method for FE approximation of the MHD system. The method is explicit in time and works for high order polynomial degrees. The divergence-free constraint is obtained using several techniques, including the projection of into a divergence-free space (Brackbill_Barnes_1980) and the hyperbolic cleaning (Dedner_et_al_2002). We should emphasize that positivity and invariant domain preservation of the method as in (Guermond_et_al_2018) is left beyond the scope of this paper.

The paper is organized as follows: in Section 2, we introduce important notations and terms that are used in the paper, such as FE spaces, meshes, and matrices. In Section 3

, we introduce the ideal MHD equations, eigenvalues, and the divergence cleaning techniques. In Section

4, we present a viscous regularization of the MHD equations using a novel residual-based shock-capturing method. Finally, in Section 5, we solve several well-known benchmark problems for ideal MHD. The method is shown to have an optimal convergence rate for smooth problems and captures discontinuous waves for non-smooth tests. Section 6 is the conclusion.

## 2 Continuous Galerkin finite element method

In this section, we introduce the finite element spaces and notations that are used in this work.

Let us consider an open bounded domain , with boundary . We denote a finite partition of into disjoint elements such that no vertex of any element is placed on the interior of an edge of another element. The union of all the closure of elements constitutes , where denotes the closure of .

We seek finite element solution in Lagrange polynomial spaces

 \calQh={v(x):v∈\calC0(Ω),v|K∈\polPk,∀K∈\calTh}, (1)

where is a space of polynomials of at most degrees. Let be the -degree Lagrange basis function in which takes value 1 at node , and 0 at any other nodes. The functions form a basis for , where is the total number of nodes in . We define the set of all nodal points contained within the support of . Furthermore, we use the following Hilbert space inner products,

We will use a mesh-function which we compute using the following projection: find such that

 ∑K∈\calTh∫Khhv\ud\bx=∑K∈\calTh∫KhKkv\ud\bx,∀v∈\calQh, (2)

where is the circumradius of the element .

## 3 Governing equations of MHD

Let us denote the time interval , where is the final time. For all , we seek solution of the MHD equations , , where is the mass density, is the momentum, is the total energy, is the magnetic field. The ideal MHD equations in conservative form with appropriate boundary conditions are defined as follows:

 \pt\bsfU+\DIV\bsfF\calE(\bsfU)+\DIV\bsfF\calB(\bsfU)=0, (3)

where the nonlinear tensor fluxes

are defined as

 \bsfF\calE:=⎛⎜ ⎜ ⎜ ⎜⎝\bbm\bbm⊗\bu+p\polI\bu(E+p)0⎞⎟ ⎟ ⎟ ⎟⎠,\bsfF\calB:=⎛⎜ ⎜ ⎜ ⎜⎝0−\bbetaa−\bu\SCAL\bbetaa\bu⊗\bB−\bB⊗\bu⎞⎟ ⎟ ⎟ ⎟⎠, (4)

with being the velocity field, and being the Maxwell stress tensor:

 \bbetaa=(−12(\bB\SCAL\bB)\polI+\bB⊗\bB),

where

denotes the identity matrix of size

. The combination of hydrodynamics and electromagnetism nature of the MHD system can be seen from the MHD equations in the form of (3), through the separated fluxes and . The thermodynamic pressure is computed from the equation of state

 p=(γ−1)e, (5)

where is the internal energy and is the adiabatic gas constant. The internal energy is computed as

 e=E−12ρ(\bu\SCAL\bu)−12(\bB\SCAL\bB).

In addition, the magnetic field should satisfy the following divergence free constraint:

 \DIV\bB=0. (6)

Note that the above equations are in non-dimensional form, as suggested by (Sitaraman_Raja_2013, Sec. 3.1).

### 3.1 Eigenvalues of the inviscid MHD

We are going to use the maximum wave speeds to construct our stabilization terms. A good approximation of the local wave speeds is obtained by eigenvalues of the inviscid MHD equations. It is well-known that eight eigenvalues are corresponding to the eight elementary waves for the ideal MHD equations, see Powell_et_al_1991; Barth_1999; Rossmanith_2006: let

be a direction vector,

be the speed of sound, and

The eigenvalues from the smallest to the largest are

 λ1,8=\bu\SCAL\be∓cf,λ2,7=\bu\SCAL\be∓b,λ3,6=\bu\SCAL\be∓cs,λ4,5=\bu\SCAL\be, (7)

where the subscripts , and correspond to fast and slow magnetosonic waves, and the Alfén waves. The characteristics correspond to entropy and divergence waves. It is clear to see that the first and eighth eigenvalues are the ones of interest, which represent the fastest moving waves. Later in Section 4, we will use the maximum wave speed

 λmax:=maxi=1,…,8|λi|

in the construction of the nonlinear viscosity term.

### 3.2 Divergence cleaning techniques

To satisfy the magnetic divergence-free condition (6), we have verified that our proposed nonlinear stabilization can be used along with popular techniques for divergence cleaning: the elliptic correction/projection method by Brackbill_Barnes_1980, and the hyperbolic correction by Dedner_et_al_2002; Tricco_2016. A numerical comparison between the mentioned techniques is included in Section 5.3.

#### 3.2.1 Elliptic correction/Projection method

In each time step, instead of using the magnetic solution from the system, the projection of onto a magnetic divergence-free space is used. Brackbill_Barnes_1980 proposed using

 \bB′=\bB−∇Ψp,

where is the solution to the Poisson equation . This yields as the non-physical part is removed from the numerical solution . After each correction in the magnetic field , the dependent variables: pressure , temperature , energy and the entropy functions are updated accordingly to ensure consistency of the discrete solution. Contrary to the seemingly ad-hoc nature of this method, Toth_2000 shows that the projection method preserves conservative properties and accuracy of the base scheme.

Pseudo time-stepping. Numerically minimizing by the projection method is considered computationally expensive. A more affordable alternative, proposed by Hayashi_2005, is solving the following time-dependent equation for ,

 ∂˜tΨp−∇2Ψp+\DIV\bB=0, (8)

where the variable is a separate time dimension. At equilibrium, the solution to (8) is the one to the Poisson equation. A semi-discretization of (8) is given by

 Ψℓ+1p−Ψℓp˜τ−∇2Ψℓ+1p+\DIV\bBℓ=0,ℓ=0,1,2,…,

where is the pseudo step size. After each pseudo time step, the magnetic solution is corrected . The divergence is gradually decreased. One can think of this technique as applying the elliptic correction multiple times. The computational cost depends on how many pseudo steps are performed.

#### 3.2.2 Hyperbolic correction

An auxiliary scalar variable is added to the system to transport and suppress the divergence violation on the magnetic field. We add a fifth equation to the MHD system

 ∂tΨl+ch\DIV\bB=−crchhhΨl,

where in 2D or in 3D (Tricco_2016), and denotes the mesh function computed in (2). The coefficient is set to be the maximum wave speed = . In addition, a source term is added to the right hand side of (3),

 \pt\bsfU+\DIV\bsfF\calE(\bsfU)+\DIV\bsfF\calB(\bsfU)=\bsfGGLM% (\bsfU), (9)

where

 \bsfGGLM(\bsfU)=−⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝00\bB\SCAL∇Ψl\DIV(Ψl\polI)0⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

## 4 Nonlinear viscosity method for MHD

The aim of this paper is to solve the inviscid MHD system (3) using continuous finite element approximations. We use the following two vector valued spaces:

 \bcalVh:=[\calQh]d,\bcalWh:=\calQh\CROSS\bcalVh\CROSS\calQh\CROSS\bcalVh, (10)

where is defined in (1). Then, we formulate the finite element approximation of the MHD system (3) as follows: find such that

 (\pt\bsfUh,\bsfVh)+ (\DIV\bsfF\calE(\bsfUh),\bsfVh)+(\DIV\bsfF\calB(\bsfUh),\bsfVh)=0,∀\bsfVh∈\bcalWh. (11)

Here the inner products are computed exactly using appropriate numerical quadrature rules. The finite element approximation of the flux terms in (11) is equivalent to central difference schemes, therefore produces spurious oscillations and cannot be used for strong discontinuities such as shocks and contact lines. Below, we introduce a nonlinear artificial viscosity approach to fix this issue.

### 4.1 New class of viscous regularization of the MHD equations

In this section, we propose using the following vanishing viscosity regularization of the MHD system. The new regularization is obtained by adding an elliptic term into the MHD system, where denotes a viscous flux tensor.

As a heuristic argument to construct

, we start from the viscous flux of the resistive MHD equations, see e.g., (Dumbser_2009),

where the temperature is a function in . The viscous shear stress tensor is defined as

where is the artificial dynamic viscosity, is the bulk viscosity of undetermined sign such that , is the artificial heat conduction, is the artificial electric resistivity of the medium. Note that the electric resistivity has the same dimension as the dynamic viscosity coefficient . We point out that the temperature and velocity field are computed at the nodal points: , , .

By considering zero magnetic field, i.e., and thus , we obtain from (3) the compressible Euler equations. For the compressible Euler equations, a frequent choice of viscous regularization is the Navier-Stokes flux, see e.g., (Nazarov_Larcher_2017), which coincides with when . The Navier-Stokes flux does not add any viscosity term to the mass equation. Therefore, positivity of density can be violated (Guermond_2014; Nazarov_Larcher_2017). In Guermond_2014, the authors have shown the importance of regularizing the mass equation and have proven that it is a key argument for keeping density and internal energy positive, and for the minimum entropy principle. Numerical validation of the viscous regularization compressible Euler equation including the mass equation in the context of high-order finite element approximation is reported in (Nazarov_Larcher_2017). They show that the new viscous flux performs better than the traditional Navier-Stokes flux in capturing complex phenomena, e.g., shocks, rarefaction, contact lines, while maintaining high order accuracy.

For MHD, we add a corresponding elliptic term scaling by a positive number to regularize the mass conservation. Therefore, the obtained viscous flux is a slightly modified Navier-Stokes viscous fluxes for the flow part and resistivity magnetic fluxes for the magnetic part:

Then, the viscous regularized finite element approximation of the MHD equations is: find such that

where is the normal vector pointing outward at every point on the boundary . In the rest of the paper, we look for a viscous solution, however, to make the notation easier, we drop the superscript sign.

We construct the viscosity coefficients , and , all are from the finite element space , and make them proportional to a nonlinear functional such that enough viscosity is added to the system in terms of stability and accuracy.

Let , , and be the low, high and dynamic viscosity coefficients at time , and .

### 4.2 First-order viscosity

At each node , we compute the first order viscosity at time as

 μL,ni=Cmaxhiλnmax,i,

where is the largest eigenvalue at time , given by (7), and is a positive parameter. A typical range of is reportedly , see (Nazarov_Hoffman_2013). In one space dimension, it can shown that leads to the first order Lax-Friedrichs scheme. We use in all the numerical examples in this article for all polynomial degrees. Note, that the mesh function is already scaled with the polynomial degree in (2).

This choice of first-order viscosity is consistent with the CFL condition discussed later. In addition, the largest eigenvalue turned out to be sufficient to construct a stable first-order viscosity for the numerical examples considered in this work. We emphasize that for a positivity preserving first-order scheme, similar to the Euler system, as in Guermond_Popov_2016

, a more accurate estimate of the maximum wave velocity is required.

Our goal is to build a high-order stabilized method. We use the first-order scheme as a safety factor and then construct a high-order extension of it.

### 4.3 High-order residual-based viscosity

The idea behind a high order scheme is to build an indicator that tracks strong shocks and discontinuities and activates the first-order method from the section above. At the same time, the indicator should decrease the viscosity where the solution is smooth. For this, we use the finite element residual of the PDE as an indicator, since it vanishes in smooth regions, , and has a jump of size in discontinuities.

We use the residual as an indicator, therefore, we need to get rid of the corresponding solution unit. Our aim is to find a normalization function such that it transforms any input function to the same arithmetic range, and is locally amplified at the regions where the input function changes drastically.

Consider a continuous scalar function and in particular the case . We propose dividing by the following function to normalize

 n(w)i=∣∣∣¯S(w)−α(maxj∈\calI(i)wj−minj∈\calI(i)wj)∣∣∣, (15)

where

indicates the global maximum variance

, and the global constant is calculated as

 α=Cl¯S(w)maxΩw−minΩw

which acts as a unit normalizer of local jumps to the global variance. As a technical detail, notice that when is a constant uniform function, we have for all . In that case, reduces to which although being zero, we do not suffer from division by zero in the calculation of . In an earlier development, the normalization function in (Striernstrom_2021) can be obtained by setting in (15). We observe that the new formulation (15) deals better with mixed-signed solutions since the two factors in are of the same unit. To demonstrate this argument, we rewrite the expression (15) as

 n(w)i=¯S(w)(1−Clmaxj∈\calI(i)wj−minj∈\calI(i)wjmaxΩw−minΩw). (16)

The parameter allows fine tuning of how sensitive this normalization affects the local discontinuities. Observe that the ratio is in between the interval , and is amplified at local solution jumps. If the input function holds both positive and negative parts, and if , it is possible that at the regions where the local jumps incidentally matches the global variance. However, once the local jump bypasses the global variance, can become uncontrollably big, which opposes to our aim that should be small at such circumstances. On the other hand, the formula (16) does not suffer this issue because the value within the bracket in (16) is strictly positive. How small can become is tuned with the parameter . For example, if , can become zero at the largest jump of . Lowering from towards would make at jumps of larger and larger but not exceeding . In smooth regions, turns into , which is the classical normalization function of our residual-based viscosity technique, see Nazarov_2013.

Next, we construct a normalized finite element residual of the MHD equations as

 ˜Rnh:=max(R(ρnh)n(ρnh),R(\bbmnh,1)n(\bbmnh,1),…,R(\bbmnh,d)n(\bbmnh,d),R(Enh)n(Enh),R(\bBnh,1)n(\bBnh,1),…,R(\bBnh,d)n(\bBnh,d)), (17)

where the normalization terms , as defined above, normalize the residual units and magnitude. It is also possible that the denominators in (17) become close to zero, in the cases when the initial solution or the analytic solution of the corresponding component is a constant or almost a constant. To avoid division by zero, instead of using , where is a scalar function, we use . We choose to be machine epsilon, i.e. approximately .

For each component , the residual can be calculate by standard finite element projection,

 ∫ΩR(qnh)φi\udx=∫Ω|BDF(qh)n+\DIV(f(qnh))|φi\udx, (18)

where corresponds to the nonlinear MHD flux in (11) for each component , and is a backward difference approximation of the time derivative at time . Note that the second-order BDF scheme is sufficient in getting 4th order accurate scheme.

For Lagrange elements, we observe that without destroying robustness, the residual can be calculated much more efficiently by approximating the above projection solution by lumping the mass matrix,

 R(qnh)≈1mi∫Ω|BDF(qh)n+\DIV(f(qnh))|φi\udx,

where for Lagrange polynomials.

Finally, we are now ready to construct the final artificial viscosity on each node at time :

 μni=min(μL,ni,μH,ni), (19)

where , is the nodal value of the residual , is a positive parameter. A typical range of is reportedly , see (Nazarov_Hoffman_2013). For the numerical simulations in this paper, we use . In the rest of this paper, we refer to using (19) in the formulation of the viscous regularized equation (14) as the “RV method”.

Once the MHD system is discretized in space using the continuous finite element method we obtain the system of ODEs (14). Let us denote this system as

 \polM\pt\bsfUh(t)=\calF(\bsfUh(t),\bmuh(t)),

where is the mass matrix, is the right-hand-side function of the system which depends on the solution and the viscosity vector is .

Next, discretize this system in time using explicit -stage Runge-Kutta methods:

 \bsfUn+1h=\bsfUnh+τn(b1\bsfK1+…+br\bsfKr), (20)

where , , , are coefficients obtained from the Butcher tableau, and the stage variables are computed as follows: for the given solution and the viscosity vector at time level , set , then let be the solution at the -th stage of the Runge-Kutta method, then compute by solving the following system:

 \polM\bsfKl=\bsfF(\bsfWl(t),\bmunh),

for all . Note that the viscosity coefficients can also be computed on the fly at every Runge-Kutta stage. However, in this work, the viscosity coefficient is constructed from the previous time level and does not change within the Runge-Kutta stages.

The time-step is computed adaptively using a CFL condition

 τn=CFLmin\bx∈Ωhh(\bx)maxj=1,…,Nh|λnmax,j|. (21)

### 4.5 Boundary conditions

In the following numerical examples we use the following boundary conditions: Dirichlet, Neumann, periodic, and slip or an impermeability boundary condition. The Neumann boundary condition is implemented as an additional boundary integral in the variational form. The periodic boundary condition is implemented by enforcing the same value of degrees of freedom (DOF) of matching nodes on two opposite boundaries in question.

The Dirichlet and slip boundary conditions are imposed strongly as a correction step after the Runge-Kutta method is applied. Assume that we want to compute solution at time level , , we impose the Dirichlet boundary condition by setting the values of the solution on the boundary nodes by its boundary data.

The slip or an impermeability boundary condition requires on the boundary in question. This condition is equivalent to . Then, the boundary value of momentum at the boundary node , is replaced by .

We refer the reader to Nazarov_Larcher_2017, where advantages of the strong implementation of the slip boundary condition versus its weak counterpart are discussed, and we refer to Guermond_et_al_2018, where the correction technique is discussed for explicit schemes in more detail.

### 4.6 Summary of the algorithm

To summarize, our stabilized finite element solution of the MHD system (3) can be obtained by the following steps. In each time step, the solution is advanced as

1. Calculate the residual using (18), combine the low order and high order viscosity to calculate the amount of artificial viscosity using (19).

2. Solve the main linear system (14) where the time derivative is discretized using a Runge-Kutta method (20) of order at least .

3. If the projection method or the pseudo time-stepping method is used for cleaning divergence of , use the techniques described in Section 3.2.1. If the hyperbolic correction method in Section 3.2.2 is used, divergence cleaning is already incorporated into the system in the previous step.

4. Apply boundary conditions strongly as described in Section 4.5.

5. Compute an approximation of the maximum local velocity using (7), and identify the next time step size using (21).

## 5 Numerical examples

In this section, we demonstrate the efficiency of our proposed stabilization method on several well-known benchmark problems. For the following tests, the residual viscosity parameters are . We observe that the parameter is not sensitive as in most circumstances, apart from enhanced smoothness at sharp shocks, any choice of would give satisfactory outcomes. We have used the classical Runge-Kutta 4 for all the tests with time step (21), where the CFL number is chosen to be . In the viscous flux (13), the viscosity coefficients , and are chosen to be given by (19).

### 5.1 Accuracy test

The first two test problems where an analytic solution can be derived are presented to show the high order property of the stabilization in the smooth regions.

#### 5.1.1 Smooth vortex problem, Wu_et_al_2018

Consider a periodic smooth vortex problem on a rectangle domain . The reference solution is a stationary flow with a vortex perturbation

 (ρ(t),\bu(t),p(t),\bB(t))=(ρ0,\bu0+δ\bu,p0+δp,\bB0+δ\bB),

where

 ρ0 =1, \bu0 =(1,1), δ\bu =μπ√2e(1−r2)/2(−r2,r1), p0 =0, δp =−μ2(1+r2)e1−r28π2, \bB0 =(0.1,0.1), δ\bB =μe(1−r2)/22π(−r2,r1),

the vortex radius is , , and the vortex strength is . The adiabatic constant is . The errors measured at final time using elements are shown in Table 1. In Table 1, we show behavior of both the Galerkin solutions and the RV solutions. For the solution, the obtained error is lower than that of the solution under the same number of degrees of freedom. However, the convergence rate of the is second order, which is suboptimal with regards to theoretical derivations. This is a known issue with continuous Galerkin approximations, which was earlier reported in, e.g., (Nazarov_Larcher_2017). Figure 1 illustrates and compares the convergence history between the Galerkin solutions and the RV solutions. All the obtained convergence rates show that the RV method does not destroy high-order accuracy of the Galerkin smooth solutions.

#### 5.1.2 Smooth wave propagation, Wu_et_al_2018

In this test case, we consider a rectangle domain . A periodic wave-like solution reads . The adiabatic constant is . The errors measured at final time using elements are shown in Table 2. Again, the obtained convergence rates are optimal with regards to the corresponding Galerkin approximations.

### 5.2 Brio-Wu MHD shock tube problem, Brio_Wu_1988

The RV method shows high-order accuracy for smooth problems considered in above. Now, we start solving problems with strong shocks and discontinuities.

We want to solve the one-dimensional Brio-Wu MHD shock tube problem, which was first introduced by Brio_Wu_1988. This is the one-dimensional Riemann problem in the domain . The problem is challenging because the solution contains waves of different types, typical for an ideal MHD.

The problem setting is the following:

 (ρ,u,p,Bx,By)={(1,0,1,0.75,1),x∈[0,0.5),(0.125,0,0.1,0.75,−1),x∈[0.5,1.0].

The adiabatic constant is . The importance of this test is to examine whether a MHD solver can represent the shocks, rarefactions, contact lines, and the whole compound structure accurately. The solution at the final time is shown in Figure 2. A well-known reference solution is provided by the Athena code (Stone_2008) with 10000 grid points.

The result in Figure 2 shows that our proposed method can capture the wave structures efficiently and accurately. We also see that solution is sharper than for the same number of degrees of freedom.

### 5.3 2D Orszag-Tang problem, (Orszag_Tang_1979)

In this section, we consider the well-known benchmark Orszag-Tang problem in 2D. The considered domain is a square, . The solution is initialized as

 (ρ,\bu,p,\bB)=(2536π,(−sin(2πy),sin(2πx)),512π,(−sin(2πy)√4π,sin(4πx)√4π)).

The adiabatic constant for this problem is . The periodic boundary condition is used in both directions. For divergence cleaning, we have used the projection method. The density solution and the artificial viscosity at time and are shown in Figure 3 using elements, and in Figure 4 using elements. At , the structure of the Orszag-Tang solution is considered turbulence (Tricco_2016).

Despite being a classical benchmark in studies of numerical schemes for MHD, the evolution of the Orszag-Tang solution after has rarely been reported or discussed in the literature. For the sake of the readers, we note that the time scale can be different in some other papers. For example, in (Toth_2000), the solution at approximately corresponds to our solution at . The authors in (Balsara_1998) pointed out that for short integration in time, whether or not a divergence cleaning procedure is involved does not affect the Orszag-Tang solution as there is no noticeable difference. We also observe this situation in our continuous finite element simulation both visually and numerically, at least up to . In fact, it is evident that the Orszag-Tang solution up to is relatively straightforward to be obtained by numerical schemes with shock-capturing capability, e.g., Balsara_1998; Guillet_2019; Toth_2000. However, long time integration would discriminate different discretizations, as remarked by Guillet_2019. Using DG method, the authors in Guillet_2019 reported that not all choices of the slope limiters and discretizations of the Powell term Powell_et_al_1991 would work well after . The authors also identified that specifically at , numerical solutions are prone to divergence blow-ups due to complex shock collisions.

For this test case, we use the FE setup same to all other test cases without changing any parameter. All the shocks are captured locally and accurately by the RV method. This could be a motivation of using entropy viscosity methods, or in particular RV methods, for shock-capturing purposes. Under the same number of degrees of freedom, the solution reveals more structural details of the shocks and turbulence. Similar to the reference (Guillet_2019) using a third-order DG scheme, we capture a twisting density upsurge in the centre of the domain using Lagrange elements.

We discuss the effect of different divergence cleaning techniques in the finite element settings in the remaining of this section.

#### Numerical comparison of the divergence cleaning techniques

We show the absolute divergence error of the Orszag-Tang test as a function of time using different divergence cleaning techniques in Figure 5. The simulation uses elements on a right triangulation mesh. We note that if a sufficiently fine mesh was used, the reference simulation without divergence cleaning would break at some time around due to large divergence error. The figure shows that our method works well with some of the most well-known divergence cleaning techniques in MHD. The divergence error remains small with the use of the projection method, pseudo time-stepping method, and hyperbolic method for divergence cleaning. Among the methods, the pseudo time-stepping method with 10 pseudo iterations performs the best. However, the method is the most expensive because 10 pseudo iterations needs to be computed in each time step. The projection method is more expensive than the hyperbolic cleaning method because a Poisson equation needs to be solved in each time step. The obtained behavior of the divergence error agrees with the reported results, e.g., Tricco_2016; Derigs_2018. We hereby conclude that the investigated divergence cleaning techniques incorporates well with our shock-capturing method, in which the divergence error remains controlled over long time simulation.

### 5.4 3D Orszag-Tang problem, (Helzel_2011)

The 3D Orszag-Tang problem is extended from the 2D Orszag-Tang problem in Section 5.3. Essentially, density and pressure are identical to the 2D problem, independent of the -coordinate. The velocity and magnetic field are also identical in the - and - directions, but zero in the -direction,

 (ρ0,\bu0,p0,\bB0)=(2536π,(−sin(2πy),sin(2πx),0),512π,(−sin(2πy)√4π,sin(4πx)√4π,0)).

We then add a small perturbation to the velocity, similar to Helzel_2011,

 δ\bu=(−ϵpsin(2πz)sin(2πy),ϵpsin(2πz)sin(2πx),ϵpsin(2πz)),

where the perturbation parameter is a small real number. We set same to (Helzel_2011). The initial solution can be written as

 (ρ(t),\bu(t),p(t),\bB(t))=(ρ0,\bu0+δ\bu,p0,\bB0+δ\bB).

The adiabatic constant is kept same to the 2D problem, . All boundaries are periodic. We use the projection method to clean the divergence. The CFL number is 0.3. A solution of the 3D Orszag-Tang problem is shown in Figure 6. We stop the simulation at to make the density plot in Figure 6(a) directly comparable with (Helzel_2011). The density cut presents the highly recognizable 2D Orszag-Tang solution. Along the -axis, the solution agrees well with the reported results in (Helzel_2011). A quarter of the solution is hidden to reveal some of the inside structures. We see that the RV method can capture the shocks efficiently in 3D with the same finite element setup in 1D and 2D. We omit the discussion of divergence error in 3D as the purpose of this paper is to develop accurate shock-capturing techniques for finite element approximations.

### 5.5 MHD Rotor, Balsara_et_al_1999

The final benchmark we present is the rotor problem, first introduced in Balsara_et_al_1999. This problem is referred to as the “first rotor problem” in Toth_2000. The initial solution is a disc being placed in the central of the domain . The ambient solution is a stationary medium at ten times lesser density along a homogeneous magnetic field,

 (ρ,\bu,p,\bB)=(1,(0,0)T,1,(5√4π,0)).

Let be the radius of the disc centered at the domain central. The disc is defined as