In this report, we address several aspects of the approximation of the MHD equations by a Galerkin Discontinuous finite volume schemes.

1 Some properties of the MHD system

Equations

The Magneto-Hydro-Dynamics (MHD) equation are a useful model for describing the behavior of a compressible conductive fluid. The unknowns are the fluid density ρ, the velocity u ∈ R 3 , the internal energy e, the pressure p and the magnetic field B ∈ R 3 . All the unknowns depend on the space variable x ∈ R 3 and the time variable t.

The equations read

    ρ ρu B Q     t + ∇ •     ρu ρu ⊗ u + (p + B•B 2 )I -B ⊗ B u ⊗ B -B ⊗ u (Q + p + B•B 2 )u -(B • u)B     = 0, Q = e + u • u 2 . (1)
The notation I stands for the 3 × 3 identity matrix. The pressure is related to the internal energy e and the density ρ by a pressure law. In this document, we shall only consider the perfect gas law with a constant polytropic exponent γ. It reads p = P (ρ, e) = (γ -1)ρe, γ > 1.

(

) 2
The previous equations are supplemented by the following divergence condition on the magnetic field

∇ • B = 0. (3)
The divergence free condition on the magnetic field is very important for physical reasons: it ensures that there is no magnetic charge. This condition is difficult to express on the numerical side. Therefore some authors [7], [START_REF] Dedner | Hyperbolic divergence cleaning for the MHD equations[END_REF] have suggested to extend the ideal MHD system in the following way

      ρ ρu B Q ψ       t + ∇ •         ρu ρu ⊗ u + (p + B•B 2 )I -B ⊗ B u ⊗ B -B ⊗ u + ψI (Q + p + B • B 2 )u -(B • u)B c 2 h ∇ • B         = 0, Q = e + u • u 2 . (4)
We have added a new unknown ψ whose role is to "clean" the divergence of the solution. Actually, the divergence perturbations are convected in the computational domain at the constant velocity c h . With adequate boundary conditions, the perturbation will be damped. The velocity c h can be chosen arbitrarily. In practice, it has to be higher than the highest wave speed of the original MHD system.

We observe that if ∇ • B = 0 and ψ =Cst, then the modified system (1) is equivalent to the MHD system [START_REF] Dedner | Hyperbolic divergence cleaning for the MHD equations[END_REF].

The two above systems can be put in a conservative form (with the space dimension d = 3)

w t + d i=1 f i (w) xi = 0 (5)
We shall make use of the Einstein summation convention on the repeated indices and also write

∂ t w + ∂ i f i (w) = 0 (6)
For any vector n = (n 1 , n 2 , n 3 ) ∈ (R 3 ), the vector

f (w, n) = d i=1 f i (w)n i (7)
is called the flux vector.

Hyperbolicity

In order to study the hyperbolicity of the system, we first write its one-dimensional form by supposing that all the data do not depend on x 2 and x 3 . The magnetic field vector can then be split into a normal component and a tangential one

B 3D = (B x , B y , B z ), b = B x , B 2D = (B y , B z ). ( 8
)
Due to the divergence free condition, its normal component b > 0 is a constant parameter. Only the tangential part B 2D is varying.

In the same way, the velocity has a normal and a tangential part u = (u x , u y , u z ), u = u x , v = (u y , u z ).

With a slight change of notations, the one-dimensional system reads ρ t + (ρu) x = 0, (ρu) t + (ρu 2 + p + 1 2 B 2 ) x = 0, (ρv) t + (ρuv -bB) x = 0,

B t + (uB -bv) x = 0, E t + ((E + p + 1 2 B 2 )u -bB • v) x = 0, E = p γ -1 + 1 2 ρ(u 2 + v 2 ) + 1 2 B 2 .
(10)

The one-dimensional conservative variables are w = (ρ, ρu, ρv, B, E),

and the primitive variables are Y = (ρ, u, p, v, B).

In the primitive variables, the system becomes

∂ t Y + C(Y)∂ x Y = 0, C(Y) =           u ρ u 1 ρ B2 ρ B3 ρ γp u u -b ρ u -b ρ B 2 -b u B 3 -b u           . ( 13
)
We note a the sound speed, c A the Alfvèn speed, c f,s the fast and slow magnetoacoustic speeds and λ i the eigenvalues sorted by increasing order.

c A = b √ ρ , a = γp ρ , c f,s = 1 2 b 2 + B 2 ρ + a 2 ± 1 4 b 2 + B 2 ρ + a 2 2 - a 2 b 2 ρ , λ 1 = u -c f λ 2 = u -c A λ 3 = u -c s λ 4 = u λ 5 = u + c s λ 6 = u + c A λ 7 = u + c f (14)
The eigenvalues are all real and thus the one-dimensional MHD system is hyperbolic. But the MHD equation are also invariant in a rotation, thus the full MHD system is also hyperbolic. It is possible to compute the eigenvectors analytically:

r 1 = (-ρ, c f , -γp, c f B 2 b 1 -c 2 f /c 2 A , c f B 3 b 1 -c 2 f /c 2 A , B 2 c 2 A /c 2 f -1 , B 3 c 2 A /c 2 f -1 ) r 2 = (0, 0, 0, B 3 , -B 2 , √ ρB 3 , - √ ρB 2 ), r 3 = (-ρ, c s , -γp, c s B 2 b (1 -c 2 s /c 2 A ) , c s B 3 b (1 -c 2 s /c 2 A ) , B 2 c 2 A /c 2 s -1 , B 3 c 2 A /c 2 s -1
), r 4 = (1, 0, 0, 0, 0, 0, 0),

r 5 = (-ρ, -c s , -γp, -c s B 2 b (1 -c 2 s /c 2 A ) , -c s B 3 b (1 -c 2 s /c 2 A ) , B 2 c 2 A /c 2 s -1 , B 3 c 2 A /c 2 s -1
),

r 6 = (0, 0, 0, B 3 , -B 2 , - √ ρB 3 , √ ρB 2 ), r 7 = (-ρ, -c f , -γp, -c f B 2 b 1 -c 2 f /c 2 A , -c f B 3 b 1 -c 2 f /c 2 A , B 2 c 2 A /c 2 f -1 , B 3 c 2 A /c 2 f -1
).

(15)

Entropy and Mock theory

The construction of the Discontinuous Galerkin approximation for the MHD system enjoys nice entropy properties. In order to state this entropy dissipation property, we first recall some notions on entropy and hyperbolic systems of conservation laws.

Shocks, entropy

It is well known that the solutions of the system

∂ t w + ∂ i f i (w) = 0 (16)
can become discontinuous in a finite time even if the initial condition is very smooth. We thus have to define a notion of weak solution. Let v be a test vector function in

D(R + × R d ) m (v is not necessarily zero at time t = 0). Let w 0 be the initial condition. A weak solution w in L ∞ (R + × R d ) m satisfies t 0,x -wv t -f i (w)∂ i v + t=0,x w 0 v = 0. ( 17
)
If w is of class C 1 in R + ×R d but on space-time surfaces, then w is a classical solution where it is smooth. On a surface of discontinuity with a normal vector n = (n x , n t ), we note [w] the jump of w. A weak solution w satisfies the Rankine-Hugoniot conditions

[w] n t + f i (w) n x i = 0 (18)
The weak solutions are generally not unique. A supplementary criterion helps to select a solution. A classical criterion is the Lax entropy criterion. To write it, we need first a supplementary conservation law that we write

∂ t S 0 (w) + ∂ i S i (w) = 0. ( 19
)
This supplementary PDE holds true when

∇S 0 • ∇f i = ∇S i (20)
If S 0 is also strictly convex then it is called a Lax entropy of the system of conservation laws. The function S i are the components of the entropy flux. We require that the weak solutions satisfies the supplementary inequality

S 0 (w) t + ∂ i S i (w) 0. ( 21
)
A particular role is then played by the so-called entropy variables

w = ∇S 0 (w). ( 22
)
We make this role more precise in the next section.

Legendre transform

A convex function S being given, the Legendre transform S is defined by

S(w) = max w (w • w -S(w)) (23)
When everything is smooth, the maximum is reached at a point w such that

w = ∇ w S(w) (24)
It defines indeed an admissible change of variables because the jacobian of the transformation is invertible (it is the hessian matrix of S)

Furthermore, it can be proved that the Legendre transform is an involution. Let us prove it in the smooth case. The gradient of S(w) is

∇ w S(w) = w • w ′ (w) + w(w) -∇ w S(w(w))w ′ (w) (25) Then ∇ w S(w) = w ⇔ w = ∇ w S(w) (26)
In other words, the gradient of S defines the inverse change of variables. We deduce

S(w) = w • w -S(w) with w = ∇ w S(w) S(w) = w • w -S(w) with ∇ w S(w) = w S(w) = w • w -S(w) with ∇ w S(w) = w (27)
In the last formula, we recognize the Legendre transform of S. Thus, we have proved that S = S.

Mock theorem

A system of conservation laws is symmetrizable if it is possible to find a change a variables, such that, in these news variables, the convection matrix is symmetric. The Mock theorem ensures that a system is symmetrizable iff it admits an entropy. Let us prove it. Suppose that the system admits an entropy S 0 . Let us note w = ∇ w S(w). According to the previous section, the inverse change of variables is given by the Legendre transform of S. Furthermore, we define an analog of the Legendre transform of the flux thanks to the formula

S k (w) := f k (w(w)) • w -S k (w(w)). ( 28
)
It is an abuse of notations because S k is not the Legendre transform of S k as defined in formula (23). We can then verify that, as for the entropy S 0 , we have

∇ w S k (w) = f k (w(w)) (29)
We deduce that in the entropy variables w the system becomes

d k=0 ∇ 2 w S k ∂ k w = 0 ( 30
)
and it is indeed a symmetric system.

In order to prove the reverse, we use the Poincaré lemma: the jacobian matrices of w(w) and g(w) = f (w(w)) are symmetric and thus w(w) and g(w) are the gradient of functions S(w) and S k (w). Using the same computations we see that S = S is indeed an entropy. The MHD system does not enter exactly this framework. We have to modify the approach and take into account the divergence condition on the magnetic field B. See Section 1.3.4

Symmetric form of the MHD system

It is possible to symmetrize the MHD system by modifying the initial equations with terms that contain only the divergence of B [START_REF] Barth | On the role of involutions in the discontinous Galerkin discretization of Maxwell and magnetohydrodynamic systems[END_REF].

For this, we introduce the (physical) entropy of the fluid

s = ln p ρ γ (31)
This entropy satisfies the following PDE

s t + u • ∇s + (γ -1) u • B p ∇ • B = 0 (32)
Combining with the mass conservation law, we find

(ρs) t + ∇ • (ρus) + (γ -1) ρ u • B p ∇ • B = 0 (33)
If the divergence of B is zero, the quantity S 0 = ρs satisfies a supplementary conservation law. It is also possible to show that S 0 is convex with respect to the conservative variables. But it is not a Lax entropy because the initial system does contain the condition ∇ • B = 0: it comes from the initial condition.

Besides, the change of variables w = ∇ w S 0 (w) does not symmetrize the MHD system.

In order to find a symmetrization, we write the system under the form

∂ t w + ∂ i f i (w) = 0, ∂ i B i = 0. ( 34
)
with the entropy condition

∂ t S 0 + ∂ i S i 0 (35)
We can add to the MHD equations some combinations of ∇ • B. Those combinations are given by a function Λ(w) where w are the symmetry variables

∂ t w + ∂ i F i (w) + ∂ i B i ∇ w Λ(w) = 0, (36)
Simple computations show that the system becomes symmetric with the change of variables

w = ∇ w S 0 (37)
The fluxes are given by

f i = ∇ w S i -B i ∇ w Λ(w) (38)
and we obtain generalized duality relations

S 0 = w∇ w S 0 -S 0 S i = w∇ w S i -S i (39)
In the considered case, we find

Λ = (γ -1)ρ u • B p ( 40
)
This function is homogeneous of degree one with respect to w, which leads to

Λ = ∇ w Λ • w (41)
This is a constructive way to write the modified MHD system of Powell [9].

It would be also interesting to find the entropies of the modified system (4).

2 Numerical resolution of the Riemann problem of the MHD

General resolution

In this section, we describe the resolution of the Riemann problem for the MHD.

We chose an approach that does not use particular properties of the MHD system and thus can be extended to other systems of conservations laws. We have first to recall some basic notions on the Riemann problem. A field i is linearly degenerated (LD) iff

∇λ i (Y) • r i (Y) = 0 ( 42
)
A field is genuinely non-linear (GNL) iff at all vector

Y ∇λ i (Y) • r i (Y) = 0. ( 43
)
It is also possible to suppose

∇λ i (Y) • r i (Y) > 0 (44)
(if it is not the case, change r i in -r i ). For the MHD system, the fields 2, 4 and 6 are LD while the other fields are GNL. the field 4 is a contact discontinuity, in which only the density ρ jumps. For the fields 2 and 4, the Riemann invariants are ρ, u, p, B 2 and the two components of ∓B + √ ρv. An important fact is that if B = 0 then the eigenvalues 1, 2 and 3 merge as the eigenvalues 5,6 and 7. The system is still diagonalizable. On the other hand, it is possible that ∇λ i (Y) • r i (Y) = 0 for some vector Y (non-convexity, defect of Genuine Non-Linearity). Thus, the Lax theorem does not apply for all data and the uniqueness of the solution to the Riemann problem is no more ensured.

To each GNL field, we can generally associate particular solutions: the shocks and the simple waves. Let us start to recall how to construct a i-simple wave. We take a left state w L and we want to join it to a state w depending on one single parameter η. For this, we consider the following ordinary differential equation

Y ′ (η) = r i (Y(η)), (45)
with the initial condition

Y ′ (η) = r i (Y(η)), (46)
The solution of this problem depends of course on the chosen normalization of the eigenvector r i . It defines a curve in the phase space R 7 . We can compute the variation of the eigenvalue along the curve by solving

ξ ′ (η) = ∇λ i (Y(η)) • Y ′ (η) = ∇λ i (Y(η)) • r i (Y(η)), ξ(η 0 ) = λ i (Y L ). ( 47
)
According to the GNL hypothesis, the change of variables ξ = ξ(η) is monotone and thus locally bijective. Let us define

Z(x, t) = Y(ξ -1 (x/t)). ( 48
)
The vector function Z is indeed a solution to Z t + A(Z)Z x = 0. We have just constructed the i-simple wave. We shall note it

Y = D i (Y L , η i ) D i (Y L , η 0,i ) = Y L . (49)
In practice, the choice of the parameter η is important. A bad choice would lead to complicated computations. For example, it is not possible to chose η = R where R is a Riemann invariant (satisfying ∇R • λ i = 0). For theoretical purposes, it is convenient to take: η = λ i , the natural normalization of the eigenvector is then ∇λ i • r i = 1. For the numerical resolution, we propose to take one component (number k) of Y (which is not a Riemann invariant). The eigenvector r i should be normalized in such way that r i,k = 1. In the following presentation we will suppose that ξ(η) is an increasing function and thus that the admissible part of the simple wave curve corresponds to η > η 0 . In practice, we can have to revert this condition, for example if the chosen parameter is the density ρ: in such a case, we must exchange the left state Y L and the right state Y R or decide that the admissible part of the curve is η < η 0 .

The shocks are more difficult to parameterize. For a shock of speed s, the Rankine-Hugoniot relations read

s(η)(w(η) -w L ) = f (w(η)) -f (w L ), w(η 0 ) = w L . (50)
Deriving this relation with respect to the parameter η, we find

s(η 0 )w ′ (Y(η 0 ))Y ′ (η 0 ) = f ′ (w L )w ′ (Y(η 0 ))Y ′ (η 0 ) (51) multiplying by w ′ (Y(η 0 )) -1 we obtain s(η 0 )Y ′ (η 0 ) = A(Y L )Y ′ (η 0 ). ( 52
)
Thus s(η 0 ) is an eigenvalue of A and Y ′ (η 0 ) a corresponding eigenvector. Locally, in the non degenerated cases, we will then find m shock curves C i , i = 1 • • • m tangent to the simple wave curves D i at η = η 0 . In addition, if the parameter η is a component of Y or the wave speed λ i (Y) the reunion of the shock and simple wave curves is of class

C 2 .
The problem is that it is difficult to order globally the shock curves because the curves C i maybe tangent at some point.

In the regular case, the usual criterion is to say that the shock belong to the i-th family iff

λ i (Y L ) > s(η) > λ i (Y(η)). ( 53
)
The interpretation is that the i-th characteristic curve coming from the left and the i-th characteristic curve coming from the right must impinge the shock.

In this way, it is generally possible to build the shock curves

Y(η) = C i (Y L , η), Y(η 0 ) = Y L . ( 54
)
We then stick those two types of solutions and introduce the mixed curves M i . The mixed curve Y = M i (Y L , η) permits to find all the states Y that can be connected to a left state Y L by a shock or a simple wave of the family i

M i (Y L , η) = D i (Y L , η) si η > η i 0 , C i (Y L , η) si η < η i 0 . (55)
The curves are of class C 2 with an adequate choice of the parameter η i .

In the situation where the field is LD there is no more distinction between the shocks and the simple waves. It is no more possible to take the wave speed as a parameter. But the method of construction is very similar. For a LD wave, we write

Y ′ (η) = r i (Y(η)), Y(0) = Y L , Y(η) = M i (Y L , η). (56)
We are now in a position to solve the Riemann problem. Starting from two states Y L and Y R . The numerical problem is to find

m parameters η 1 , • • • , η m such that M m (• • • M 2 (M 1 (Y L , η 1 ), η 2 ) • • • , η m ) = Y R . (57)
Let us note that it is not a numerical trivial problem because in some cases it is possible to find several solutions satisfying the Lax criterion. See [11] for the theory and [10] for the numerical consequences.

Shock curves construction

The Rankine-Hugoniot relations for the MHD can be written

m 2 [τ ] + p + 1 2 B 2 = 0, m [v] -b [B] = 0, m 2 [τ B] -b 2 [B] = 0, γ γ -1 pτ + 1 2 m 2 τ 2 + τ - b 2 2m 2 B 2 = 0, τ = 1/ρ. ( 58
)
The Rankine-Hugoniot relations can also be rewritten

m 2 1 ρ - 1 ρ L + p -p L + B 2 -B 2 L 2 = 0, m 1 ρ - 1 ρ L -b (B 2 -B 2,L ) = 0, m 1 ρ - 1 ρ L -b (B 3 -B 3,L ) = 0, m 2 B 2 ρ - B 2,L ρ L -b 2 (B 2 -B 2,L ) = 0, m 2 B 3 ρ - B 3,L ρ L -b 2 (B 3 -B 3,L ) = 0, γ γ -1 p ρ - p L ρ L + 1 2 m 2 1 ρ 2 - 1 ρ 2 L + B 2 ρ - B 2 L ρ L - b 2 2m 2 B 2 -B 2 L = 0, m = ρ(u -s) = ρ L (u L -s). (59)
In order to construct the shock curve i, we suppose that we know the left state Y L and the density ρ of the right state Y = (ρ, u, p, v 2 , v 3 , B 2 , B 3 ). The unknowns are then the 6 remaining components of Y and (m, s) which gives 8 unknowns. Solving the eight equations in (59) should permit to express the unknowns as functions of ρ and Y L . It appears that the system can be put in a polynomial form. Recently, new algorithms have been designed to rigorously solve this kind of system in a formal way [5]. They are implemented for example in Maple. Using this formal resolution, we find that the shock speed s is a root of a polynomial P (s) of degree 6 whose coefficients depend on the components of Y L . We do not give their expressions here because they are rather complicated.

The other unknowns can be expressed as explicit (but complicated) functions of s. The method of resolution is then the following

• Let Y L and ρ be given;

• Compute all the real roots of P (s);

• For each real root, compute the full vector Y;

• Compute the ith wave speed associated to Y L and Y;

• Check if the constructed shock satisfies the Lax characteristic condition for the field i, i.e.

if λ i (Y L ) > s > λ i (Y).
• If it is the case, return the corresponding Y.

It might happen that several solutions are found... The parametrization of the simple waves or the LD waves is simply obtained by solving (46) with a classical fourth order Runge-Kutta algorithm. We could have also used a more algebraical approach based on the Riemann invariants.

Depending on the chosen parameter η i for the ith wave, we have to normalize accordingly the eigenvector. For practical reasons, our chosen parameters are

η = (ρ 1 , α 2 , ρ 3 , ρ 4 , ρ 5 , α 6 , ρ 7 ) (60)
The parameters α 2 and α 6 correspond to the Alfvèn waves. Recall that through these waves, the transverse magnetic field has a constant norm and only rotates around the x axis. It is thus natural that the two parameters α 2 and α 6 are the angle of the rotation of the magnetic field around this axis. Consequently, the eigenvectors 1, 3, 4, 5 and 7 are normalized in such way that their first component is one (their expressions are given in (15)). The eigenvectors 2 and 6 have to be divided by -√ ρ.

Numerical resolution of the non-linear system

The non-linear system (57) has to be solved by a fixed-point algorithm. It is not easy to find the good guess that ensures convergence. In order to improve the robustness, we first write a fixed point method that approaches the Newton algorithm when the two states Y L and Y R are close to each other. For this, we define

F(η) = M m (• • • M 2 (M 1 (Y L , η 1 ), η 2 ) • • • , η m ) -Y R , η = (η 1 • • • η m ). ( 61
)
in such a way that the non-linear system becomes

F(η) = 0. (62)
In order to implement a quasi-Newton method, we have to find an approximation of the jacobian matrix of F. If we suppose that Y L ≃ Y R , it is reasonable to suppose also that η ≃ η 0 . We then find

F ′ (η) ≃ J(η) = ( dM 1 (Y L , η 1 ) dη 1 • • • dM m (Y L , η m ) dη m ) (63) because DM i (Y, η i ) DY ≃ I (64)
where I denotes the m × m identity matrix. If η i < η 0 i is on the side of a simple wave, we find

dM i (Y, η i ) dη i = r i (Y ). ( 65
)
We could use the same formula in the case of a shock wave η i > η 0 i because the left and right states are close and thus η is close to η 0 . But for more generality, we also describe the way to compute dMi(Y,ηi) dηi in the case of a shock wave. For this, we differentiate the Rankine-Hugoniot relations

s(η)(w(Y(η)) -w L ) = f (w(Y(η))) -f (w L ). ( 66
)
We find

s ′ (w -w L ) + sw ′ Y ′ = f ′ w ′ Y ′ (67)
But we also know that one component of Y is equal to η (we can suppose that it is the first). We can then rewrite the previous linear system in the form

sw ′ (Y) -f ′ (w)w ′ (Y) w -w L 1, 0 • • • 0 0 Y ′ s ′ =      0 . . . 0 1      (68)
Its resolution by a standard LU method provides Y ′ (η) and s ′ (η). We are now able to compute the approximation of the jacobian matrix of F given in (63). The quasi-Newton algorithm is

η (n+1) = η (n) -J(η (n) ) -1 F(η (n) ). ( 69
)
In practice, the initial vector guess η (0) has to be chosen carefully or the iterative method does not converge. One way to improve the robustness would be to enrich the method by a continuation method.

Numerical application

We verify that our rough method can reproduce the results found in [11]. We first fix the physical constants to

γ = 5 3 , b = 1.5. ( 70
)
We consider the following Riemann problem

Y L =           3 0 3 0 0 1 0           , Y R =           1 0 1 0 0 cos(1.5) sin(1.5)           (71)
The initial guess is

η (0) = (ρ 1 , α 2 , ρ 3 , ρ 4 , ρ 5 , α 6 , ρ 7 ) = (2.2, 0.7, 2.1, 1.5, 1.1, 0.5, 1) (72)
We obtain a convergence of the algorithm towards the correct parameter vector found in [11].

- --------------------------------------------- -----------------------------------------------Error: 8.68676031519186201E-010

Let us emphasize that more sophisticated approaches exist but they are necessarily based on fine properties of the MHD system. Our numerical approach can be extended to many systems of conservation laws where the Rankine-Hugoniot relations can be written under a polynomial form.

3 Multiple solution of a simplified MHD system Non uniqueness is known to arise in the case of coplanar initial conditions, i.e. when the transverse magnetic field has opposite orientation on each side of the initial condition.

To study this phenomenon, the idea is to consider a simple 3 × 3 system, derived from the MHD model (the derivation is explained in [8]).

w t + F(w) x = 0, w =   u v w   , F(w) =   cu 2 + v 2 + w 2 2uv 2uw   . ( 73
)
This model allows us to get a qualitative description of the interactions between Alfvèn and magnetoacoustic waves. The vector (v, w) stands for the transverse velocity or magnetic field and u is a thermodynamic parameter of the fluid.

Hyperbolicity

The jacobian matrix of the fluxes reads:

A(w) = F ′ (w) =   2cu 2v 2w 2v 2u 0 2w 0 2u   ( 74
)
Because A is symmetric, it is diagonalizable and its eigenvalues are all real and thus the system is hyperbolic. The computation of the eigenvectors and their associated eigenvalues gives:

λ s = (c + 1)u -r λ a = 2u λ f = (c + 1)u + r r s =   1 2 λ s -u v w   r a =   0 w -v   r f =   1 2 λ f -u v w   (75) with r = (c -1) 2 u 2 + 4(v 2 + w 2 ).
Assuming c = 1, it leads to

λ s = 2u -2r λ a = 2u λ f = 2u + ur r s =   -r v w   r a =   0 w -v   r f =   r v w   (76)
Using definitions introduced in section 2.1, we can observe that the field 2, linked to the eigenvalue λ a is linearly degenerate, while fields 1 and 3, linked to eigenvalues λ s and λ f , are genuinely nonlinear.

To study the entropy of the system, we can easily notice that the system is symmetrizable with the variables:

w =   u v w   (77)
Thanks to Mock theorem (see 1.3.3) the entropic variables are thus the conservative set of variables.

We are now interesting in the determination of the entropy and the entropy fluxes associated. Using previous properties,

w = ∇ w S 0 (w) F (w) = ∇ w S 1 (w) S 1 (w) = F (w) • w -S 1 (w) (78)
we get:

S 0 = 1 2 u 2 + v 2 + w 2 (79)
S * 1 = c u 3 3 + uv 2 + uw 2 (80)
that is,

S 1 = 2u cu 2 3 + v 2 + w 2 (81)
(S 0 , S 1 ) stands for the entropy couple of the system.

Riemann problem solution

In order to study the Riemann problem of the system (73), it is easier to use polar coordinates.

Y =   u r θ  
The system becomes:

Y t + B(Y)Y x Y =   u r θ   B(Y) =   2cu 2r 0 2r 2u 0 0 0 2u   ( 82
)
The eigenvalues of the system are unchanged and the associated eigenvectors are:

r ′ a =   0 0 1   r ′ f,s =   ±1 1 0   (83) t x W L W R W 1 W 2 λ - 1 λ 1 + λ 2 λ - 3 λ + 3 V 1 ( x t ) V 3 ( x t ) -→ n -→ t
∇R i • r i = 0 (84)
In particular an i-Riemann invariant is constant along the corresponding rarefaction wave.

Remarks:

The rarefaction waves can be computed solving the ODE:

V ′ (ξ) = r i (ξ) with ξ = x t (85)
Deriving the Riemann invariant, we get:

d dt r i (w(ξ)) = ∇ w R i • w ′ (ξ) = ∇ w R i • r i (ξ) = 0 (86)
In a linearly degenerate field, the eigenvalue is a Riemann invariant and this property doesn't depend on the set of coordinates we use.

For the field 1, using the associated eigenvectors, we can determine two independent Riemann invariants:

R 1 s = θ R 2 s = u + r (87)
In the same way, for fields 2 and 3, we get:

R 1 a = u R 2 a = r R 1 f = θ R 2 f = u -r (88)
A rarefaction wave corresponds to λ 1 (w L ) < λ 1 (W 1 ), while a shock arises when λ 1 (w L ) > λ 1 (W 1 ). In this case, the solution is discontinue. For shocks we can write the Rankine-Hugoniot relations:

n t [w] + n x [F (w)] = 0, (89)
where n stands for the normal vector, that is :

n t (w 1 -w L ) + n x (F (w 1 ) -F (w L )) = 0 (90)
The tangential and normal vectors are respectively given by:

-→ σ = σ 1 1 and -→ n = 1 -σ 1 (91)
The Rankine Hugoniot relation reads thus:

σ 1 (w 1 -w L ) = F (w 1 ) -F (w L ) (92)
Developing the first term, we get:

σ(u 1 -u L ) = cu 2 1 + v 2 1 + w 2 1 -(cu 2 L + v 2 L + w 2 L ) = c(u 2 1 -u 2 L ) + (v 2 1 -v 2 L ) + (w 2 1 -w 2 L ) = c(u 1 -u L )(u 1 + u L ) + (v 1 -v L )(v 1 + v L ) + (w 1 -w L )(w 1 + w L ) = 2cu(u 1 -u L ) + 2v(v 1 -v L ) + 2w(w 1 -w L ) Setting ã = 1 2 (a 1 -a L )
, the relation becomes:

σ(u 1 -u L ) = 2ũ(u 1 -u L ) + 2u(u 1 -u L ) (93)
And, in the same way, we have:

σ(v 1 -v L ) = 2ṽ(u 1 -u L ) + 2u(v 1 -v L ) (94) σ(w 1 -w L ) = 2 w(u 1 -u L ) + 2u(w 1 -w L ) (95)
Using the fluxes expression, we finally obtain:

σ(w 1 -w L ) = F ′ ( w)(w 1 -w L ) (96)
The Riemann invariants are also shocks invariants. Indeed, concerning the first wave, the Riemann invariant is given by:

R 1 s = θ R 2 s = u + r (97)
Using the two last Rankine-Hugoniot relations 94,95, it becomes: If σ = 2ũ, we finally get:

(σ -2ũ)(v 1 -v L ) = 2ṽ(u 1 -u L ) (98) (σ -2ũ)(w 1 -w L ) = 2 w(u 1 -u L ) (99) (100) w w R w L V v L v R ( v w ) θ
v 1 -v L = ṽ w (w 1 -w L ) (101)
that is:

v 1 -v L w 1 -w L = ṽ w (102)
Because the ratio is constant the 1-Riemann invariant θ is still constant across the shock wave.

To study the other invariants we use polar coordinates. The state

w L =   u L v L w L   (103)
becomes:

w L =   u L r L cos θ L r L sin θ L   ( 104
)
Because θ is constant along rarefaction and across shock waves, we have that:

w 1 =   u 1 r 1 cos θ L r 1 sin θ L   (105)
The averaged state is given thus by:

w =   1 2 (u 1 + u L ) 1 2 (r 1 + r L ) cos θ L 1 2 (r 1 + r L ) sin θ L   (106)
Writing the Rankine-Hugoniot relations in polar coordinates:

σ 1 (u 1 -u L ) = 2u(u 1 -u L ) + 2 r L + r 1 2 cos θ((r 1 -r L ) cos θ L ) + 2 r L + r 1 2 sin θ((r 1 -r L ) sin θ L ) = 2u(u 1 -u L ) + 2r(r 1 -r L ) (107) σ 1 (r 1 -r L ) = 2u(r 1 -r L ) + 2r(u 1 -u L ) (108)
with σ 1 = 2u -2r.

Rewriting 107, we deduce that:

(2u -2r)(u 1 -u L ) = 2u(u 1 -u L ) + 2r(r 1 -r L ) (109)
That is:

u 1 = u L -r 1 + r L (110)
In the same way, the second equation 108 gives a similar result:

u R = u L + r R + r L (111)
The Riemann invariants are thus constant across shocks. We can then deduce the parametrization of the waves curves:

Ms(r) =   uL -r + rL r θL   , Ma(θ) =   uL rL θ   , M f (r) =   uL + r -rL r θL   (112)
Let us write the Lax conditions for shocks 1 and 3. According to the Lax characteristic condition, the speed of a i-shock connecting a left and a right state verifies:

λ i,L > σ > λ i,R (113)
In a 1-shock: σ = 2u L -2r. So

r L < r R and u L > u R (114)
In a 3-shock: σ = 2u L + 2r, i.e.

r L > r R and u L < u R (115)
On the other hand, the Lax entropy condition reads:

σ [S 0 ] ≥ S 1 (116)
Using previous results the inequality becomes:

(2u L -2r 1 ) 1 2 [u 2 + v 2 + w 2 ] -[2u( cu 2 3 + v 2 + w 2 )] ≥ 0 (u L -r 1 )[u 2 + r 2 ] -[2u( cu 2 3 + v 2 + w 2 )] ≥ 0 1. ∀K ∈ T , K is an open set; 2. ∀(K, L) ∈ T × T K ∩ L = ∅; 3. K∈T K = R d .
The thickness of the mesh can be measured by the following parameter

h = sup K∈T |K| |∂K| , (122)
where |K| denoted the volume of the cell K and |∂K| the surface of the cell boundary ∂K.

We are looking for an approximation of the solution w that is polynomial in each cell K. More precisely, let P k (K) be a linear space of polynomials of degree ≤ k defined on the cell K. We denote by P k (K) = (P k (K)) m the corresponding vector space. The approximation is thus discontinuous at the cell boundaries ∂K (and it justifies the name of the method). The approximation space is then

E h = w ∈ L 2 (R d ) m , ∀K ∈ T , w| K ∈ P k (K) (123)
The test functions v are taken in this vector space E h . We multiply the conservation laws by v, integrate on a cell K and sum over all the cells. This leads naturally to the introduction of the following form

B(w, v) = K∈T ∂K f (w L , w R , n)v L - K f i ∂ i v. ( 124
)
The form is linear with respect to v and would be bilinear if the conservation system were linear.

It is necessary to introduce the numerical flux f (w L , w R , n)v L because the solution w and the test function v may be discontinuous at the cell interfaces ∂K.

The approximation consists then in finding an element w in C 1 ([0, T ], E h ) such that for all elements v in E h

x∈R 3 ∂ t w • v + B(w, v) = 0. ( 125
)
The numerical flux has to satisfy

f (w, w, n) = f • n = f i n i (consistence) f (w L , w R , n) = -f (w R , w L , -n) (conservation) (126)
The simplest example is the centered flux

f (w L , w R , n) = 1 2 f i (w L ) + f i (w R ) n i (127)
But this choice leads to oscillations in discontinuous solutions, even if a proper time integration gives a linearly stable scheme. In the next section, we propose other numerical fluxes that lead to better approximation. The approximation can be stated in a more precise way. For this, we consider a basis (e K,i ) of the space P k (K). As a convention, we extend these functions by zero outside K.

We then take v = e K,i in the weak form (125). We obtain that for all the cells L w ′ L,j (t)

L e L,j • e L,i + R∈V (L) ∂L∩∂R f (w L,j e L,j , w R,j e R,j , n L,R ) • e L,i - L f k • ∂ k e L,i = 0. ( 128
)
Is this formula, we have used the Einstein summation convention. We also denote by n L,R the normal vector oriented from cell L to cell R along the boundary ∂L of the cell L (we take as a convention that the Left cell is on the side of -n L,R and the Right cell on the side of n L,R ). The set of the neighboring cells R to the cell L is V (L).

The term L e L,j • e L,i corresponds to a mass matrix term that can be inverted once at the beginning of the computation. If the chosen basis on P k (K) is orthonormal, the mass matrix is diagonal. All the integrals on the cells or their boundaries are computed with a Gauss integration, which we do not describe here in order to avoid heavy notations.

Then, the approximation system is transformed into a first order differential equations system and can be solved by any standard integration algorithm (as Runge-Kutta, Adams, etc.)

Numerical flux

Rusanov flux

A more stable but still very simple numerical flux is the Rusanov flux, which reads λ max = max

0 ξ 1 max 1 j m |λ j (w(w(ξ)))| w(ξ) = ξw L + (1 -ξ)w R f (w L , w R , n) = f (w L ) + f (w R ) 2 • n - λ max 2 (w R -w L ) ( 129
)
where we have noted the entropy variables w and the wave speeds at the state w λ j (w). By taking all the components of v to 1 in the Galerkin weak formulation, we see that the integral of w over the whole space is constant with respect to time, thanks to the conservation property of the flux. It is also possible to state a discrete entropy dissipation property of the scheme. Thanks to the Lax-Wendroff theorem, this property ensures that the scheme converges to entropy solution (when it converges). We would like that

d dt K S 0 + ∂K S(w L , w R , n) 0 ( 130
)
where S(w L , w R , n) is a numerical entropy flux consistant with the entropy flux S i n i . But taking v = w in the Galerkin formulation, we find

d dt K ∇ w S 0 • ∂ t w + ∂K f (w L , w R , n)w L - K ∇ w S i ∂ i w = 0 d dt K ∂ t S 0 + ∂K f (w L , w R , n)w L - K ∂ i S i = 0 d dt K ∂ t S 0 + ∂K f (w L , w R , n)w L -S i (w L )n i = 0 (131)
It is then natural to split the entropy flux into a conservative part and a nonconservative part

f (w L , w R , n)w L -S i (w L )n i = S(w L , w R , n) + D(w L , w R , n) S(w L , w R , n) = f (w L , w R , n) w L + w R 2 - S i (w L ) + S i (w R ) 2 • n i D(w L , w R , n) = - 1 2 f (w L , w R , n)(w R -w L ) -S i (w R )n i -S i (w L )n i (132) Using the fact that f i = ∇ w S i (133)
we have also the following expression for the numerical entropy dissipation

D(w L , w R , n) = - 1 2 ξmax ξmin w ′ (ξ) • (f (w L , w R , n) -f (w(w(ξ))))dξ (134)
where ξ → w(ξ) is an arbitrary parametric curve joining w L and w R

w(ξ min ) = w L w(ξ max ) = w R (135)
A sufficient condition for the scheme to satisfy an entropy condition is thus

1 0 w ′ (ξ) • (f (w L , w R , n) -f (w(w(ξ))))dξ 0 (136)
It can be verified that the Rusanov flux is entropy dissipative (the parametric curve can be here a straight line).

HLLD flux

The HLLD approximate Riemann solver, described in [6], is less dissipative than the Rusanov scheme.

It is based on the same assumption as for HLLC scheme, that is the normal velocity is constant across the Riemann fan. This scheme resolves exactly not only contact discontinuities, as HLLC scheme, but also all isolated discontinuities formed in the MHD system.

Because the HLLD scheme corresponds to the HLLC one when the magnetic field vanishes, all properties of HLLC scheme are preserved and in particular the conservation of positivity.

t x v L v R v * L v * R v * * S L S R S * L S * R
S L = min(u L , u R ) -max(c f L , c f R ) S R = max(u L , u R ) + max(c f L , c f R ) (137)
where c f L and c f R denote the fast magnetoacoustic speeds for left and right states respectively.

S M is computed as the averaged normal velocity from the HLL average.

S M = (S R -u R )ρ R u R -(S L -u L )ρ L u L -p T R + p T L (S R -u R )ρ R -(S L -u L )ρ L (138)
Assuming that the normal velocity and the total pressure are constant across the Riemann fan, we get :

u * L = u * * L = u * * R = u * R = S M p * T L = p * * T L = p * * T R = p * T R = p * T
Using these equalities and (138), the average total pressure is thus given by:

p * T = (S R -u R )ρ R p T L -(S L -u L )ρ L p T R + ρ L ρ R (S R -u R )(S R -u R )(u R -u L ) (S R -u R )ρ R -(S L -u L )ρ L
Then applying the jump condition across S L and S R , the first intermediate states can be evaluated:

ρ * α = ρ α S α -u α S α -S M v * α = v α -B x B yα S M -u α ρ α (S α -u α )(S α -S M ) -B 2 x w * α = w α -B x B zα S M -u α ρ α (S α -u α )(S α -S M ) -B 2 x B * yα = B yα ρ α (S α -u α ) -B 2 x ρ α (S α -u α )(S α -S M ) -B 2 x B * zα = B zα ρ α (S α -u α ) -B 2 x ρ α (S α -u α )(S α -S M ) -B 2 x e * α = e α (S α -u α ) -p T α u α + p * T S M + B x (v α B α -v * α B * α ) S α -S M
We are now considering the inner states. Because of the jump conditions for the continuity equation and for the normal momentum,

ρ * * α = ρ * α p * * T α = p * T α (139)
The appropriate speed of Alfvèn waves are thus given by :

S * L = S M -|Bx| √ ρ * L S * L = S M + |Bx| √ ρ * R (140)
Because of the jump conditions across s * L and S * R waves, the inner states can be determined:

v * * L = v * * R = ρ * L v * L + ρ * R v * R + (B * yR -B * yL )sgn(B x ) ρ * L + ρ * R (141) w * * L = w * * R = ρ * L w * L + ρ * R w * R + (B * zR -B * zL )sgn(B x ) ρ * L + ρ * R (142) B * * yL = B * * yR = ρ * L By * R + ρ * R B * yL + ρ * L ρ * R (v * R -v * L )sgn(B x ) ρ * L + ρ * R (143) B * * zL = B * * zR = ρ * L By * z + ρ * R B * zL + ρ * L ρ * R (w * R -w * L )sgn(B x ) ρ * L + ρ * R (144) e * * α = e * α ∓ ρ * α (v * α B * α -v * * B * * ) sgn(B x ) (145)
Finally the numerical fluxes of the solver are deduced from the integral of conservation laws over the Riemann fan, in the same way as for Rusanov or HLL schemes.

F HLLD =                F L if S L > 0 F * L if S L ≤ 0 ≤ S * L F * * L if S * L ≤ 0 ≤ S M F * * R if S M ≤ 0 ≤ S * R F * R if S * R ≤ 0 ≤ S R F R if S R < 0 (146) t x v L v R v * L v * R u + c ρ u -c ρ u Figure
F * α = F α + S α U * α -S α U α F * * α = F α + S * α U * * α -(S * α -S α )U * α -S L U L (147) 4.2

.3 Multiwave approximate Riemann solver using relaxation

In [START_REF] Bouchut | A multiwave approximate riemann solver for ideal mhd based on relaxation ii -numerical implementation with 3 and 5 wave[END_REF], an approximate Riemann solver for one-dimensional ideal MHD is described derived from a relaxation system. This solver satisfies entropy inequality, preserves the positivity of the density and the internal energy. For more simplicity we present only the 3-wave solver that well resolves fast waves and material contacts.

ρ t + (ρu) x = 0 (ρu) t + (ρu 2 + π) x = 0 (ρu ⊥ ) t + (ρuu ⊥ + π ⊥ ) x = 0 E t + [(E + π)u + π ⊥ • u ⊥ ] x = 0 (B ⊥ ) t + (B ⊥ u -B x u ⊥ ) x = 0 ( 148
)
where the relaxation pressures π and π ⊥ evolve according to

(ρπ) t + (ρπu) x + (|b| 2 + c 2 b )u x -c a b • (u ⊥ ) x = 0 (ρπ ⊥ ) t + (ρπ ⊥ u) x -c a bu x + c 2 a (u ⊥ ) x = 0 (149)
Initially the Riemann problem starts with the relaxation pressures at equilibrium

π = p + 1 2 |B ⊥ | 2 - 1 2 B 2 x π ⊥ = -B x B ⊥
w * α = (ρ α , u * , u * ⊥ , e α , B ⊥α , π * , π * ⊥ ) (150)
where α denotes the l or R and we have:

u * = c L * u L + c R * u R + π L -π R c L + c R u * ⊥ = c L * u ⊥L + c R * u ⊥R + π ⊥L -π ⊥R c L + c R π * = c R * π L + c L * π R -c L c -R(u R -u L ) c L + c R π * ⊥ = c R * π ⊥L + c L * π ⊥R -c L c -R(u ⊥R -u ⊥L ) c L + c R (151)

Time integration

We now address the problem of the time approximation. The Runge-Kutta method is a standard approach that is not described here. We will concentrate on an Adams approach, which has some advantages (the possibility to use very easily different time steps) and some drawbacks (a sometimes more limiting CFL condition).

Adams time integration

In order to obtain the Adams scheme, we suppose that the solution is approximated at some times

t 0 < t 1 < • • • < t n < • • • (152)
by a sequence of elements

w n (•) ≃ w(t n , •) ∈ E h (153)
By an integration in time of (125) we see that the semi-discrete solution satisfies

w L,j (t n+1 ) -w L,j (t n ) L e L,j • e L,i + t n+1 t=t n R∈V (L) ∂L∩∂R f (w L,j e L,j , w R,j e R,j , n L,R ) • e L,i - t=t n+1 t=t n L f k • ∂ k e L,i = 0. ( 154
)
In order to approximate the time integration, we first set

F L/R,i (t) = ∂L∩∂R f (w L,j (t)e L,j , w R,j (t)e R,j , n L,R ) • e L,i S L,i (t) = L f (w L,j (t)e L,j ) k • ∂ k e L,i (155)
in such a way that the weak form can also be written

w L,j (t n+1 ) -w L,j (t n ) L e L,j • e L,i + R∈V (L) t n+1 t=t n F L/R,i (t) - t=t n+1 t=t n S L,i (t) = 0. ( 156
)
We also define the discrete flux and source term at time t n

F n L/R,i = ∂L∩∂R f (w n L,j e L,j , w n R,j e R,j , n L,R ) • e L,i S n L,i = L f (w n L,j e L,j ) k • ∂ k e L,i (157)
We then construct the time interpolation polynomials F L/R,i (t) and S L,i (t) of F n L/R,i and S n L,i by using r + 1 interpolation points. In other word, the interpolation polynomials are defined by

F L/R,i (t n-l ) = F n-l L/R,i and S L,i (t n-l ) = S n-l L,i , l = 0 • • • r (158)
The time integration of the boundary terms and the source terms are then obtained by an exact integration of F L/R,i (t) and S L,i (t) on the interval [t n , t n+1 ]

t n+1 t=t n ∂L∩∂R f (w L,j e L,j , w R,j e R,j , n L,R ) • e L,i ≃ t n+1 t=t n F L/R,i (t)dt, t=t n+1 t=t n L f k • ∂ k e L,i ≃ t=t n+1 t=t n S L,i (t)dt. ( 159
)
The method requires to store the flux terms on the cell edges and the source terms in the cell at the r + 1 previous times. It is also necessary to initialize the scheme, for example by r steps of a Runge-Kutta algorithm.

Multi time steps approach

An advantage of the Adams approach is that it is quite easy to adapt it to the case where the time step are different from one cell to one another. This is useful when small cells are mixed with big cells in order to reduce the computational cost. The first step is to attribute to each cell what we call a CFL level, which is only based on a geometric criterion. We define

h K = |K| |∂K| h min = min K∈T h K h max = max K∈T h K (160)
With a standard one-step time method, the time step is fixed by the smallest

cell ∆t = CF L × h min λ max (161)
where λ max is the highest wave speed in the mesh. We shall say that a cell K is of level n (and we note level

(K) = n) if 2 n h min h K < 2 n+1 h min (162)
In this way, the smallest cells are of level n = 0 and the biggest cells are of level A second order scheme is applied with CF L = 0.3 and t f inal = 0.5s..

N = log 2 h max h min (163)
While the computation elapses in 112.81 s. usually, it lasts 34.29s. using the multi time step approach.

Theoretical stability study

Slope limiter

Unfortunately, despite its clean construction, and like other numerical methods for discontinuous solutions, the Discontinuous Galerkin approximation may suffer from oscillations in the shocks and contacts. It is thus necessary to add to the method a step of slope limiting before the computation of the fluxes. Several approaches exist. In this work, we use a very simple method that is entropy dissipative. It is very satisfying for one-dimensional problems. However, it is sometimes too robust and dissipative for higher space dimensions.

The principle is classical: we have to avoid the apparition of a local extremum. For this, we consider in each cell L a vector of limiters which is set initially to

β i L = 1, i = 1 • • • m (167)
(no limitation). The maximum limitation is achieved when all the components are equal to zero. After the slope limiting step, each component i of the conservative variables vector will be changed to the following convex combination between the non-limited values and the mean value on the cell of the approximation

w i L,j (t) = β i L w i L,j (t) + (1 -β i L )w i L (168)
where w i L denotes the mean value of the ith component of w L . We observe that when all the limiter components β i L = 0, we recover the mean value of the approximation and the scheme degenerates to first order. It is also easy to see that the slope limiting step is necessarily conservative. Finally, thanks to the Jensen's inequality, it is also entropy dissipative. Now, in order to compute β, we use the following algorithm. First, we loop on all the cell edges L/R. We then loop on all the Gauss points G of the edge and compute the value of w L (G) and w R (G). If for a component i, w i L (G) (respectively w i R (G)) is not between w i L and w i R then decrease the limiter β i R (respectively β i L ) accordingly. More precisely, we take as new limiter value

β i L ← max(min(β i L , w i R -w i L w i (G) -w i L ), 0) (169)
(and the symmetric formula for β i R ). At the end of this algorithm, the values of w L in the cell L are replace by the values given in (168).

6 Numerical results

Simplified MHD

The system is solved using a Galerkin Discontinuous approximation technique (see section 4) and a standard Godunov scheme is employed for fluxes approximation. The following graphs Fig. 6 and Fig. 7 present the results obtained for the initial condition above:

W L =   1 1 0   , W R =   1/2 1 0   . (170)
As we can see for MHD system, non uniqueness of the Riemann problem solution can occur. The initial condition above leads to such phenomenon:

W L =   1 1 0   , W R =   1/2 -1 0   . (171)
The left state can be linked to the right state by an entropy satisfying shock wave. Indeed for a shock speed σ = 3/2 Rankine Hugoniot conditions hold. In

σ [S 0 ] -[S 1 ] = 49 48 > 0, (172)
The shock satisfies entropy condition. But using previous results and classical techniques we can also point out an other entropy satisfying shock. The Riemann problem allows thus several solutions.

As shown on Fig. 8 and Fig. 9 he numerical results obtained for this initial state exemplify the property.

For MHD non uniqueness arises when the initial conditions are coplanar, that is to say when the transverse magnetic field has a different orientation from right to left. For the simplified system the phenomenon is the same and the numerical scheme seems to converge toward an entropy satisfying solution different from the chosen exact solution obtained from the theory given in section 3. We give below a measure of the L 1 norm of the error, evaluated on a test case of a sinusoidal function convected on a domain with periodic boundaries.

The plotting of L 1 error, see Fig. 20, corresponds to density rho. The straight line corresponds to first order, the dashdot line to the second order and finally the dashed line refers to third order.

2D academic test cases

Convergence test case

As for 1D test case, we give a measure of the L 1 norm of the error. The test case employed is again a sinusoidal function, of period 2, convected on the whole domain. The plot 21 presents the rate of convergence for first and second order.

Figure 1 :

 1 Figure 1: Schematic structure of the Riemann fan with four intermediate states

Figure 2 :

 2 Figure 2: Representation of the ratio ṽ w

Figure 3 :

 3 Figure 3: Schematic structure of the Riemann fan with four intermediate states

4 :

 4 Figure 4: Schematic structure of the Riemann fan with two intermediate states

  The parameters c a , c b and b read √ ρ|B x |,ρ √ p ′ and sgn(B x ) √ ρB ⊥ respectively. Assuming b = 0 and c a = c b = c, the solver is reduced to 3 waves and 2 intermediate states. The 2 intermediate states read:

Figure 5 :

 5 Figure5: 2D mesh and its 7 levels of CFL and a zoom of the tiny cells region

Figure 6 :

 6 Figure 6: Results of the planar test case for the simplified MHD : u

Figure 9 :Figure 10 :Figure 11 :Figure 12 :Figure 13 :Figure 14 :Figure 15 :Figure 16 :Figure 17 :Figure 18 :

 9101112131415161718 Figure 9: Results of the coplanar test case for the simplified MHD : v

Figure 19 :Figure 20 :Figure 21 :

 192021 Figure 19: All seven waves : B y at time t = 0.1 with a resolution ∆x = 0.001, with slope limiter

This work has been initiated during a CEMRACS project in July and August 2008 in Luminy. The project was entitled GADMHD (It has been supported by the INRIA CALVI project

Thanks to relation (110), we compute :

Finally thanks to (S 0 , S 1 ) definition and relation (110), the inequality becomes:

The characteristic condition u L > u R is recovered and the Lax characteristic and entropy conditions are equivalent in this framework. The same inequality is obtained for the 3-wave with a change of sign.

To solve the Riemann problem we still have to determine the functions V 1 ( x t ) and V 3 ( x t ) in the 1 and 3 rarefaction waves.

In the 1-rarefaction, the Riemann invariants are constants i.e.:

Using this property and solving the EDO (85), we get a set of equations:

The same kind of parametrization is obtained for the 3-wave:

4 Discontinuous Galerkin approximation

The Discontinuous Galerkin (DG) approximation technique is a generalization of the finite volume approach in order to achieve higher order. It is well suited to hyperbolic systems of conservation laws. Its application to the MHD equations is studied for example by Barth in [START_REF] Barth | On the role of involutions in the discontinous Galerkin discretization of Maxwell and magnetohydrodynamic systems[END_REF]. In order to make the presentation simpler, we first present the space semi-discrete version of the scheme. The time integration will be studied later on.

Space approximation

We are interested in an approximation of the following system

stated in the whole space R d (the boundary conditions problematic is not addressed in this document). Let us consider a mesh T of R d made of cells K satisfying

We also define a level for the edge

Let ∆t be the time step associated to the biggest cells (we call it the macro time step)

According to the previous definitions, this time step satisfies

The time loop algorithm is then the following for i =1 to 2 N do let j be the biggest integer such that 2 j divides i for all the edges L/R of level ≤ j do compute the integral of the flux term F L/R on a time interval of length ∆t j-N and distribute it to the two neighboring cells end for for all the cells L of level ≤ j do compute the integral of the source term S L on a time interval of length ∆t j-N and distribute it to the corresponding cell. end for Update only the cells of level ≤ j end for With this algorithm, the fluxes are computed more times on the small edges than on the big edges but are always distributed on the two sides of the edge in order to keep a conservative scheme. The time integration is always performed by the Adams approach: the interpolation polynomial is calculated from the r more recent fluxes or sources evaluations. At the end of a macro time step of size ∆t (when i = 2 N in the algorithm) all the cells are updated together.

If the number of small cells is small and if the other cells have almost the same size, the gain is almost of 2 N .

To illustrate this assertion we give some examples for one and two dimension problems.

1D test case:

We consider the transport of a sinusoidal function through the computational domain. A transformation is applied to the mesh which consists in map x to x 2 . We obtain thus 8 levels of CFL. We use a second order scheme with CF L = 0.2 and t f inal = 1s.

Without the multi time stepping approach, the computation lasts 51 s. while it takes 2.9 s. using the algorithm.

2D test case:

The second test case is a simple Sod test applied to compressible Euler equations. The domain is discretized as follow in Fig. 5:

One-dimensional test cases

The MHD equations possess seven eigenvalues, some of which may coincide depending on the direction and the strength of the magnetic field. As pointed out in [START_REF] Brio | An upwind differencing scheme for the equations of ideal magnetohydrodynamics[END_REF] the MHD system is thus non-strictly hyperbolic and possess non-convexity. As a consequence a solution of the Riemann problem may be composed not only of ordinary shock and rarefaction waves but also other waves as compound waves.

Thanks to the solver presented in Section 2 or with Torrilhon's Riemann solver, available at http://wwwmath.ethz.ch/~matorril/mhdsolver, exact solutions can be computed in most cases.

For all one-dimensional test cases γ = 5

3 , b = B x is constant and x 0 = 0 denotes the position where the initial discontinuity is applied.

Compound shocks

We propose to study first a compound shocks test case. As shown in Fig. 10, Fig. 11 and Fig. 12, this problem shows the formation of a left-going slow compound wave with a weak right-going slow shock and a contact discontinuity.

Left state

Right state ρ

w R = 0. p L = 1. p = 0.1 B x,L = 0.75 B x = 0.75 B y,L = 1.

B y = -1. B z,L = 0.

B z = 0.

All seven waves

This problem shows the formation of all seven possible MHD waves. All magnetosonic shock are weak shocks. See