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 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 spacial dimensions. Sharp shocks and discontinuity resolutions are obtained.


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.1 arXiv:2112.08885v1[math.NA] 16 Dec 2021 Finite element approximation of MHD 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 B 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 [1,2].Existing numerical methods for compressible Euler equations are hard to extend to approximate the MHD system.If one reason is solving d-additional equations, where d is the space dimension, other two main reasons are (i) developing high order robust shockcapturing techniques, (ii) satisfying the divergence-free constraint for the magnetic field, i.e., ∇• B = 0. 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 [3][4][5][6][7][8][9] 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 e.g., [10][11][12].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 nonunique Riemann solution [13].To overcome this difficulty, numerical methods that do not use Riemann solvers have been developed.For instance central DG schemes reported in [14][15][16] 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 [17,18] 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 timedependent.Only a few works were done to simulate resistive MHD using GLS stabilization, see for instance [19][20][21] 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, [22][23][24][25] 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 [17].In [26], we proved that the FE approximation of residualbased 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 e.g., [27][28][29][30][31].For other approaches to constructing the nonlinear artificial viscosity methods, we refer to the work of [32][33][34], 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 B into a divergence-free space [35] and the hyperbolic cleaning [36].We should emphasize that positivity and invariant domain preservation of the method as in [37] 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.

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 Ω ∈ R d , d = 1, 2, 3 with boundary Γ.We denote T h a finite partition of Ω into disjoint elements K 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 K constitutes Ω, where Ω denotes the closure of Ω.
We seek finite element solution in Lagrange polynomial spaces where P k is a space of polynomials of at most k degrees.Let Φ i be the kdegree Lagrange basis function in Q h which takes value 1 at node i, and 0 at any other nodes.The functions {Φ i } N h i=1 form a basis for Q h , where N h is the total number of nodes in Q h .We define I(i) the set of all nodal points contained within the support of Φ i .Furthermore, we use the following Hilbert space inner products, We will use a mesh-function h h (x) ∈ Q h which we compute using the following projection: find where h K is the circumradius of the element K.

Governing equations of MHD
Let us denote the time interval [0, t ], where t is the final time.For all (x, t) ∈ D := Ω × [0, t ], we seek solution of the MHD equations , where ρ is the mass density, m is the momentum, E is the total energy, B is the magnetic field.The ideal MHD equations in conservative form with appropriate boundary conditions are defined as follows: where the nonlinear tensor fluxes F E (U), F B (U) are defined as with u := m/ρ being the velocity field, and β being the Maxwell stress tensor: where I denotes the identity matrix of size d × d.The combination of hydrodynamics and electromagnetism nature of the MHD system can be seen from the MHD equations in the form of (3.1), through the separated fluxes F E (U) and F B (U).The thermodynamic pressure p is computed from the equation of state where e is the internal energy and γ is the adiabatic gas constant.The internal energy is computed as In addition, the magnetic field should satisfy the following divergence free constraint:

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 e.g., [4,38,39]: let e ∈ R d be a direction vector, a := γp/ρ be the speed of sound, b := B•e/ρ 1 2 and .
The eigenvalues from the smallest to the largest are where the subscripts f , s and a correspond to fast and slow magnetosonic waves, and the Alfén waves.The characteristics λ 4,5 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 in the construction of the nonlinear viscosity term.

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

Elliptic correction/Projection method
In each time step, instead of using the magnetic solution B from the system, B the projection of B onto a magnetic divergence-free space is used.[35] proposed using where Ψ p is the solution to the Poisson equation ∇ 2 Ψ p − ∇• B = 0.This yields ∇• B = 0 as the non-physical part ∇Ψ p is removed from the numerical solution B. After each correction in the magnetic field B, the dependent variables: pressure p, temperature T , energy E and the entropy functions are updated accordingly to ensure consistency of the discrete solution.Contrary to the seemingly ad-hoc nature of this method, [41] shows that the projection method preserves conservative properties and accuracy of the base scheme.
Pseudo time-stepping.Numerically minimizing |∇• B | by the projection method is considered computationally expensive.A more affordable alternative, proposed by [42], is solving the following time-dependent equation for Ψ p , where the variable t is a separate time dimension.At equilibrium, the solution to (3.6) is the one to the Poisson equation.A semi-discretization of (3.6) is given by where τ is the pseudo step size.After each pseudo time step, the magnetic solution B is corrected B +1 = B −∇Ψ +1 p .The divergence ∇• B 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.

Hyperbolic correction
An auxiliary scalar variable Ψ l 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 where c r = 0.3 in 2D or c r = 1.0 in 3D [40], and h h (x) denotes the mesh function computed in (2.2).The coefficient c h is set to be the maximum wave speed In addition, a source term is added to the right hand side of (3.1), where

Nonlinear viscosity method for MHD
The aim of this paper is to solve the inviscid MHD system (3.1) using continuous finite element approximations.We use the following two vector valued spaces: where Q h is defined in (2.1).Then, we formulate the finite element approximation of the MHD system (3.1) as follows: find Here the inner products are computed exactly using appropriate numerical quadrature rules.The finite element approximation of the flux terms in (4.2) 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.

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 ∇• F V (U) into the MHD system, where F V (U) denotes a viscous flux tensor.
As a heuristic argument to construct F V (U), we start from the viscous flux of the resistive MHD equations, see e.g., [43], where the temperature T h is a function in Q h .The viscous shear stress tensor τ h is defined as where µ h ≥ 0 is the artificial dynamic viscosity, λ h is the bulk viscosity of undetermined sign such that λ h +2µ h > 0, κ h ≥ 0 is the artificial heat conduction, η h ≥ 0 is the artificial electric resistivity of the medium.Note that the electric resistivity η h has the same dimension as the dynamic viscosity coefficient µ h .We point out that the temperature and velocity field are computed at the nodal points: By considering zero magnetic field, i.e., B ≡ 0 and thus F B (U) ≡ 0, we obtain from (3.1) the compressible Euler equations.For the compressible Euler equations, a frequent choice of viscous regularization is the Navier-Stokes flux, see e.g., [25], which coincides with F 0 V (U h ) when B ≡ 0. The Navier-Stokes flux does not add any viscosity term to the mass equation.Therefore, positivity of density can be violated [25,44].In [44], 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 [25].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 ν h 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 where n 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 ν h , κ h , and η h , all are from the finite element space Q h , and make them proportional to a nonlinear functional such that enough viscosity is added to the system in terms of stability and accuracy.

First-order viscosity
At each node i ∈ {1, 2, . . ., N h }, we compute the first order viscosity at time t = t n as µ L,n i = C max h i λ n max,i , where λ n max is the largest eigenvalue at time t = t n , given by (3.5), and C max is a positive parameter.A typical range of C max is reportedly [0.15, 0.5], see [27].In one space dimension, it can shown that C max = 0.5 leads to the first order Lax-Friedrichs scheme.We use C max = 0.5 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.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 [45], 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.

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, i.e., ∼ O(h k+1 h ), and has a jump of size ∼ O(h −1 h ) 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 w : R d → R and in particular the case w ∈ Q h .We propose dividing w by the following function to normalize w n(w) i = S(w) − α max j∈I(i) where S indicates the global maximum variance S(w) = w − mean(w) L ∞ (Ω) , and the global constant α is calculated as which acts as a unit normalizer of local jumps to the global variance.As a technical detail, notice that when w is a constant uniform function, we have max Ω w − min Ω w = max j∈I(i) w j − min j∈I(i) w j = 0 for all i ∈ {1, 2, . . ., N h }.In that case, n(w) reduces to (1 − C l ) S(w) which although being zero, we do not suffer from division by zero in the calculation of n(w).In an earlier development, the normalization function n(w) in [30] can be obtained by setting α = 1 in (4.6).We observe that the new formulation (4.6) deals better with mixed-signed solutions since the two factors in n(w) are of the same unit.To demonstrate this argument, we rewrite the expression (4.6) as The parameter C l ∈ (0, 1) allows fine tuning of how sensitive this normalization affects the local discontinuities.Observe that the ratio (max j∈I(i) w j − min j∈I(i) w j )/(max Ω w−min Ω w) is in between the interval [0, 1], and is amplified at local solution jumps.If the input function w holds both positive and negative parts, and if α = 1, it is possible that n(w) i = 0 at the regions where the local jumps incidentally matches the global variance.However, once the local jump bypasses the global variance, n(w) can become uncontrollably big, which opposes to our aim that n(w) should be small at such circumstances.On the other hand, the formula (4.7) does not suffer this issue because the value within the bracket in (4.7) is strictly positive.How small n(w) can become is tuned with the parameter C l .For example, if C l = 1, n(w) can become zero at the largest jump of w.Lowering C l from 1 towards 0 would make n(w) at jumps of w larger and larger but not exceeding S(w).In smooth regions, n(w) turns into S(w), which is the classical normalization function of our residual-based viscosity technique, see [26].
Next, we construct a normalized finite element residual of the MHD equations as  n(B n h,1 ), . . ., n(B n h,d ), as defined above, normalize the residual units and magnitude.It is also possible that the denominators in (4.8) 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 R(q h )/n(q h ), where q h is a scalar function, we use R(q n h ) n(q n h ) n(q n h ) 2 + .We choose to be machine epsilon, i.e. approximately 2.2e −16 .

For each component q
where f (q n h ) corresponds to the nonlinear MHD flux in (4.2) for each component q n h , and BDF(q h ) n is a backward difference approximation of the time derivative ∂ t q h at time t = t n .Note that the second-order BDF scheme is sufficient in getting 4th order accurate scheme.
For P 1 , P 3 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, where m i = Ω ϕ i dx > 0 for P 1 , P 3 Lagrange polynomials.Finally, we are now ready to construct the final artificial viscosity on each node i at time t = t n : where [27].For the numerical simulations in this paper, we use C R = 1.In the rest of this paper, we refer to using (4.10) in the formulation of the viscous regularized equation (4.5) as the "RV method".

Adaptive time-stepping
Once the MHD system is discretized in space using the continuous finite element method we obtain the system of ODEs (4.5).Let us denote this system as M∂ t U h (t) = F(U h (t), µ h (t)), where M ∈ R (2d+2)N h ×(2d+2)N h is the mass matrix, F(U h (t), µ h (t)) is the right-hand-side function of the system which depends on the solution U h (t) and the viscosity vector is µ h (t) := (ν h (t), µ h (t), κ h (t), η h (t)).
Next, discretize this system in time using explicit r-stage Runge-Kutta methods: . ., r are coefficients obtained from the Butcher tableau, and the stage variables K i are computed as follows: for the given solution U n h and the viscosity vector µ n h := µ h (t n ) at time level t n , set W 0 := U n h , then let W l be the solution at the l-th stage of the Runge-Kutta method, then compute K l by solving the following system: for all l = 0, . . ., r.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 τ n is computed adaptively using a CFL condition We refer the reader to [25], where advantages of the strong implementation of the slip boundary condition versus its weak counterpart are discussed, and we refer to [37], where the correction technique is discussed for explicit schemes in more detail.

Summary of the algorithm
To summarize, our stabilized finite element solution of the MHD system (3.1) can be obtained by the following steps.In each time step, the solution is advanced as 1.Calculate the residual using (4.9), combine the low order and high order viscosity to calculate the amount of artificial viscosity using (4.10).2. Solve the main linear system (4.5)where the time derivative is discretized using a Runge-Kutta method (4.11) of order at least k + 1. 3. If the projection method or the pseudo time-stepping method is used for cleaning divergence of B h , 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.

Compute an approximation of the maximum local velocity using (3.5), and
identify the next time step size using (4.12).

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 C max = 0.5, C R = 1.0,C l = 0.4.We observe that the parameter C l is not sensitive as in most circumstances, apart from enhanced smoothness at sharp shocks, any choice of C l ∈ (0, 1) would give satisfactory outcomes.We have used the classical Runge-Kutta 4 for all the tests with time step (4.12),where the CFL number is chosen to be 0.3.In the viscous flux (4.4), the viscosity coefficients ν h , κ h , and η h are chosen to be µ h given by (4.10).

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.where

Smooth vortex problem, [46]
the vortex radius is r = r 2 1 + r 2 2 , (r 1 , r 2 ) = (x, y) − u 0 t, and the vortex strength is µ = 5.389489439.The adiabatic constant is γ = 5  3 .The errors measured at final time t = 0.05 using P 1 , P 2 , P 3 elements are shown in Table 1.In Table 1, we show behavior of both the Galerkin solutions and the RV solutions.For the P 2 solution, the obtained error is lower than that of the P 1 solution under the same number of degrees of freedom.However, the convergence rate of the P 2 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., [25].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.
Fig. 1: Accuracy test on the smooth vortex problem: convergence history for velocity u h against the exact solution.The RV method does not destroy the accuracy of high-order solutions.

Brio-Wu MHD shock tube problem, [3]
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 [3].This is the one-dimensional Riemann problem in the domain Ω = [0, 1].The problem is challenging because the solution contains waves of different types, typical for an ideal MHD.
The result in Figure 2 shows that our proposed method can capture the wave structures efficiently and accurately.We also see that P 3 solution is sharper than P 1 for the same number of degrees of freedom.
5.3 2D Orszag-Tang problem, [48] In this section, we consider the well-known benchmark Orszag Fig. 2: Comparison of P 1 and P 3 RV-solutions to the Brio-Wu problem.The reference solution is obtained from the Athena code [47] with 10000 grid points.
Under the same number of degrees of freedom, the P 3 solution captures the compound structure more accurately than the P 1 solution.
The adiabatic constant for this problem is γ = 5 3 .The periodic boundary condition is used in both directions.For divergence cleaning, we have used the projection method.The density solution ρ h and the artificial viscosity µ h at time t = 0.5 and t = 1.0 are shown in Figure 3 using P 1 elements, and in Figure 4 using P 3 elements.At t = 1.0, the structure of the Orszag-Tang solution is considered turbulence [40].
Despite being a classical benchmark in studies of numerical schemes for MHD, the evolution of the Orszag-Tang solution after t = 0.5 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 [41], the solution at t = 3.14 approximately corresponds to our solution at t = 0.5.The authors in [49] 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 t = 0.5.In fact, it is evident that the Orszag-Tang solution up to t = 0.5 is relatively straightforward to be obtained by numerical schemes with shock-capturing capability, e.g., [41,49,50].However, long time integration would discriminate different discretizations, as remarked by [50].Using DG method, the authors in [50] reported that not all choices of the slope limiters and discretizations of the Powell term [4] would work well after t = 0.5.The authors also identified that specifically at 0.75 ≤ t ≤ 0.85, 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 P 3 solution reveals more structural details of the shocks and turbulence.Similar to the reference [50] using a third-order DG scheme, we capture a twisting density upsurge in the centre of the domain using P 3 Lagrange elements.
We discuss the effect of different divergence cleaning techniques in the finite element settings in the remaining of this section.We note that if a sufficiently fine mesh was used, the reference simulation without divergence cleaning would break at some time around t = 0.7 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., [40,51].We hereby conclude that the investigated divergence cleaning techniques incorporates well with our shockcapturing method, in which the divergence error remains controlled over long time simulation.The red, yellow, and purple lines respectively corresponds to the simulations using the projection method, the pseudo time-stepping method, and the hyperbolic cleaning method to clean the divergence.

3D Orszag-Tang problem, [52]
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 z-coordinate.The velocity and magnetic field are also identical in the x-and y-directions, but zero in the z-direction, We then add a small perturbation to the velocity, similar to [52], δu = (− p sin(2πz) sin(2πy), p sin(2πz) sin(2πx), p sin(2πz)), where the perturbation parameter p is a small real number.We set p = 0.2 same to [52].The initial solution can be written as (ρ(t), u(t), p(t), B(t)) = (ρ 0 , u 0 + δu, p 0 , B 0 + δB).
The adiabatic constant is kept same to the 2D problem, γ = 5 3 .All boundaries are periodic.We use the projection method to clean the divergence.The CFL number is 0.3.A P 3 solution of the 3D Orszag-Tang problem is shown in Figure 6.We stop the simulation at t = 0.5 to make the density plot in Figure 6(a) directly comparable with [52].The density cut z = 0 presents the highly recognizable 2D Orszag-Tang solution.Along the z-axis, the solution agrees well with the reported results in [52].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.The final benchmark we present is the rotor problem, first introduced in [6].This problem is referred to as the "first rotor problem" in [41].The initial solution is a disc being placed in the central of the domain Ω = [0, 1] × [0, 1].The ambient solution is a stationary medium at ten times lesser density along Finite element approximation of MHD a homogeneous magnetic field, (ρ, u, p, B) = 1, (0, 0) T , 1, 5 √ 4π , 0 .
Let r be the radius of the disc centered at the domain central.The disc is defined as (ρ, u, p) = 10, u0 r0 (0.5 − y), u0 r0 (x − 0.5) if r < r 0 1 + 9f, f u0 r (0.5 − y), f u0 r (x − 0.5) if r 0 ≤ r < r 1 , where f = (r 1 −r)/(r 1 −r 0 ), the inner radius r 0 = 0.1 and the outer radius r 1 = 0.115.The adiabatic constant is γ = 1.4.Due to imbalance in the centrifugal force, the disc is rotated as it proceeds in time, generating a twisting moment in the magnetic field.The difficulties of this test include: (i) resolving the torsional Alfvén waves, and (ii) during the rotation, the pressure can get close to zero which is challenging for numerical methods.[41] reported that many schemes fail to solve this problem owing to negative pressure.Our density and pressure solution at t = 0.15 are plotted in Figure 7.The projection method is used for cleaning the divergence.It can be seen that the signature circle contour lines within the middle of Figure 7(c) showing the rotational evolution, see e.g., [41], are finely captured.

Conclusion
The main purpose of this paper has been to numerically investigate a new viscous regularization of the ideal MHD equations.The proposed viscous regularization is inspired by the work of [44] on the Euler equations combining with the elliptic terms of the resistive MHD equations.The discretization consists of a high-order residual-based viscosity finite element method in space and the classical explicit Runge-Kuta methods in time.The results have shown that the proposed method can capture different structures of MHD discontinuities while achieving high-order accuracy, experimented up to fourth-order, in shock-free regions.Future extensions of this paper may concern continuous analysis of the proposed viscous flux regarding e.g., positivity-preserving properties, Galilean invariance, and entropy principles.Discrete invariant-domain preserving schemes can be built upon such continuous properties using, for instance, the viscous operator by [53].

Statements and Declarations
Funding.This research is funded by Uppsala University, Sweden.
Competing interests.The authors have no conflicts of interest to disclose.
Data availability statement.Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

.12) 4 . 5
Boundary conditionsIn 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 t n+1 , U n+1 h , we impose the Dirichlet boundary condition by setting the values of the solution on the boundary nodes N j by its boundary data.The slip or an impermeability boundary condition requires u n+1 h •n = 0 on the boundary in question.This condition is equivalent to m n+1 h •n = 0.Then, the boundary value of momentum at the boundary node N j , m n+1 h (N j ) is replaced by m n+1 h (N j ) − m n+1 h (N j )•n.

Fig. 5 :
Absolute divergence error Ω |∇• B h | dx as a function of time in the 2D Orszag-Tang test.The reference simulation without divergence cleaning corresponds to the blue dot-dashed line.

Fig. 7 :
Fig. 7: P 3 solution of the Rotor problem, 300 × 300 nodes Φ i be the low, high and dynamic viscosity coefficients at time t = t n , and h h

Table 1 :
Smooth vortex problem.L 1 , L 2 error and convergence rates of velocity u and magnetic field B. P 1 , P 2 , P 3 solutions.t = 0.05.