Jonatan Núñez-De La Rosa

Claus-Dieter Munz
email: munz@iag.uni-stuttgart.de

Hybrid DG/FV schemes for magnetohydrodynamics and relativistic hydrodynamics

Keywords: discontinuous Galerkin spectral element methods, magnetohydrodynamics, shock capturing, divergence cleaning, high-order methods

This paper presents a high order hybrid discontinuous Galerkin/finite volume scheme for solving the equations of the magnetohydrodynamics (MHD) and of the relativistic hydrodynamics (SRHD) on quadrilateral meshes. In this approach, for the spatial discretization, an arbitrary high order discontinuous Galerkin spectral element (DG) method is combined with a finite volume (FV) scheme in order to simulate complex flow problems involving strong shocks. Regarding the time discretization, a fourth order strong stability preserving Runge-Kutta method is used. In the proposed hybrid scheme, a shock indicator is computed at the beginning of each Runge-Kutta stage in order to flag those elements containing shock waves or discontinuities. Subsequently, the DG solution in these troubled elements and in the current time step is projected onto a subdomain composed of finite volume subcells. Right after, the DG operator is applied to those unflagged elements, which, in principle, are oscillation-free, meanwhile the troubled elements are evolved with a robust second/third order FV operator. With this approach we are able to numerically simulate very challenging problems in the context of MHD and SRHD in one, and two space dimensions and with very high order polynomials. We make convergence tests and show a comprehensive one-and two dimensional testbench for both equation systems, focusing in problems with strong shocks. The presented hybrid approach shows that numerical schemes of very high order of accuracy are able to simulate these complex flow problems in an efficient and robust manner.

Introduction

High order numerical methods for solving partial differential equations are becoming the new standard in computational science and engineering because of the advantages provided by this family of numerical schemes. Among the high-order methods, the discontinuous Galerkin schemes are the most renowned candidates and are playing a central role because of all the attractive properties for solving the equations of fluid dynamics on complex geometries, and in general, nonlinear hyperbolic systems of conservation laws [1][START_REF] Canuto | Spectral Methods: Fundamentals in Single Domains[END_REF][START_REF] Kopriva | Implementing Spectral Methods for Partial Differential Equations[END_REF]. In these schemes, high-order accuracy is obtained through the approximation of the solution with a higher degree polynomial. The method has the capacity to handle complicated geometries, and due to its locality, it is highly parallelizable. In 1971, the discontinuous Galerkin method was originally introduced by Nitsche [START_REF] Nitsche | [END_REF] for the approximation of elliptic equations. In 1973, Reed and Hill [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF] proposed the discontinuous Galerkin method for solving the steady-state neutron transport equation. The method was recovered by Cockburn and Shu [START_REF] Cockburn | [END_REF] to solve time dependent non-linear conservation laws. In the subsequent years, in a series of papers, the general framework of the discontinuous Galerkin for nonlinear systems of conservation laws, even in multiple space dimensions was established [7][8][9].

High-order methods are not as robust as low-order methods. Because of hyperbolic conservation laws allow discontinuities in their solution, these will produce spurious oscillations in the case the solution is computed with a high-order scheme. This phenomenon is called the Gibbs phenomenon [10,11]. As a consequence, the oscillations will create unphysical states, and destroy conservation. Shock capturing methods have to be constructed and adapted to a high order numerical scheme in order to overcome this difficulty. In the case of finite difference and finite volume methods, the most used shock capturing strategies are those that employ slope/flux limiters [12,13], and an adaptive reconstruction operator [14,15]. For discontinuous Galerkin schemes, several approaches have been developed, for instance the generalized slope limiter [START_REF] Cockburn | [END_REF], the moment limiter [16,17], and the artificial diffusion [START_REF] Persson | Proceedings of the AIAA 44th Aerospace Sciences Meeting and Exhibit, AIAA-2006-0112[END_REF][START_REF] Casoni | [END_REF]. The family of WENO reconstruction methods used as limiters in the discontinuous Galerkin methods has been well studied in the last years [20][21][22]. The main drawback of the like-WENO limiters is that they have a huge computational overhead when high-order approximations are used.

A very promising approach that allows to solve conservation laws even when shocks/discontinuities take place consists in hybridizing high order numerical schemes with robust low order methods, and in this way we take advantage of the best features of each of them. For instance, a hybrid discontinuous Galerkin/finite volume scheme is an excellent approach to tackle complex flows with shocks [START_REF] Sonntag | Finite volumes for complex applications VII-Elliptic, parabolic and hyperbolic problems[END_REF][24][START_REF] Núñez-De La | High-Order Methods for Computational Astrophysics[END_REF][START_REF] Sonntag | [END_REF]. The main idea behind these schemes is to use a high order discontinuous Galerkin scheme in flow regions where shocks and discontinuities are absent, and to employ a robust finite volume method in flow regions with such shocks/discontinuities. A more complete strategy consists in efficiently detecting troubled regions (with shocks, discontinuities, violation of positivity of some quantities of interest like pressure and density) and then applying an adaptive mesh refinement algorithm in those regions along with a robust finite volume scheme. Furthermore, this finite volume scheme can also have an adaptive reconstruction operator that lower its order in the extreme case it cannot obtain a physically correct solution. As last option in this adaptation, the algorithm will use a first order Godunov scheme. In this work we follow a simple, efficient and robust strategy based on the works of [24,[START_REF] Sonntag | [END_REF]27], and it is essentially a combination of a high-order DG method with a robust FV-WENO scheme on quadrilateral meshes. In the proposed hybrid scheme, an oscillations indicator is used to mark those elements containing discontinuities. The DG solution in these troubled elements and in the current time step is projected in to a subdomain made of finite-volume subcells, which are then evolved with a robust FV method. The cells without oscillations are evolved with a high-order DG method.

Several high resolution shock-capturing and high order numerical formulations have been constructed in the last years to solve the MHD and SRHD equations. For the MHD equations the finite difference methods [28,29], and the finite volume methods [30][31][32][33] have been the most widely used. For the SRHD equations several high resolution shock capturing schemes based on finite difference and finite volume schemes have been successfully developed [34][35][36], and also high order schemes [37][38][39][40]. The hybrid DG/FV scheme presented in this work will be used to solve these equations in an efficient and robust manner, and using very high order polynomials. In the case of the MHD equations, besides the shock capturing strategy, and additional stabilization procedure has to be developed in order to maintain the solenoidal constraint [41,[START_REF] Powell | An approximate Riemann solver for magnetohydrodynamics[END_REF]. The constraint ∇ • B = 0 is not satisfied by the numerical scheme, even for high order methods. The most used approaches to keep it are the divergence cleaning and the constrained transport algorithms. The main idea of the divergence cleaning is to remove the numerical errors produced in the computation of ∇•B. The Hodge projection [41], the 8-wave formulation [START_REF] Powell | An approximate Riemann solver for magnetohydrodynamics[END_REF], and the generalized Lagrange multiplier (GLM) approach [START_REF] Dedner | [END_REF] belong to the family of divergence cleaning. The constrained transport schemes preserve the solenoidal constraint of the magnetic field with machine accuracy, but introducing a staggered mesh [44]. In order to keep as simple as possible the numerical method, in this work we will use the GLM approach for controlling the solenoidal constraint of the magnetic field. In this approach a new scalar field ψ is introduced and couples the solenoidal constraint with the evolution equation of the magnetic field. The divergence errors produced during the numerical simulation are propagated to the boundary of the computational domain. This scheme is rather simple to implement, and the produced errors in the computation are small enough to guarantee conservation of all physical variables. In our algorithm we use the mixed hyperbolic/parabolic divergence cleaning as it is reported in [45].

The structure of this paper is as follows: in section 2 we briefly describe the MHD and SRHD equations written in conservation form. The GLM divergence cleaning procedure is also taken into consideration in this section. In section 3 the high order hybrid DG/FV method for conservation laws is discussed in detail. In section 4 several test problems are shown and analyzed. This includes convergence tests, shock tube problems and two-dimensional flow configurations with shocks for the MHD and SRHD equations. In section 5 the conclusions of this work are presented.

Governing equations

MHD equations

Magnetohydrodynamics is a macroscopic theory concerned with the study of magnetic fields interacting with electrically conducting fluids, ranging from laboratory to space plasmas. The ideal MHD equations are usually derived as a combination of the Euler equations of the hydrodynamics and the Ampère's law, but neglecting the displacement current. As a system of conservation laws, the ideal MHD equations consist of the conservation of mass, the conservation of momentum, the conservation of energy, and the induction equations:

∂ρ ∂t + ∇ • S = 0, (1a)
∂S ∂t + ∇ • S ⊗ v + P = 0, (1b)
∂E ∂t + ∇ • Ev + P • v = 0, (1c)
∂B ∂t + ∇ • B ⊗ v -v ⊗ B = 0. ( 1d
)
The conservative variables are the mass density ρ, the momentum density S, the energy density E and the magnetic field B. The magnetic field is subject to the solenoidal constraint, that is

∇ • B = 0. (2)
The pressure tensor in equations (1b) and (1c) includes the hydrodynamic and the magnetic pressure

P = p + 1 2 |B| 2 I -B ⊗ B. (3)
Finally, an equation of state (EOS) is used to close the system. In this work we only consider the ideal gas EOS with adiabatic exponent γ

p = ρ (γ -1), (4)
where is the specific internal energy. This EOS can be written as a function of the total energy E, and the primitive variables, ρ, v, and

B p = (γ -1) E - 1 2 ρ |v| 2 - 1 2 |B| 2 .
(5)

Divergence cleaning

The solenoidal constraint given in equation ( 2) is a condition that the numerical method has to maintain in order to obtain physically correct solutions. In this work we apply the GLM divergence cleaning because this is a simple and reliable method [START_REF] Dedner | [END_REF]46]. The mixed hyperbolic/parabolic divergence cleaning is used in this work [START_REF] Dedner | [END_REF]. This scheme introduces a scalar field ψ in order to couple the solenoidal constraint (2) with the induction equation (1d). As result, the following equations are obtained

∂B ∂t + ∇ • B ⊗ v -v ⊗ B + ψI = 0, (6a)
∂ψ ∂t + ∇ • c 2 h B = - c 2 h c 2 p ψ. (6b)
The MHD-GLM system consists in equations (1a), (1b), (1c), (6a), and (6b). Observe that this system is conservative (except for the equation for the unphysical scalar field ψ).

Eigenvalues of the MHD-GLM system

We solve the MHD equations with an explicit numerical scheme. Therefore, the computation of the eigenvalues is required in order to compute the current time step (see section 3.3). From the Jacobian matrices for the MHD-GLM equations we can get the eigenvalues of the MHD-GLM system. This matrices, denoted by A α , are defined by

A α = ∂f α (u) ∂u , (α = x, y, z), (7)
with u the state vector of conservative variables (along the scalar field ψ) and f ≡ [f , g, h] the tensor of physical fluxes. The fluxes f , g, h are respectively the flux in the x-direction, in the y-direction, and in the z-direction. The eigenvalues of the Jacobian matrices in the x-direction, A x , in nondecreasing order, are

λ 1 = -c h , λ 2 = v x -c f , λ 3 = v x -c a , λ 4 = v x -c s , λ 5 = v x , λ 6 = v x + c s , λ 7 = v x + c a , λ 8 = v x + c f , λ 9 = c h , (8)
where c f , c a , and c s are the fast, Alfvén, and slow characteristic speeds. The cases for the y and z directions can be easily derived from symmetry. In equation ( 8) we can see that there are three MHD waves families associated with the characteristic speeds and an entropy mode. The three characteristic speeds are given by

c a = |b x | , c s = 1 2 a 2 + b 2 - a 2 + b 2 2 -4a 2 b 2 x , c f = 1 2 a 2 + b 2 + a 2 + b 2 2 -4a 2 b 2 x , (9)
where the following abbreviations were used

a 2 = γp ρ , b 2 = |B| 2 ρ , b 2 x = B 2 x ρ . (10)
Two new eigenvalues appear in the MHD-GLM equations, namely λ 1,9 = ±c h . These eigenmodes transport the divergence errors to the boundaries at speed c h , and are also damped at a rate c 2 h /c 2 p . Following [START_REF] Dedner | [END_REF], the value of the constant c p is chosen after setting c r ≡ c 2 p /c h = 0.18, and the magnitude of the wave speeds c h is set to the maximum allowed speed in the ideal MHD equations (without divergence cleaning)

c h = max α max |λ 2,α |, |λ 8,α | , α = x, y, z. (11)

SRHD equations

Taking into account the effects of the special relativity, the SRHD equations are a generalization of the Euler equations. A relativistic fluid may have as main features at least one of the following properties: (1) the fluid velocity is close to the speed of light, and (2) the internal energy of the fluid is greater or approximate to its rest mass density [START_REF] Rezzolla | Relativistic Hydrodynamics[END_REF]. The SRHD equations, written as a system of conservation laws [START_REF] Mignone | [END_REF]49], and taking the speed of light c = 1, are given by

∂D ∂t + ∇ • (Dv) = 0, (12a)
∂S ∂t + ∇ • (S ⊗ v + pI) = 0, (12b)
∂E ∂t + ∇ • S = 0. ( 12c
)
The conserved quantities D, S and E are the mass density, the momentum density, and the total energy density, respectively. They are all measured in the laboratory frame, and are related to quantities in the local rest frame of the fluid, the so-called primitive variables (the rest-mass density ρ, the gas pressure p, and the specific enthalpy h), through

D = ρΓ, (13a)
S = ρhΓ 2 v, (13b)
E = ρhΓ 2 -p, (13c)
where the Γ is the Lorentz factor defined by

Γ = 1 √ 1 -v 2 . ( 14
)
The SRHD system ( 12) is closed with an EOS: h = h(p, ρ), also written as p = p(ρ, ε). The SRHD system is hyperbolic for causal EOS [START_REF] Rezzolla | Relativistic Hydrodynamics[END_REF][START_REF] Anile | Relativistic Fluids and Magnetofluids[END_REF], i.e., for those where the local sound speed satisfies c s < 1, where c s is defined by

hc 2 s = -ρ ∂h ∂ρ ρ ∂h ∂p -1 -1 . (15)
In this work, we employ an ideal gas EOS,

p = (γ -1)(ε -ρ), or h = 1 + γ γ -1 p ρ . (16)
Here γ = c p /c v is the ratio of specific heats and ε is the sum of the internal and rest-mass energy densities in the local frame and is related to the specific enthalpy as

h = ε + p ρ . ( 17
)
With this choice of EOS, the speed of sound is given by

c 2 s = γp ρh . (18)

Eigenvalues of the SRHD system

The eigenvalues of the SRHD equations are also required for the explicit numerical scheme for computing the time step. The Jacobian matrices A α for the SRHD equations are defined by

A α = ∂f α (u) ∂u , (α = x, y, z), ( 19
)
where u is the state vector of conservative variables, and f = [f , g, h] the tensor of physical fluxes. Assuming an ideal gas EOS, the eigenvalues of the matrix A x are given by

λ 1 = 1 1 -v 2 c 2 s v x (1 -c 2 s ) -c s (1 -v 2 ) 1 -v x v x -(v 2 -v x v x )c 2 s , (20)
λ 2 =v x , λ 3 = v x , λ 4 = v x , (21)
λ 5 = 1 1 -v 2 c 2 s v x (1 -c 2 s ) + c s (1 -v 2 ) 1 -v x v x -(v 2 -v x v x )c 2 s . ( 22
)
The cases for the matrices A y and A z easily follows from symmetry. The eigenvalues of A x represent the five characteristic speeds associated with two sound wave modes (λ 1 and λ 5 ) and three entropy modes (λ 2 , λ 3 , λ 4 ). The eigenvalues satisfy the following relation:

λ 1 < λ 2 = λ 3 = λ 4 < λ 5 .

Conservative to primitive variables conversion

In the SRHD system, the relation between the primitive and conservative variables is quite complex. In this work, two strategies for expressing the primitive variables in terms of the conservative ones are followed: the first one, by solving iteratively an equation for the pressure and then determining the other variables from this [START_REF] Mignone | [END_REF]. The second one, which is more robust, consists in solving a quartic equation on the velocity v (as it is outlined in [START_REF] Schneider | [END_REF]52], and [49]). We combine both strategies in our algorithm, using the iterative scheme as base scheme, and when this fails, the robust approach is used. The details of both algorithms can be found in [START_REF] Mignone | [END_REF]49].

Numerical methods

Discontinuous Galerkin methods

Let us consider the following system of conservation laws

∂u ∂t + ∂f (u) ∂x + ∂g(u) ∂y + ∂h(u) ∂z = 0, (23)
also written as

∂u ∂t + ∇ x • f(u) = 0, (24)
where u ∈ R m is the vector of conserved quantities, and f = f , g, h is the tensor of physical fluxes, with f , g, h ∈ R m being the flux functions in the directions x, y, and z, respectively. Next let us consider the test function φ = φ(x) = φ(x, y, z), which is an arbitrary differentiable function defined in a subset Ω k ⊂ Ω ⊂ R 3 . We assume that the subset Ω k corresponds to a cell in the discontinuous Galerkin scheme. The variational formulation consists in projecting the residual of the partial differential equation ( 24) on the test function φ(x) and demanding that this projection is zero

Ω k ∂u ∂t + ∇ x • f(u) φ(x) dx = 0. ( 25
)
The weak formulation of the equation ( 24) is obtained by integrating by parts the flux term in the equation ( 25)

Ω k ∂u ∂t φ(x) dx + ∂Ω k f • n φ(x) ds - Ω k f • ∇φ(x) dx = 0. (26)
So far, no approximation has been done on the vector u or the fluxes f , g, h. Depending in the discretization ansatz on these functions, we have different flavors of the discontinuous Galerkin method. In this work we will focus on the very efficient discontinuous Galerkin spectral element method (DGSEM) developed by Kopriva and Kolias [53], Kopriva [54]. For an outstanding description and derivation of the method and how it is implemented, we recommend the book by Kopriva [START_REF] Kopriva | Implementing Spectral Methods for Partial Differential Equations[END_REF]. In a recent work the DGSEM has been applied to solve the Navier-Stokes equations [55]. In the DGSEM framework, a mapping from the physical domain, where we want the solution, to the reference domain, where actually the equation is solved, is constructed. With the help of this transformation, we map the equation to the reference domain and finally we solve the transformed equation in the reference domain. We stress that in the general formulation of the discontinuous Galerkin, it is usual to solve the original equation in the physical space, where the involved integrals are mapped to the reference element. In our case we map first the equation to the reference element and then we solve it there. Therefore, we apply the discontinuous Galerkin formulation [START_REF] Sonntag | [END_REF] not to the original equation (24), but onto the mapped equation.

Coordinate transformation

Let us start by considering the domain Ω ⊂ R 3 with boundary ∂Ω. We decompose Ω in K non-overlapping hexahedral elements Ω k , which are further mapped onto the reference element [-1, 1] 3 (see figure 1). A point in the physical space is described by x = (x, y, z) = (x 1 , x 2 , x 3 ), and a point in the reference element by ξ = (ξ, η, ζ) = (ξ 1 , ξ 2 , ξ 3 ). The mapping between the physical and the reference element is defined as x = X(ξ). Let us consider a function f (x) = (x, y, z). Under a coordinate transformation, the derivative on the reference space is ∇ ξ f = J∇ x f , where J is the Jacobi matrix of the transformation x = X(ξ). Because we want to map the equations from the physical space to the reference space, in practice, the coordinate transformation we need is the inverse, namely, ξ = X -1 (x), which is represented in matrix form by J -1 .

Any point in the physical space can be represented by using two different basis vectors: The covariant basis vectors e i , for i = 1, 2, 3, and the contravariant basis vectors e i , for i = 1, 2, 3. By definition, the covariant vectors are tangent to the coordinate lines, while the contravariant are normal the coordinate lines. That is, they are defined by the following expressions Vectors can be written either in covariant basis or in contravariant basis. That is possible because these bases form a base of the three-dimensional space. Now, coming back to the inverse of the Jacobi matrix, we use the definition of covariant basis vectors and we then introduce them in the expression for J, yielding

e i = ∂x ∂ξ i , i = 1, 2, 3, (covariant vectors), (27a)
e i = ∇ξ i , i = 1, 2, 3, (contravariant vectors). (27b) ξ η ζ X -1 (x)
J = e 1 , e 2 , e 3 . ( 28
)
Using the general formula for the inverse of a 3 × 3 matrix expressed in terms of the cross product and triple product, we get

J -1 = 1 det(J)           e 2 × e 3 e 3 × e 1 e 1 × e 2           , (29)
where the Jacobian J = det(J) is given by

J = det(J) = e i • e j × e k . (30)
From this last equation, we can easily see that the contravariant vectors e i = ∇ξ i are related with the covariant vectors

e i = ∂x ∂ξ i through Je i = e j × e k , i, j, k are cyclic. (31)

Conservation laws under mappings

The divergence operator in the conservation law (24) has been applied in the physical space. This operator under the coordinate transformation x = X(ξ) is given by

∇ x • f = 1 J 3 i=1 ∂ ∂ξ i Je i • f := 1 J ∇ ξ • f, (32)
with f = f 1 , f 2 , f 3 the tensor of contravariant fluxes and it is defined through its components by f i ≡ Je i • f . The conservation law in the reference element takes the form

∂u ∂t + 1 J ∇ ξ • f = 0. ( 33
)
Observe that u = u(ξ, t), f = f u(ξ, t) , and J = J(ξ). The variational formulation of the transformed conservation law is, after multiplying with a test function φ = φ(ξ),

Ω ∂u ∂t + 1 J ∇ ξ • f φ dx = 0. ( 34
)
Transforming the integral to the reference element E = [-1, 1] 3 , and reminding that dx = J(ξ) dξ, we obtain the variational formulation of the transformed equation on the reference element

E J(ξ) ∂u ∂t + ∇ ξ • f φ(ξ) dξ = 0. (35)
Finally, the weak formulation of the equation ( 33) is obtained by integrating by parts the flux term in the equation ( 35)

E J(ξ) ∂u ∂t φ(ξ) dξ - E f • ∇ ξ φ(ξ) dξ + ∂E f • ñ * φ(ξ) dσ = 0, ( 36
)
where ñ is the unit normal vector to the surface at the faces of the reference element. Observe that for a given direction in the reference space (e.g. ñ can be the normal unit vector at the face of the reference element in the direction +ξ i ) we have that

f • ñ = f • n = (f • n) n , ( 37
)
where n is the normal vector at the surface in the physical space, which is computed from the metric terms evaluated at the surfaces, and n is the unit normal vector n = n/ n ; the components of the vector n are calculated from

n l = 3 i=1 Je i l ñi , l = 1, 2, 3. ( 38
)
Its norm represents the area of the face element in the physical space and it is given by

σ = n = 3 l=1 Je i l 2 . (39)
The superindex i means that n is associated with ñ in the direction of the coordinate ξ i as we pointed out before. We will not use indices for this denomination and henceforth it will be understood from the context which direction it corresponds. Therefore

f • ñ * = f • n * = f • n * σ. ( 40
)
At cell interfaces the solution is double-valued (one value is obtained from the data inside the current cell and the other value corresponds to the value calculated by the neighboring cell), just like the in the finite volume method. This discontinuity at cell interfaces implies that the flux has to be calculated by introducing the so-called numerical flux, and we write that in the variational formulation as f • n * . In the finite volume community these numerical fluxes are computed usually with a Riemann solver. There is a large variety of Riemann solvers available for different systems of conservation laws, see for example the book by Toro [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] for an exhaustive treatment of this topic. We have to mention that the numerical flux in the normal direction is defined as a function of the left and right states at the element interface:

f * = f * (u L , u R ; n) = f • n * ,
and it is computed by solving a one-dimensional Riemann problem for the states in the Gauss-Legendre quadrature points on the surface of the element. The Riemann solver used in our computations is the very well known Rusanov numerical flux [START_REF] Rusanov | [END_REF]. The Rusanov numerical flux is given by

f * (u L , u R ) = 1 2 f (u L ) + f (u L ) -|λ max |(u L -u R ) , (41)
where λ max is the largest local wave speed, which guarantees the stability of the scheme. Finally, the weak form of the equation ( 33) reads

E J(ξ) ∂u ∂t φ(ξ) dξ Time Derivative Integral - E f • ∇ ξ φ(ξ) dξ Volume Integral + ∂E f • n * σφ(ξ) dσ Surface Integral = 0. ( 42
)
x y

-1 +1 -1 +1 ξ η X -1 (x)

Discretization on hexahedral elements

The discontinuous Galerkin spectral element method is a discretization method specially designed for the approximate solution of conservation laws on quadrilateral/hexahedral computational domains. The partial differential equation is expressed in weak form and the solution and fluxes are approximated by using high-order Lagrangian interpolants. We stress that the polynomials are represented in a nodal form. The integrals are approximated by quadrature (in our case, Gauss quadrature), and the nodal points, where the solution is computed, are the Gauss-Legendre quadrature points. Along the boundaries, the nodes are chosen to be also Gauss-Legendre quadrature points in order to compute the boundary integrals. In the Figure 2 these nodal points are depicted for the two-dimensional case. The ansatz discussed before will be used in the following for the purpose of discretizing each integral in the equation [START_REF] Powell | An approximate Riemann solver for magnetohydrodynamics[END_REF].

In every hexahedral element, we approximate the vector of conserved variables and the contravariant fluxes by polynomials, which basis functions are obtained from the tensor-product of one-dimensional Lagrange polynomials of degree N. The state is given by

u = N i, j,k=0 ûi jk ψ i jk (ξ), (43)
with ûi jk the nodal degrees of freedom. The components of the contravariant fluxes are written as

f m u(ξ) = N i, j,k=0 f m i jk ψ i jk (ξ), m = 1, 2, 3. (44)
The coefficients f m i jk are computed in terms of the covariant fluxes, as it was pointed out before

f m i jk = f m i jk u(ξ i jk ) = 3 l=1 Je m l (ξ i jk )f l u(ξ i jk ) . (45)
The basis functions ψ i jk (ξ) used before have a tensor-product form

ψ i jk (ξ) = i (ξ 1 ) j (ξ 2 ) k (ξ 3 ), (46)
with i the Lagrange polynomials

j ξ = N i=0 i j ξ -ξ i ξ j -ξ i , (47)
with the property

j ξ i = δ i j , i, j = 0, . . . , N. (48)
Time derivative integral. In the time derivative integral in the equation ( 42), we insert the approximation for u given in [START_REF] Dedner | [END_REF]. Additionally, we assume that the test function is the same as the basis functions, φ = ψ i jk (this is the Galerkin ansatz), and finally we approximate the integrals by Gaussian quadrature

∂ ∂t E J(ξ) u(ξ) φ(ξ) dξ = ∂ ∂t 1 -1 1 -1 1 -1 J(ξ) N l,m,n=0 ûlmn ψ lmn (ξ) ψ i jk (ξ) dξ = ω i ω j ω k J(ξ i jk ) d ûi jk dt ,
where we have made use of the property (48) of the Lagrange interpolating polynomials. The ω i are the weights in the Gaussian quadrature.

Volume integral. The volume integral in the equation ( 42) involves the tensor f. This integral can be written as

E f • ∇ ξ φ(ξ) dξ = 3 d=1 E f d u(ξ) ∂φ(ξ) ∂ξ d dξ = ω j ω k N λ=0 D λi f 1 λ jk ω λ + ω i ω k N µ=0 D µ j f 2 iµk ω µ + ω i ω j N ν=0 D νk f 3 i jν ω ν ,
where we have introduced the differentiation matrix D λi , which is defined as

D i j = d j (ξ) dξ ξ=ξ i , i, j = 0, . . . , N. (49)
Surface integral. The last integral we will compute is the surface integral in the equation [START_REF] Powell | An approximate Riemann solver for magnetohydrodynamics[END_REF]. Observe that this surface integral is the sum of the surface integrals over each face of the reference element 1) . The vector ξ l means that the points are on the surface ∂E l , and σ l is the area of the element face l. We will compute only the integrals in the direction ±ξ 1 . Remember that the sign in the directions is included in the normal vector. The numerical flux f * = f • n * on the faces ±ξ 1 of the reference element E, as we pointed out before, are calculated by using the states u L and u R . The state u L is obtained from the extrapolation of the inner state to the element boundaries and the state u R from the extrapolation in the neighboring element (or in the case that the face to be evaluated is located on the physical boundary, from the boundary conditions imposed in the problem). This extrapolation is carried out on the Gauss-Legendre quadrature points on the face of the reference element. Therefore, the numerical fluxes on the faces ±ξ 1 (i.e. the faces at ξ 1 = ±1) are approximated by

∂E f • n * σφ(ξ) dσ = 6 l=1 ∂E l f • nl * σ l φ(ξ l ) dσ l = ω j ω k f * σ +ξ 1 jk i (+1) + f * σ -ξ 1 jk i (-1) + ω i ω k f * σ +ξ 2 ik j (+1) + f * σ -ξ 2 ik j (-1) + ω i ω j f * σ +ξ 3 i j k (+1) + f * σ -ξ 3 i j k (-
f • nl * σ l l=±ξ 1 = f * (±1, ξ 2 , ξ 3 ) σ ±ξ 1 = N j,k=0 f * σ ±ξ 1 jk j (ξ 2 ) k (ξ 3 ). ( 50
)
Semi-discrete formulation. Once discretized all integrals appearing in the equation [START_REF] Powell | An approximate Riemann solver for magnetohydrodynamics[END_REF], we are able to write down the semi-discrete formulation of the DGSEM on hexahedral elements. Before we write the final equation, we define some useful quantities, which can be precomputed and stored in memory. The operators to be defined are

Di j = -D ji ω j ω i , ˆ i ±1 = i ±1 ω i , i, j = 0, . . . , N. (51)
The semi-discrete formulation has the following form

d ûi jk dt = - 1 J i jk N λ=0 Diλ f 1 λ jk + f * σ +ξ 1 jk ˆ i (+1) + f * σ -ξ 1 jk ˆ i (-1) - 1 J i jk N µ=0 D jµ f 2 iµk + f * σ +ξ 2 ik ˆ j (+1) + f * σ -ξ 2 ik ˆ j (-1) - 1 J i jk N ν=0 Dkν f 3 i jν + f * σ +ξ 3 i j ˆ k (+1) + f * σ -ξ 3 i j ˆ k (-1) . (52)
The semi-discrete formulation given by equation ( 52) consists of one-dimensional DGSEM operators (the quantities in square brackets), each of them applied along a direction.

Discretization of the GLM divergence cleaning

The equations for B x and ψ are decoupled from the rest of the MHD-GLM system. In fact, cleaning a vector quantity B with divergence errors results in the linear system

∂ ∂t B x ψ + 0 1 c 2 h 0 ∂ ∂x B x ψ =        0 - c 2 h c 2 p ψ        . (53)
The numerical flux is then derived as the solution of the local Riemann problem with left-hand state (B x,l , ψ l ) T and right hand state (B x,r , ψ r ) T as

B x,m = 1 2 (B x,r + B x,l ) - 1 2c h (ψ r -ψ l ), (54)
ψ m = 1 2 (ψ r + ψ l ) - c h 2 (B x,r -B x,l ). ( 55
)
It is possible to employ for the quantities B x and ψ the Riemann solver used in the full MHD system, but in this work we use the exact solution (54) of the linear Riemann problem (53). Dedner et al. [START_REF] Dedner | [END_REF] suggest to use the solution (54) as input for the Riemann solver used in the solution of the other conserved quantities.

There are basically two possible ways to deal with the source term in equation (6b). The first one consists in adding this source term to the hyperbolic update, that is to the semi-discrete scheme (52) (see [58]). The second one is based on an operator-splitting approach. Following the idea presented in [START_REF] Dedner | [END_REF], we first solve the homogeneous MHD-GLM system in a so-called hyperbolic step, and then we consider the source term in the source step. The scalar field ψ is then

ψ (∆t) = ψ (0) exp -α p c h ∆h/∆t , with α p = ∆h c h c 2 p , (56)
where ψ (0) has been computed in the hyperbolic step, and ∆h = min (∆x, ∆y, ∆z) is the minimum mesh size. This approach is very simple to implement and is unconditionally stable [START_REF] Dedner | [END_REF].

Time discretization

The semi-discrete scheme ( 52) is solved explicitly using the so-called method of lines. The method of lines acquires the accuracy order of the integrator used to solve the system of ordinary differential equations, under the condition that the spatial discretization is of the same order of accuracy or higher [START_REF] Schiesser | The Numerical Method of Lines: Integration of Partial Differential Equations[END_REF]. The system of ordinary differential equations ( 52) is solved with the family of strong-stability preserving Runge-Kutta Methods (SSPRK) [START_REF] Shu | [END_REF][61][62]. The main idea of these methods resides in assuming that the first-order forward Euler method is strongly stable under the total variation norm (and in general, any given norm) and a suitable time step restriction. From this, the aim is to construct a higher-order time discretization that preserves strong stability under the TVD norm, and maybe, with a different time step restriction [START_REF] Shu | [END_REF]. The SSPRK methods are desirable in problems with discontinuities and strong shocks, because they guarantee that, as part of the time-integration process, no additional oscillations are introduced.

The Courant-Friedrichs-Lewy (CFL) condition is a necessary condition for stability of any explicit one-level numerical scheme [63]. It establishes that the domain of dependence of the solution is contained in the numerical domain of dependence of the numerical method. This condition is applied to constraint the time step. For the discontinuous Galerkin method, Cockburn and Shu [1] gave an approximate CFL condition for linear stability, and it has the following form

∆t = C CFL 2N + 1 min k ∆x |λ k,x | , ∆y |λ k,y | , ∆z |λ k,z | . ( 57
)

Shock capturing

In the discontinuous Galerkin method is common that some unphysical states are produced in the intermediate stages of the simulation (e.g., during the interpolation step, or the conservative to primitive variables conversion), especially when the flow involves shocks waves or strong rarefactions, and very high-order are used (interpolation order larger than second order). These unphysical states come from the loss of positivity of the density or the pressure, or the generation of super-luminal speeds. Several shock capturing schemes have been considered in the past, like the moment limiter [16,17], the artificial viscosity [START_REF] Persson | Proceedings of the AIAA 44th Aerospace Sciences Meeting and Exhibit, AIAA-2006-0112[END_REF][START_REF] Casoni | [END_REF], and the WENO limiter [20][21][22]. It is very well known that in the vicinity of a discontinuity, shock-capturing schemes will reduce the solution accuracy to first-order [64]. This is because the location of the shock is proportional to the mesh size [65]. Taking this into account, it is natural to construct a robust shock-capturing by considering only two stages: the shock detection and a further order reduction of the scheme in the vicinity of the shock. This kind of schemes are called the fallback approaches [32,38,66]. For robustness and efficiency reasons, the hybrid DG/FV approach [START_REF] Sonntag | [END_REF]27], which is based on the fallback approach (see table 1), will be used as the main building block for the nonlinear stabilization of the solution when shocks take place. The hybrid DG/FV is constructed in such a way that, in regions of smooth flows, the DGSEM method is employed, and those parts of the flow having shocks, the DGSEM elements are interpreted as quadrilateral/hexahedral subdomains. In each of these subdomains, the nodal DG solution values are used to build a new local domain composed now of finite volume subcells, which are evolved with a robust finite volume method with third order WENO reconstruction. In the following we will discuss with more detail the hybrid DG/FV algorithm.

Hybrid DG/FV shock capturing

The hybrid DG/FV shock capturing method used in this work consists basically in two steps: the first step is detecting the elements in the DG computational domain that contain oscillating polynomials, which cause instabilities in the computation and make the code to blow up. The second step is replacing the troubled elements with subdomains made of (N + 1) d finite volume subcells each, with d the space dimension. Then the subdomains are evolved with a robust finite volume method which uses a second or third order reconstruction operator, typically MUSCL (second order) or WENO3 (third order).

Let us represent the solution in the DG reference element E at time level t n by u n = N i, j,k=0 ûn i jk ψ i jk (ξ). In the context of DGSEM schemes, the reference element E has (N + 1) 3 degrees of freedom. The DG solution is computed in the (N + 1) 3 quadrature points. The reference element E is then decomposed into (N + 1) 3 equidistant FV subcells E i jk , in such a way that E = i jk E i jk , see figure 3. While in the DG reference element E the solution is represented by Table 1: Shock capturing strategy for the higher-order discontinuous Galerkin methods discussed in this work [27,32,38,39,66].

Checking of physical quantities with restrictions

Positivity of density Positivity of pressure Maximum propagation speed (only in the relativistic case)

Detection of regions with strong shocks

Shocks indicators Oscillations indicators Marking of troubled points/cells and direct neighbors

Special treatment in troubled regions

Employment of robust Riemann solvers Employment of robust first-order Godunov scheme Employment of robust second/third order numerical scheme Figure 3: Hybrid DG/FV. DG element with (N + 1) 2 inner nodal points (the blue filled circles) where the solution is computed, and which is approximated by a tensor-product of Lagrange polynomials of degree N (left). Corresponding FV subdomain made of (N + 1) 2 equidistant subcells (right), with the barycenter of every subcell shown (the blue filled circles), and the points where the numerical fluxes are computed with the finite volume scheme (the white filled squares), and the borders of the subcells (dashed lines). a polynomial, in the FV subdomain the solution at time level t n is represented by (N + 1) 3 piecewise constant subcell averages v n i jk . These values are obtained from the DG polynomial itself by using L 2 projection, that is

v n i jk = P FV (u n ) = 1 E i jk E i jk u(x, t n ) dx. ( 58
)
The inverse operator P DG = P -1 FV has to satisfy the following constraint: The integral conservation of the cell average over the DG element

E u(x, t n ) dx = E v(x, t n ) dx. ( 59
)
We solve then the conservation law for the subcell cell averages

∂v ∂t + 1 J ∇ ξ • f = 0, (60)
by using a second/third order FV method, see [START_REF] Sonntag | [END_REF]33]. Summarizing, given the DG state u n i at time level n, the oscillation indicator operator OI is applied to u n i , ∀i, and only the troubled elements are flagged for shock capturing. For those not flagged elements, the state u n is evolved with the DG solver operator DG, but if the element is marked as a troubled cell, then we use the operator P FV to project the state u n onto a FV subcells subdomain of piecewise cell averages, yielding v n ; then we evolve the states v n in the FV subdomain with the FV solver operator FV, to get v n+1 . The solution v n+1 is then projected back onto the DG element with the projection operator P DG = P -1 FV . Since the stencil of the reconstruction operator used in the FV subdomain needs at least two neighboring cell data (one left and one right to the subcell of interest), those subcells at edges of the FV subdomain require boundary values, which are obtained from the neighboring elements (see figure 4). A flow diagram for the DGSEM is depicted in the figure 5.

Finite volume schemes

For the sake of completeness, we briefly describe the foundations of the finite volume method with third-order WENO reconstruction. Let us consider the following system of conservation laws on structured Cartesian meshes

∂u ∂t + ∂f (u) ∂x + ∂g(u) ∂y + ∂h(u) ∂z = 0. (61)
After the integration of the equation ( 61) over the computational cell

Ω i jk = x i-1 2 , x i+ 1 2 × y j-1 2 , y j+ 1 2 × z k-1 2 , z k+ 1 2 , we get the semi-discrete scheme du i jk dt = - fi+ 1 2 , jk -fi-1 2 , jk ∆x - ĝi, j+ 1 2 ,k -ĝi, j-1 2 ,k ∆y - ĥi j,k+ 1 2 -ĥi j,k-1 2 ∆z . ( 62
)
The quantity u i jk is the spatial average of u in the cell Ω i jk at time t

u i jk = 1 Ω i jk Ω i jk u(x, y, z) dz dy dx, (63)
Reading with Ω i jk = ∆x∆y∆z, and fi± 1 2 , jk , ĝi, j± 1 2 ,k , and ĥi j,k± 1 2 are spatial averages of the physical fluxes over the cell faces x i± 1 2 , y j± 1 2 , and z k± 1 2 , respectively, at time t

fi± 1 2 , jk = 1 σ jk σ jk f u(x i± 1 2 , y, z) dz dy, ĝi, j± 1 2 ,k = 1 σ ik σ ik g u(x, y j± 1 2 , z) dz dx, ĥi j,k± 1 2 = 1 σ i j σ i j h u(x, y, z k± 1 2 ) dy dx, (64)
with the surfaces elements defined by

σ i j = x i-1 2 , x i+ 1 2 × y j-1 2 , y j+ 1 2 , σ jk = y j-1 2 , y j+ 1 2 × z k-1 2 , z k+ 1 2 , and σ ik = x i-1 2 , x i+ 1 2 × z k-1 2 , z k+ 1 2 .
The area of the faces are then σ i j = ∆x∆y, σ jk = ∆y∆z, and σ ik = ∆x∆z. The semidiscrete scheme ( 62) is an exact relation as well as the averaged quantities (63) and (64). For the purpose of this work, we approximate these integral as follows fi± 1 2 , jk =

f u(x i± 1 2 , y j , z k ) ∆y∆z , ĝi, j± 1 2 ,k = g u(x i , y j± 1 2 , z k ) ∆x∆z , ĥi j,k± 1 2 = h u(x i , y j , z k± 1 2 ) ∆x∆y . ( 65
)
The point-wise values u(x i± 1 2 , y j , z k ), u(x i , y j± 1 2 , z k ), u(x i , y j , z k± 1 2 ) are reconstructed from the cell averages u i jk with a third-order WENO scheme. Observe that two sets of reconstructed values at a given face are present: those obtained through the use of the cell Ω i jk as the main cell in the reconstruction procedure and those associated with the neighboring cell. These are known as the left and right values at the face interface: u L and u R . Finally the fluxes are evaluated by replacing f , g, and h by a monotone flux. We employ the same Riemann solver used in the DGSEM scheme, that is the Rusanov numerical flux [START_REF] Rusanov | [END_REF].

WENO reconstruction

The basic idea of the WENO schemes is based on an adaptive reconstruction procedure to obtain a higher-order approximation on smooth regions while the scheme remains non-oscillatory near discontinuities [15]. In each stencil, a polynomial is reconstructed from the cell averages ūi of the solution and later a weighted combination of all these polynomials is constructed. The nonlinear weights are computed taking into account the smoothness of every polynomial in its respective stencil.

Let us approximate the value of the function u(x) at the points x i-1 2 , and x i+ 1 2 , by using polynomials of degree N = 1. For the stencil S k we get the corresponding values u (k) i- 1 2 , and u (k)

i+ 1 2
. Let us write down the stencils

S 1 = {I i-1 , I i } , S 2 = {I i , I i+1 } , (66)
where

I i = x i-1 2 , x i+ 1 2 .
The corresponding interpolated values for the point x i-1 2 are

u (1) i-1 2 = 1 2 ūi-1 + ūi , u (2)
i-

1 2 = 1 2 3ū i -ūi+1 , (67)
and the corresponding interpolated values for the point x i+ 1 2 are u (1)

i+ 1 2 = 1 2 -ū i-1 + 3ū i , u (2)
i+ 1 2 = 1 2 ūi + ūi+1 . (68)
If we choose the large stencil

S = {I i-1 , I i , I i+1 } , (69)
which is the union of all 2 stencils S k , then we are able to find an interpolation polynomial p(x) of degree at most 2, satisfying p(x j ) = u j for i -N ≤ j ≤ i + N and giving the approximations u i-1 2 ≡ p(x i-1 2 ), u i+ 1 2 ≡ p(x i+ 1 2 ),

u i-1 2 = 1 6 2ū i-1 + 5ū i -ūi+1 , u i+ 1 2 = 1 6 -ū i-1 + 5ū i + 2ū i+1 , (70)
provided that the function is smooth in the large stencil S . The WENO idea consists in choosing the final approximation as a convex combination of the approximations u (k)

i+ 1 2 u i+ 1 2 = N+1 k=1 ω k u (k) i+ 1 2 , (71)
where ω k ≥ 0, N+1 k=1 ω k = 1. The approximations calculated by using the large stencils can be written as a linear convex combination of the approximations u (k)

i+ 1 2 based on the small stencils S i u i+ 1 2 = N+1 k=1 γ k u (k) i+ 1 2 , ( 72
)
where the constants γ j satisfy N+1 j=1 γ j = 1, are usually referred to as the linear weights in the WENO literature. For the third-order WENO scheme, the linear weights for the points x i- 1 2 , and x i+ 1 2 , are

x i-1 2 : γ 1 = 2 3 , γ 2 = 1 2 , x i+ 1 2 : γ 1 = 1 2 , γ 2 = 2 3 . ( 73
)
The choice of the nonlinear weights ω k relies on the smoothness indicator β k , which measures the relative smoothness of the function u(x) in the stencil S k . The larger this smoothness indicator β k , the less smooth the function u(x) is in the stencil S k . The smoothness indicators for the third-order WENO scheme are

β 1 = (ū i-1 -ūi ) 2 , β 2 = (ū i -ūi+1 ) 2 . ( 74
)
Notice that these smoothness indicators are quadratic functions of the values of u(x) in the relevant stencils. Equipped with these smoothness indicators, we can now define the nonlinear weights as

ω k = ωk N i=1 ωi , with ωk = γ k (ε + β k ) 2 . ( 75
)
Here ε is a small positive number used to avoid the denominator becoming zero and is typically chosen to be ε = 10 -6 , but in our calculations we have used ε = 10 -24 .

Oscillations indicators

In our a priori shock capturing strategy for discontinuous Galerkin methods, the first step consists in determining which elements contain oscillating polynomials, that is, those elements where shock waves or discontinuities are producing these kind of behavior in the approximating polynomials. The oscillation indicator is applied to every element at each stage of the explicit Runge-Kutta time discretization. Once the troubled element is flagged, the next step is to apply the hybrid DG/FV approach, where the DG element is used to build a new block made of finite volume subcells, and which is evolved using a robust second/third order finite volume scheme. In this work we mainly use the Persson indicator. [START_REF] Persson | Proceedings of the AIAA 44th Aerospace Sciences Meeting and Exhibit, AIAA-2006-0112[END_REF], and it measures cell-wise the maximum energy decay in the highest and second highest degrees of freedom of the DG polynomial for the state variable u:

Persson indicator. The Persson indicator was developed by Persson and Peraire

η(u) = log 10                        max                                               N DOF (N) j=N DOF (N-1)+1 û j 2 N DOF (N) j=1 û j 2                        ,                        N DOF (N-1) j=N DOF (N-2)+1 û j 2 N DOF (N-1) j=1 û j 2                                                                      , ( 76
)
where N DOF (N) denotes the number of degrees of freedom for a DG polynomial of degree N and û j the j th degree of freedom of the polynomial representation of u. Observe that this indicator takes into account the influence of odd/even order effects in the DG polynomial. The typical variable used by the Persson indicator in order to determine if a cell suffers of spurious oscillations is the pressure, since shock waves imply a jump in the pressure and not in the density. With this is distinguished shock waves and contact discontinuities.

Numerical computations

In this section we will present several numerical calculations of the MHD and SRHD equations with the hybrid DG/FV scheme. Convergence tests, and one-and two-dimensional flow problems involving shock waves are discussed. We emphasize on the employment of very high order polynomials for the discretization of the vector of conservative quantities as well as the fluxes. In the two-dimensional computations the solution has been plotted in a hybrid manner, that is we plot the DG polynomial solution along the FV subdomains with their corresponding subcell distribution.

Table 2: Convergence rates for the MHD equations with initial condition given by the smooth circularly polarized Alfvén wave. Results for L 2 norm error of the x-component of magnetic field are given. The computations were performed with the pure DGSEM, and the hybrid DG/FV with all elements flagged as troubled. The schemes use polynomial degrees N = 3, 5, 7, that is for orders of accuracy O(4), O(6) and O(8).

DGSEM

Hybrid

Propagation of a smooth circularly polarized Alfvén wave

In order to check the order of accuracy of the algorithm with high order polynomials, we need to know a smooth exact solution of the MHD equations. The propagation of a smooth circularly polarized Alfvén wave is an example of such exact solution. Therefore it is suitable for performing a convergence test. Following Mignone et al. [45], the simulation parameters are set as follows: the two-dimensional computational domain has dimensions L x = 1.0 and L y = 0.5, and n x = 2n y , where n x is the number of elements in x-direction, and n y is the number of elements in y-direction. The Alfvén wave travels along the diagonal of the domain, and forms an angle α = arctan(0.5) ≈ 26.6 • with respect to the x-axis. The primitive variables are set in the following way: the density is set to ρ = 1, and the pressure p = 0.1. Let us define the quantities v 0 = 0.1, B 0 = 0.1, and x = (x cos α + y sin α). The components of the velocity and magnetic field perpendicular to the wave vector are given by v ⊥ = v 0 sin(2πx ) and B ⊥ = B 0 sin(2πx ), and the velocity and magnetic field in z-direction are v z = v 0 cos(2πx ) and B z = B 0 cos(2πx ). The perpendicular components v ⊥ and B ⊥ can be written as linear combination of B x and B y through B ⊥ = B y cos α -B x sin α, and B = B x cos α + B y sin α. An ideal gas EOS with adiabatic exponent γ = 5/3 is used. Periodic boundary conditions are prescribed everywhere. The errors and the rate of convergence are reported in the table 2, for computations at simulation time t = 1, and for polynomials of degree N = 3, 5, 7. From this table we can observe that the experimental order of convergence (EOC) is achieved by the pure DGSEM, but the hybrid DG/FV scheme with all elements flagged as troubled is third order accurate, the order of the WENO3 reconstruction.

One-dimensional Riemann problems

Here, we run a series of one-dimensional Riemann problems, which have been extensively used to check the robustness of numerical schemes for conservation laws. We analyze the test problems reported in [28,31,67]. These problems are characterized by the appearance of different shock patterns after the breaking of the membrane. In order to clarify, a Riemann problem is an initial-value problem for a conservation law defined as The membrane is localized in the point x m = 0.5. An ideal gas EOS with adiabatic index γ = 5/3 is assumed for all Riemann problems. Dirichlet boundary conditions are imposed at domain boundaries. The computational domain is decomposed into 500 elements, for all polynomials degrees. The hybrid DG/FV scheme used for solving these problems employs polynomials of degree N = 2, 4, 6. The four initial conditions are given in table 3, and the exact solutions have been obtained with the exact Riemann solver developed by Torrilhon [68], and they are plotted as red lines (see figure 6). The description of the flow patterns is presented below together with some comments regarding the capabilities of the hybrid DG/FV scheme with WENO3 reconstruction in the FV solver.

∂u ∂t + ∂f (u) ∂x = 0, u(x, 0) =        u L (x) for x < x m , u R (x) for x > x m , (77)
Riemann problem 1 (RP-1). The flow pattern features, after the breaking of the membrane, two fast shocks, one weak slow shock, one slow rarefaction and one contact discontinuity. In the figure 6 (top) are depicted plots of the density, thermal pressure, and the y-component of the magnetic field at simulation time t = 0.08. The discontinuities are very well captured. Very small oscillations are observed for all polynomials degrees.

Riemann problem 2 (RP-2). This problem comprises multiple weak discontinuities, namely, two fast shocks, two slow shocks, two rotational discontinuities and one contact discontinuity. The shocks travel from each side of the contact discontinuity. In the figure 6 (second row, from top) are depicted plots of the density, thermal pressure, and the y-component of the magnetic field at simulation time t = 0.2. Thanks to the subcell resolution of the hybrid DG/FV scheme with WENO3 reconstruction, we observe a very good agreement with the exact solution.

Riemann problem 3 (RP-3). In this problem are present a switch-on slow rarefaction wave. The following flow structures are developed: a hydrodynamic rarefaction, a switch-on slow rarefaction, a contact discontinuity, a slow shock, a rotational discontinuity, and a fast rarefaction. In the figure 6 (second row, from bottom) are depicted plots of the density, thermal pressure, and the y-component of the magnetic field at time t = 0.16. All structures mentioned before are present in the figure, from left to right. The discontinuities have been very well captured and they are oscillations-free.

Riemann problem 4 (RP-4). This problem is characterized by the formation of the so-called compound structures. These structures involve a shock and a rarefaction wave traveling together. Besides, the test involves a left-going slow compound wave. In the figure 6 and it has been discretized by using a mesh of 120 × 120 elements.

Current sheet

In this test problem, a region is uniformly filled with a gas at rest. The initial setup of the magnetic field is such that this switches signs at the slices x = +0.25 and x = -0.25. This configuration is perturbed with a sinusoidal velocity function in y, which generates nonlinear, linearly polarized Alfvén waves. These Alfvén waves turn into magnetosonic waves due to the magnetic pressure does not remain constant. The two current sheets at x = ±0.25 originate magnetic reconnection. The magnetic reconnection drives highly over-pressurized regions, which launch magnetosonic waves transverse to the field, causing magnetic energy to be transformed into thermal energy [70]. Large magnetic field gradients are produced in the proximity of the points where the magnetic reconnection takes place. This problem is useful to check if the numerical scheme is capable of handling magnetic reconnection and shock waves.

The computational domain is Ω = [-0.5, +0.5] × [-0.5, +0.5]. The density and the pressure are uniform in the whole computational domain, with ρ = 1.0, and p = 0.5β. The parameter β = 0.1 is the ratio of gas pressure to magnetic energy density. The components of the velocity are v y = v z = 0, and v x = A sin(2πy), where A = 0.1 is a parameter that is typically used to test the robustness of the algorithm. The magnetic field components are B x = B z = 0, and B y = 1 for |x| > 0.25 and B y = -1 otherwise. An ideal gas EOS with adiabatic index γ = 5/3 is assumed. Periodic boundary conditions are set at all four boundaries of the domain. The problem runs up to final time t = 10. In the figure 8 are depicted the density at times t = 5.0 and t = 7.5. At time t = 5 the solution is plotted along with the mesh. The FV subcells are located in those regions with strong shocks. Magnetic islands emerge and merge each other. The numerical results are similar to those computed with a high-order finite volume scheme [33]. Due to the complex structures present in this problem along with the large gradients in the flow variables, it is usual that a high-order method without an appropriate shock-capturing strategy crashes. We stress the robustness of the hybrid DG/FV scheme for this kind of problems, for long time runs.

Orszag-Tang vortex

The Orszag-Tang vortex was first analyzed by Orszag and Tang [71] for the incompressible MHD equations. Nowadays this test is widely used to assess the capacity of the numerical scheme for handling the formation and the interactions of MHD shocks in the context of the compressible MHD equations [69,[72][73][74]. Initially, the flow profile consists in the superposition of a velocity vortex with a magnetic vortex. This configuration is rather unstable, which originates a broad range of MHD waves, interacting with each other, making a transition towards turbulence.

The density and the pressure are uniform in the whole domain, with ρ = γ 2 and p = γ. With this choice of the Density, t = 0.5 density and pressure, the sound speed is c s = 1. The initial velocity profile is given by

v x = -sin(2πy), v y = + sin(2πx), v z = 0, (80)
and the magnetic field is

B x = -sin(2πy), B y = + sin(4πx), B z = 0. (81)
An ideal gas EOS with adiabatic index γ = 5/3 is assumed. The computational domain is Ω = [0, 1] × [0, 1], and it has been decomposed into 120 × 120 elements. Periodic boundary conditions are set at the boundaries of the domain. The simulation time is t = 1.0, and in the figure 9 are depicted the density and gas pressure at time t = 0.5. For the plot of the density, the mesh is also presented. The FV subdomains are evolved with a third-order finite volume WENO scheme. The hybrid DG/FV method is quite stable and robust even for this kind of supersonic flow configurations.

Observe that the contours are not smooth. This is because we have chosen a small threshold for deciding where to use a FV solver and where a DG solver. We wanted to use as much DG cells as possible for this problem. The polynomial degree we employed is N = 5. The divergence cleaning of Dedner et al. [START_REF] Dedner | [END_REF] was used, with c r = 0.18, and c h determined by the maximum propagation speed in the system. Additionally, in the figure 10 is plotted the pressure along the slices y = 0.4277 and y = 0.3125 at time t = 0.5. Computations were made with a very high-order finite difference WENO5 scheme (mesh of 600 2 cells) and the hybrid DG/FV method with N = 5 (mesh of 120 2 elements).

The profiles match very well, except at boundaries, where a slight shift is observed. It is also noticeable the lack of oscillations in these slices for the hybrid DG/FV scheme, considering the finite difference as reference solution.

Cylindrical blast wave

This problem features a cylindrical region located in the center of a domain. The cylindrical region is filled with a magnetized overpressured gas. After the system is released, strong shock waves move outwards [72,75,76]. The computational domain is Ω = [0, 1] × [0, 1]. The density and pressure are uniform in the whole domain, with ρ = 1 and p = 0.1. In the cylindrical region (xx c ) 2 + (yy c ) 2 < R, with (x c , y c ) = (0.5, 0.5), and R = 0.1, the pressure is p = 10.0. The initial velocity is zero, v x = v y = v z = 0, and the magnetic field is B x = B y = B z = 1/ √

2. An ideal gas EOS with adiabatic index γ = 5/3 is employed. Periodic boundaries are also assumed. The simulation time is t = 0.5, and in the figure 11 are depicted the density and pressure at time t = 0.1. The computational domain is decomposed into 120 × 120 elements. The spatial discretization is the hybrid DG/FV scheme with WENO3 reconstruction operator. The polynomial used for the simulation has degree N = 5. The divergence cleaning of Dedner et al. [START_REF] Dedner | [END_REF] was employed, with c r = 0.18. The FV solver is able to capture the shocks arising after the release of the cylindrical region. The solution of this very difficult problem with the hybrid scheme confirms the robustness of the shock-capturing approach. Very small oscillations can be observed in the contour lines of the pressure plot.

Rotor problem

This problem consists in a high-density, rapidly spinning fluid in a low-density fluid. Initially, the rotor and the ambient medium are subject to an uniform magnetic field. The rotating fluid launches torsional Alfvén waves into the fluid at rest [77]. The problem setup is as follows: in the ambient medium, the initial density and pressure are uniform, with ρ = 1.0 and p = 1.0, and the velocity is zero, v x = v y = v z = 0. In the interior of the cylindrical rotor (0.0 ≤ (xx c ) 2 + (yy c ) 2 ≤ 0.1, with (x c , y c ) = (0.5, 0.5)), the fluid density is ρ = 10.0, and the pressure is the same as in the ambient fluid. The fluid inside rotor has an angular velocity ω such that v = ωr = 1 at r = 0.1. We applied a linear taper to the velocity and the density. This taper is applied only in the region 0 ≤ r ≤ 0.105 in such a way that the density and the velocity have the same values that those of the ambient fluid. Finally, the components of the magnetic field are B x = 5.0, B y = 0, and B z = 0 in the whole computational domain.

The computational domain is Ω = [0, 1] × [0, 1]. An ideal gas EOS with γ = 7/5 is employed. Dirichlet boundary conditions are assumed. The computational domain is discretized using 120 × 120 elements. The polynomial used in the calculations has degree N = 5. For the treatment of the solenoidal constraint, the divergence cleaning of Dedner et al. [START_REF] Dedner | [END_REF] has been employed. The simulation time is t = 0.5, and in the figure 12 are depicted the density (along the mesh), pressure, Mach number, and the magnitude of the magnetic field at time t = 0.25. In the figure 13 is plotted the pressure along the slices x = 0.5 and y = 0.5 at time t = 0.25. The computations were done with a very high-order finite difference WENO5 scheme and the hybrid DG/FV method with N = 5. Observe the good agreement with the solution obtained with the high-order finite difference scheme. A lack of oscillations is also noticeable. We point out the difference at boundaries in this plot. This is an effect of the boundary conditions implemented in the hybrid DG/FV scheme. A similar behavior also occurs for other computations with open boundaries in this work.

Relativistic hydrodynamics 4.5.1. Convergence test

We start by performing a convergence test of the hybrid DG/FV scheme for the SRHD equations. We consider a exact solution of the SRHD equations proposed by He and Tang [78]. The problem consists in a wave propagating in the physical domain [0, 2/ √ 2] × [0, 2] at an angle θ = 30 • relative to the horizontal axis. The initial profile is with A = 0.2, and v 0 = 0.2. This problem has the exact solution

ρ = 1 + A sin (φ(0)) , v x = v 0 , v y = 0, p = 1, (82)
ρ = 1 + A sin (φ(t)) , v x = v 0 , v y = 0, p = 1. (83)
The function φ(t) is defined as follows

φ(t) = 2π (x cos θ + y sin θ) -(v x cos θ + v y sin θ)t . ( 84
)
The discretization of the computational domain satisfies N y = 2N x , where N x , and N y are, respectively, the number of elements in x-and y-direction. Periodic boundary conditions are imposed. In the table 4 are shown the convergence rates of the rest-mass density for the SRHD equations solved with the DGSEM and the hybrid DG/FV methods. The simulation time is t = 1, and for polynomials of degree N = 3, 5, 7. Observe that the EOC is achieved by the pure DGSEM, but the hybrid DG/FV scheme with all elements flagged as troubled is third order accurate, the order of the WENO3 reconstruction.

One-dimensional Riemann problems

In this section we will present four one-dimensional Riemann problems for the SRHD equations. These problems comprise the standard testbench for SRHD [79,80]. Like for the MHD equations, for all test cases in this section we will assume that the one-dimensional domain is the interval [0, 1]. The membrane is localized at x m = 0.5. The numerical solution computed with the hybrid DG/FV scheme with polynomials of degree N = 2, 4, 6 is compared with the exact solution obtained from the exact Riemann solver by Rezzolla and Zanotti [81]. The left and right vector states in primitive variables are written as w L = ρ, v x , v y , p L and w R = ρ, v x , v y , p R , respectively.

Riemann problem 1 (RP-1). The initial condition is given by the following left and right states w L = 10, 0, 0, 40/3 , w R = 1, 0, 0, 0 .

The EOS corresponds to the ideal gas EOS with adiabatic index γ = 5/3. This test problem is also known as the mildly relativistic blast wave [35,82,83]. The zero pressure in the right state is approximated by p ≈ 2/3 × 10 -6 for numerical reasons. Due to this very small value of the pressure, any oscillation generated in the approximate solution of the pressure will break up the simulation. A very robust shock capturing scheme is then required for handling such kind of problems. Plots of the rest-mass density, pressure, and the x-component of the velocity are depicted in the figure 14 (top). The domain has been discretized with 500 elements, and a polynomial of degree N = 5 is used in the The fluid behind the shock wave is traveling with a lightly relativistic speed v = 0.72c to the right. In this dense shell behind the shock, the fluid is compressed. The fluid is thermodynamically relativistic because it is heated in such a way that the internal energy is much larger than the rest-mass energy, but mildly relativistic dynamically. The calculations are oscillations-free. This shows the robustness of the shock-capturing method.

Riemann problem 2 (RP-2). This problem, also called the highly relativistic blast wave, has the following initial state The ideal gas EOS with adiabatic index γ = 5/3 has been employed. The interval is discretized with 500 elements. Dirichlet boundary conditions have been utilized at domain edges. The final simulation time is t = 0.4. Plots of the rest-mass density, pressure, and the x-component of the velocity are depicted in the figure 14 (second row, from top to bottom). After the breakup of the membrane, a dense intermediate state is formed, and it is located between the right-traveling shock wave and the rarefaction wave moving to the left. Observe that the shock wave and the contact discontinuity are very close to each other. The very thin and dense shell is rather well resolved thanks to the subcell resolution of the hybrid DG/FV method. It is noticeable from the plots the lack of oscillations.

w L = 1, 0, 0, 10 3 , w R = 1, 0, 0, 10 -2 . ( 86
Riemann problem 3 (RP-3). The transverse blast wave problem has initial conditions similar to the RP-2, except that the transverse velocity in the right state is non-zero. The initial left and right states are

w L = 1, 0, 0, 10 3 , w R = 1, 0, 0.99, 10 -2 . ( 87
)
We have employed the ideal gas EOS with adiabatic index γ = 5/3. The domain has 500 elements, and the final simulation time is t = 0.4. Plots of the rest-mass density, pressure, and the x-component of the velocity are depicted in the figure 14 (second row, from bottom to top). It worths to mention that in Newtonian hydrodynamics the transverse momentum is not coupled with the longitudinal one and it thus is simply advected. The momentum equations in SRHD are coupled to each other through the Lorentz factor. This introduces new physical effects, especially with non-zero transverse velocities [84,85]. Observe that the solution is in very good agreement with the exact solution, and the computed solution is oscillation-free.

Riemann problem 4 (RP-4). In the planar shock reflection an ideal cold fluid collides a wall. Afterwards a shock wave propagates backwards, leaving the gas behind at rest. The initial state is given by w = 1, 0.99999, 0, 0.01

The reflecting wall is located in x = 1, and at x = 0 transmissive boundary conditions are imposed. An ideal gas EOS with adiabatic index γ = 4/3 is used. The domain is discretized with 500 elements. The final simulation time is t = 1.5. The exact solution of this problem was first obtained by Blandford and McKee [86]. Plots of the restmass density, pressure, and the x-component of the velocity computed with the hybrid DG/FV scheme are depicted in the figure 14 (bottom). Close to the wall (x = 1), the numerical solution of the rest-mass density does not show an undershooting, as it occurs in the solution computed with finite-volume or finite difference methods [36]. The undershooting is due to the wall heating phenomenon [87]. The shock has been very well captured with the hybrid scheme, but some small oscillations are still visible around it.

Two-dimensional Riemann problems

Two-dimensional Riemann problems are now considered. The initial states have been taken from [38,88]. The computational domain is Ω = [-1, 1]×[-1, 1], and it has been divided into four quadrants. The separating membranes are located along the axes x and y. In these Riemann problems are present all essential features typically found in a two-dimensional flow, for instance, shock reflections, shock interactions, vortices, etc. We write the initial states for every quadrant (Q1, Q2, Q3, and Q4) Q1 := (x, y) ∈ [-1, 1] 2 x ≥ 0, y ≥ 0 , Q2 := (x, y) ∈ [-1, 1] 2 x < 0, y ≥ 0 , Q3 := (x, y) ∈ [-1, 1] 2 x < 0, y < 0 , Q4 := (x, y) ∈ [-1, 1] 2 x ≥ 0, y < 0 .

(89)

For all Riemann problems, an ideal gas EOS with adiabatic index γ = 5/3 is used, and transmissive boundary conditions are used in all faces of the computational domain. Riemann problem 5 (RP-5). The domain is subdivided into four quadrants, where the states in primitive variables are given by w Q1 = 0.035145216124503, 0, 0, 0.162931056509027 , w Q2 = 0.1, 0.7, 0, 1 , w Q3 = 0.5, 0, 0, 1 , w Q4 = 0.1, 0, 0.7, 1 .

(90)

The final simulation time is t = 0.8. In the figure 15 are depicted contour plots of the rest-mass density (along the mesh) and pressure. This problem features contact discontinuities and shock waves. After the breakup of the membranes, two contact discontinuities are originated on the left and bottom of the domain, and in the first quadrant two curved front shocks appear. In the third quadrant a like-jet structure moving in south-west direction is present [88]. The hybrid DG/FV scheme is quite robust and allows to compute in a multidimensional fashion such complex flow patterns. The approximated solution computed with the hybrid DG/FV scheme using polynomials of degree N = 5 on a mesh of 120 × 120 elements in good agreement with the reported in the literature [40,88].

Riemann problem 6 (RP-6). The initial condition is given by the following states in primitive variables w Q1 = 0.5, 0.5, -0.5, 5 , w Q2 = 1, 0.5, 0.5, 5 , w Q3 = 3, -0.5, 0.5, 5 , w Q4 = 1.5, -0.5, -0.5, 5 .

(91)

The final simulation time is t = 0.8. In the figure 16 are shown contour plots of the rest-mass density and the Lorentz factor. The mesh with the FV subcells is also shown. From the plots can be seen the interaction of four vortex sheets, which form a spiral with very low rest-mass density in the center of the domain. The hybrid DG/FV scheme using polynomials of degree N = 5 on a mesh of 120 × 120 elements has an excellent behavior in the very low density region, which means that the scheme does not produce negative densities or pressures even in regions with cavitation. Riemann problem 7 (RP-7). In this two-dimensional Riemann problem, the initial condition is given as follows w Q1 = 1, 0, 0, 1 , w Q2 = 0.5771, -0.3529, 0, 0.4 , w Q3 = 1, -0.3529, -0.3529, 1 , w Q4 = 0.5771, 0, -0.3529, 0.4 .

(92)

The final simulation time is t = 0.8. In the figure 17 are shown contour plots of the rest-mass density and the Mach number computed with a hybrid DG/FV scheme using polynomials of degree N = 5 on a mesh of 120 × 120 elements. In this problem, the planar rarefaction waves originated after the breakup of the membranes interact each other. This interaction generates two symmetric shock waves in the zone where the rarefaction waves have interplayed. Observe the FV subdomains in and around the regions where the rarefaction waves are present.

Cloud-shock interaction

In this problem, we follow the initial setup given by He and Tang [78]. An ideal gas EOS with adiabatic index γ = 5/3 is employed in the simulation. The computational domain is Ω = [0, 2]×[0, 1]. Dirichlet boundary conditions are imposed. The center of the cylindrical cloud is located in the point (1.4, 0.5) and its radius is r = 0.15. The initial state for the cloud in primitive variables is w c = (3.1538, 0, 0, 0.05). The left traveling shock wave at time t = 0 is located at x = 1.6, with left and right states given by w L = (1, 0, 0, 0.05), w R = (1.86522508063118, -0.19678110737829, 0, 0.15).

(93)

We discretize the computational domain with 160 × 80 elements. In the figure 18 are depicted contour plots of the rest-mass density and Mach number at simulation time t = 3.0. The mesh and the FV subdomains are also plotted. The computations were done with the hybrid DG/FV scheme with polynomials of degree N = 5. The profile of the bubble after the interaction with the shock wave is in good agreement with the reported in the literature [36,78]. The computation was stable and negative densities or pressures were not generated.

hybrid DG/FV scheme is an excellent candidate as a numerical solver for equations systems typical in computational astrophysics. Among other outstanding properties of the hybrid scheme we can mention: the method is highly parallelizable [START_REF] Sonntag | [END_REF], allows unstructured quadrilateral/hexahedral meshes, and it is very robust when strong shocks take place in the simulation. Further research is focused on hybrid DG/FV schemes with quad-/octree based adaptive mesh refinement in order to simulate in a more efficient manner small structures in magnetohydrodynamic and relativistic turbulence simulations. Furthermore, the proposed method can also be used in the context of polygon/polyhedron computational meshes. Some similar work has been done in that direction by other researchers [89]. This will be useful for some CFD codes for astrophysics which make use of DG schemes on moving unstructured grids defined as the Voronoi tessellation of a set of mesh-generating points, which will require such efficient shock-capturing schemes.

Figure 1 :

 1 Figure 1: Mapping from the physical element Ω i to the reference element [-1, 1] 3 .

Figure 2 :

 2 Figure2: Mapping from the physical element Ω i to the reference element [-1, 1] 2 in a two-dimensional domain. The blue-filled circles depict the two-dimensional Gauss-Legendre quadrature points in the reference element, where the solution is computed, and the red-filled squares represent the nodes where the boundary fluxes are calculated, and they are located at Gauss-Legendre quadrature points along the boundary.

Figure 4 :

 4 Figure4: Hybrid DG/FV. Flagged element surrounded by DG elements (left). The DG polynomial in the troubled cell is then projected onto a FV subdomain made of (N + 1) 2 subcells. For those subcells lying at the FV subdomain edges, neighboring elements are also projected onto a FV subdomain in order to use the cell averages at edges for computing the reconstructed values at subcell interfaces (right).

Figure 5 :

 5 Figure 5: Flow diagram for the hybrid DG/FV method.

  (bottom) are depicted plots of the density, thermal pressure, and the y-component of the magnetic field at simulation time t = 0.1. The flow pattern in the figures, from left to right, are a left-going fast

Figure 6 :Figure 8 :

 68 Figure6: Plots for the Riemann problems RP-1, RP-2, RP-3, and RP-4. The initial states have been provided in the table 3. An ideal gas EOS with adiabatic index γ = 5/3 is used. The hybrid DG/FV scheme used for solving these problems employs polynomials of degree N = 2, 4, 6.

Figure 9 :

 9 Figure 9: Orszag-Tang vortex problem. Contour plots of the density (left) and the pressure (right) at time t = 0.5. The calculations were performed with a hybrid DG/FV scheme, with polynomial degree N = 5. The computational domain is Ω = [0, 1] × [0, 1], and it has been discretized by using a mesh of 120 × 120 elements.

xFDx

  Pressure, t = 0.5, y = 0.3125 Pressure, t = 0.5, y = 0.4277 FD-W5-600 2 DG-N5-1202

Figure 10 :

 10 Figure 10: Orszag-Tang vortex problem. Plots of the pressure along the slices y = 0.3125 (left) and y = 0.4277 (right) at time t = 0.5. The computations were done with a very high-order finite difference WENO5 scheme and the hybrid DG/FV method with N = 5.

Figure 11 :

 11 Figure 11: Cylindrical blast wave. Contour plots of the density (left) and pressure (right) at time t = 0.1 obtained with a hybrid DG/FV scheme, with polynomial degree N = 5. The computational domain is Ω = [0, 1] × [0, 1], and it is discretized by using a mesh of 120 × 120 elements.

x 5 FDy 5 FD

 55 Magnetic Field Y, t = 0.25, y = 0.Magnetic Field X, t = 0.25, x = 0.

Figure 13 :

 13 Figure13: Rotor problem. Plots of the y-component of the magnetic field along the slice y = 0.5 (left) and plots of the x-component of the magnetic field along the slice x = 0.5 (right) at time t = 0.25. The computations were done with a very high-order finite difference WENO5 scheme and the hybrid DG/FV method with N = 5.

)Figure 14 :

 14 Figure 14: One-dimensional Riemann problems. Plots for the Riemann problems RP-1, RP-2, RP-3, and RP-4. The computations were performed with the hybrid DG/FV scheme with polynomials of degree N = 2, 4, 6, and with WENO3 reconstruction operator in the FV solver.

Figure 15 :

 15 Figure 15: Two-dimensional Riemann problem RP-5. Contour plots of the rest-mass density (left) and pressure (right) at t = 0.8. The solution was obtained with a hybrid DG/FV scheme with polynomials of degree N = 5. The domain Ω = [-1, 1] × [-1, 1] was decomposed into 120 × 120 elements.

Figure 16 :

 16 Figure 16: Two-dimensional Riemann problem RP-6. Contour plots of the rest-mass density (left) and Lorentz factor (right) at t = 0.8. The solution was obtained with a hybrid DG/FV scheme with polynomials of degree N = 5. The domain Ω = [-1, 1] × [-1, 1] was decomposed into 120 × 120 elements.

Table 3 :

 3 Riemann problems for the magnetohydrodynamics equations. The initial states have been taken from Ryu and Jones[28]. The adiabatic index is set to γ = 5/3 and the computational domain is the interval [0, 1]. For the problems 1 and 2, the magnetic field components have to be divided by m is the position of the initial discontinuity, and u L and u R are two constant states written in conservative variables. Typically, the initial condition in a hyperbolic boundary-value problem is provided in terms of the primitive variables[START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. Therefore, in this work, the left state will be represented by w L = ρ, v x , v y , v z , p, B x , B y , B z L and the right state by w R = ρ, v x , v y , v z , p, B x , B y , B z R . These vectors are written in primitive variables, in contrast with the vectors u L and u

	√
	4π.
	Test	State	ρ	v x	v y	v z	p	B x	B y	B z
	RP-1	w L w R	1.00 1.00	10.00 -10.00	0.00 0.00	0.00 0.00	20.00 1.00	5.00 5.00	5.00 5.00	0.00 0.00
	RP-2	w L w R	1.08 1.00	1.20 0.00	0.01 0.00	0.50 0.00	0.95 1.00	2.00 2.00	3.60 4.00	2.00 2.00
	RP-3	w L w R	1.00 0.30	0.00 0.00	0.00 0.00	0.00 1.00	1.00 0.20	0.70 0.70	0.00 1.00	0.00 0.00
	RP-4	w L w R	1.00 0.125	0.00 0.00	0.00 0.00	0.00 0.00	1.00 0.10	0.75 0.75	1.00 -1.00	0.00 0.00
	where x

R . For the four Riemann problems discussed, the one-dimensional domain is the closed interval [0, 1].

Table 4 :

 4 Convergence rates for the SRHD equations with initial condition given by the smooth flow problem. Results for L 2 norm error of the rest-mass density are given. The computations were done with the pure DGSEM, and the hybrid DG/FV with all elements flagged as troubled. The schemes use polynomial degrees N = 3, 5, 7, that is for orders of accuracy O(4), O(6) and O(8). FV scheme. Dirichlet boundary conditions are set at both sides of the interval. The final simulation time is t = 0.4. The waves present in this problem are a left-going transonic rarefaction wave, a contact discontinuity and a right-going shock wave.

	DGSEM	Hybrid DG/FV

Acknowledgments

This research was funded by the Deutscher Akademischer Austauschdienst under grant A0873943. Further support was provided by the Deutsche Forschungsgemeinschaft, grant MU 1319/12-2. We thank Michael Dumbser for the very careful reading of the manuscript and their valuable suggestions, which improved the quality of this paper. Many thanks to Luciano Rezzolla, and Olindo Zanotti for providing the exact solution of the one-dimensional SRHD Riemann problems discussed in this work.

rarefaction, a left-going slow compound wave, a contact discontinuity, a right-going slow shock wave, and a rightgoing fast rarefaction wave. The hybrid DG/FV scheme with WENO3 reconstruction produces a solution that is in good agreement with the exact solution, except in the compound wave and contact discontinuity.

Multidimensional test problems 4.4.1. Magnetic field loop advection

The high-order DGSEM is not able by itself to maintain the solenoidal constraint ∇ • B = 0. It requires an additional procedure to tackle this issue. Here we employ the divergence cleaning of Dedner et al. [START_REF] Dedner | [END_REF], with c r = 0.18, and c h determined by the maximum propagation speed in the system. The magnetic field loop advection is designed to assess the capability of the numerical scheme to keep the solenoidal constraint. The problem consists in a cylindrical current distribution, which is advected along some direction of the computational domain. The main feature of the loop is that it remains in magnetostatic balance, and therefore, after some periods, its profile should be kept. Without a divergence cleaning algorithm, the magnetic loop will smear over the time [69].

The problem setup is as follows: the computational domain is Ω = [-1.0, +1.0] × [-0.5, +0.5]. The density and the pressure are initially set to ρ = 1 and p = 1 in the whole domain. The initial velocity profile is given by

where v 0 = √ 5 is the magnitude of the velocity, and n x = 2/ √ 5 and n y = 1/ √ 5 are the components of the unit vector in the loop travel direction. The magnetic field is uniform everywhere, except for the loop structure of radius R = 0.3. For r ≤ R we have

where r = x 2 + y 2 , and B 0 = 10 -3 is the magnitude of the magnetic field, and it is chosen in such a way that the magnetic pressure is smaller than the gas pressure. An ideal gas EOS with adiabatic index γ = 5/3 is used. Regarding the discretization, the computational domain is made of 180 × 90 elements. A polynomial degree N = 5 is employed in the hybrid DG/FV scheme. We set periodic boundary conditions at the boundaries of the computational domain.

The simulation time is t = 2, and in the figure 7 is depicted the magnetic pressure at times t ∈ 0.0, 2.0 . Observe that the loop profile is only kept when the divergence cleaning is used in the algorithm.

Conclusions

In this work we have presented the buildings blocks of a high-order hybrid DG/FV scheme for MHD and SRHD. The hybrid scheme is based on a discontinuous Galerkin spectral element on quadrilateral/hexahedral meshes as main spatial discretization operator, and it is employed basically in regions with smooth flows. For those regions containing discontinuities, a robust second/third order finite volume scheme is used as shock capturing approach. In order to handle discontinuities, an oscillations indicator is used to detect and to mark the elements containing them. These troubled elements are projected onto a subdomain made of FV subcells, and then evolved with the robust FV method. In this way the very high-order hybrid method is capable of solving problems with shock waves. Time discretization is performed with a fourth-order SSPRK method. The solenoidal constraint is maintained with the GLM divergence cleaning. From our experience in the area of very high-order FV-WENO schemes for computational astrophysics [33,36], the DGSEM domains are evolved using the conservative variables, while in the FV subdomains the solution is computed by employing the primitive variables. This choice provides more stability to the computations.

A comprehensive testbench for MHD and SRHD was discussed. One-and two-dimensional test problems were presented. Convergence tests have been also discussed in order to check the order of accuracy of the hybrid scheme for both equations systems. The computations show that the hybrid scheme converges to the theoretical order of accuracy. For problems with shocks, the hybrid scheme captures very well the discontinuities, even in very cases with extreme density or pressure conditions. Polynomials up to sixth order degree were used. These results demonstrated that the