François Bouchut
email: francois.bouchut@u-pem.fr

Xavier Lhébrard
email: xavier.lhebrard@u-pem.fr

A 5-wave relaxation solver for the shallow water MHD system

Keywords: Shallow water magnetohydrodynamics, approximate Riemann solver, relaxation, contact discontinuities, entropy inequality Mathematics Subject Classification: 76W05, 76M12, 35L65

come

A 5-wave relaxation solver for the shallow water MHD system

Introduction

The shallow water magnetohydrodynamic (SWMHD) system has been introduced in [START_REF] Gilman | Magnetohydrodynamic "shallow water" equations for the solar tachocline[END_REF] to describe the thin layer evolution of the solar tachocline. It is written in 2d in the tangent plane approximation as

∂ t h + ∇ • (hu) = 0, (1.1)
∂ t (hu) + ∇ • (hu ⊗ u -hb ⊗ b) + ∇(gh 2 /2) = 0, (1.2)
∂ t (hb) + ∇ • (hb ⊗ u -hu ⊗ b) + u∇ • (hb) = 0, (1.3)
where g > 0 is the gravity constant, h ≥ 0 is the thickness of the fluid, u = (u,v) is the velocity, b = (a,b) is the magnetic field, and the notation ∇ • (b ⊗ u) is for the vector with index i given by j ∂ j (b i u j ). The system should be complemented with Coriolis force and topography, but these sources will not be considered in this paper.

The system (1.1)-(1.3) is endowed with an entropy (energy) inequality

∂ t 1 2 h|u| 2 + 1 2 gh 2 + 1 2 h|b| 2 +∇ • 1 2 h|u| 2 + gh 2 + 1 2 h|b| 2 u -hb(b • u) ≤ 0, (1.4)
that becomes an equality in the absence of shocks.

The induction equation (1.3) implies, by taking its divergence, that

∂ t ∇ • (hb) + ∇ • u∇ • (hb) = 0, (1.5)
meaning that ∇ • (hb) is transported at velocity u. In particular, ∇ • (hb) remains identically zero if it vanishes initially. This situation ∇ • (hb) = 0 (that cancels the last term in (1.3)) is indeed the physically relevant one, but for numerical purposes it is convenient to relax this constraint and consider general data. The particular form (1.3) has been introduced in [START_REF] Dellar | Hamiltonian and symmetric hyperbolic structures of shallow water magnetohydrodynamics[END_REF] for the SWMHD, and in [START_REF] Janhunen | A Positive conservative method for magnetohydrodynamics based on HLL and Roe methods[END_REF] for the full MHD system. It enables to use one-dimensional solvers in two dimensions, indeed this is why the term u∇•(hb) has been added in (1.3).

Multidimensional simulations of the SWMHD system have been performed in [START_REF] Rossmanith | A constrained transport method for the shallow water MHD equations[END_REF][START_REF] Rossmanith | A wave propagation algorithm for hyperbolic systems on curved manifolds[END_REF]. The system is closely related to the MHD system, to which many works have been devoted. An important issue in multidimensional simulations is to minimize the numerical viscosity by using accurate solvers, in particular on contact discontinuities; while being robust, see for example [START_REF] Fuchs | Approximate Riemann solvers and robust high-order finite volume schemes for multi-dimensional ideal MHD equations[END_REF][START_REF] Waagan | A robust numerical scheme for highly compressible magnetohydrodynamics: Nonlinear stability, implementation and tests[END_REF][START_REF] Balsara | Multidimensional HLL and HLLC Riemann solvers for unstructured meshes with application to Euler and MHD flows[END_REF].

If dependency is only in one spatial variable x, the system simplifies to

∂ t h + ∂ x (hu) = 0, (1.6)
∂ t (hu) + ∂ x (hu 2 + P ) = 0, (1.7) ∂ t (hv) + ∂ x (huv + P ⊥ ) = 0, (

∂ t (ha) + u∂ x (ha) = 0, (1.9)

∂ t (hb) + ∂ x (hbu -hav) + v∂ x (ha) = 0, (1.10)
with

P = g h 2 2
ha 2 , P ⊥ = -hab, (1.11) and the energy inequality (1.4) becomes

∂ t 1 2 h(u 2 + v 2 ) + 1 2 gh 2 + 1 2 h(a 2 + b 2 ) +∂ x 1 2 h(u 2 + v 2 ) + gh 2 + 1 2 h(a 2 + b 2 ) u -ha(au + bv) ≤ 0.
(1.12)

The eigenvalues of the system (1.6)-(1.10) are u, u ± |a|, u ± a 2 + gh. The associated waves are called respectively material (or divergence) waves, Alfven waves and magnetogravity waves, see [START_REF] Desterck | Hyperbolic theory of the shallow water magnetohydrodynamics equations[END_REF][START_REF] Zeitlin | Remarks on rotating shallow-water magnetohydrodynamics[END_REF]. Some of these waves will have the same speed when a or h vanishes, hence the system is nonstrictly hyperbolic.

The system has three types of contact discontinuities corresponding to linearly degenerate eigenvalues: the material contacts associated to the eigenvalue u, the left Alfven contacts associated to u -|a|, and the right Alfven contacts associated to u + |a|. The jump relations associated to these contact discontinuities are as follows. Across a material contact, the quantities u, v, g h 2 2ha 2 , hab are constant. Across an Alfven contact, the quantities h, u, a are constant, and moreover for a left Alfven contact b sgn av is constant, while for a right Alfven contact b sgn a + v is constant.

The system (1.6)-(1.10) is nonconservative in the variables ha, hb. However, ha jumps only through the material contacts, where u and v are continuous. Therefore, there is indeed no ambiguity in the non conservative products u∂ x (ha) and v∂ x (ha), that are well-defined.

A finite volume scheme for the quasilinear system (1.6)-(1.10) can be classically built following Godunov's approach, considering piecewise constant approximations of U = (h,hu,hv,ha,hb), (1.13) and invoking an approximate Riemann solver at the interface between two cells, see for example [START_REF] Harten | On upstream differencing and Godunovtype schemes for hyperbolic conservation laws[END_REF] or [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF]Section 2.3]. A difficulty is however that the system is not conservative. The energy is nevertheless obviously convex with respect to U . The SWMHD system is closely related to the compressible MHD system, for which several entropy schemes are known [START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF][START_REF] Bézard | An entropic solver for ideal Lagrangian magnetohydrodynamics[END_REF][START_REF] Després | Lagrangian systems of conservation laws. Invariance properties of Lagrangian systems of conservation laws, approximate Riemann solvers and the entropy condition[END_REF][START_REF] Bouchut | A multiwave approximate Riemann solver for ideal MHD based on relaxation II -Numerical implementation with 3 and 5 waves[END_REF]. In this paper we apply the relaxation approach of [START_REF] Bouchut | Entropy satisfying flux vector splittings and kinetic BGK models[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF][START_REF] Bouchut | A multiwave approximate Riemann solver for ideal MHD based on relaxation I -Theoretical framework[END_REF][START_REF] Coquel | Relaxation of fluid systems[END_REF] to the SWMHD system, in order to get an approximate Riemann solver that is entropy satisfying, ensuring robustness, while being exact on isolated Alfven contacts. The relaxation system is of Suliciu type as introduced in [START_REF] Suliciu | Some stability-instability problems in phase transitions modelled by piecewise linear elastic or viscoelastic constitutive equations[END_REF], and the approximate Riemann solver belongs to the family of HLLC solvers, as in [START_REF] Harten | On upstream differencing and Godunovtype schemes for hyperbolic conservation laws[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF][START_REF] Berthon | An HLLC scheme to solve the M1 model of radiative transfer in two space dimensions[END_REF][START_REF] Sangam | An HLLC scheme for ten-moments approximation coupled with magnetic field[END_REF][START_REF] Berthon | A local entropy minimum principle for deriving entropy preserving schemes[END_REF][START_REF] Fuchs | Approximate Riemann solvers and robust high-order finite volume schemes for multi-dimensional ideal MHD equations[END_REF][START_REF] Balsara | Multidimensional HLL and HLLC Riemann solvers for unstructured meshes with application to Euler and MHD flows[END_REF]. The paper is organized as follows. In Section 2 we describe the relaxation approximate solver and its entropy property. In Section 3 we derive explicit optimal choices of the speeds that enable to obtain stability and accuracy. In Section 4 we state our main theorem giving the properties of our relaxation approximate Riemann solver. Finally, in Section 5 we perform numerical tests.

Approximate Riemann solver 2.1 Relaxation approach

In order to get an approximate Riemann solver for (1.6)-(1.10), we use a standard relaxation approach, used for example in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF] for the Euler equations, in [START_REF] Bouchut | A multiwave approximate Riemann solver for ideal MHD based on relaxation II -Numerical implementation with 3 and 5 waves[END_REF] for the MHD equations and in [START_REF] Bouchut | A new model for shallow viscoelastic fluids[END_REF] for shallow elastic fluids. An abstract general description can be found in [START_REF] Coquel | Relaxation of fluid systems[END_REF], and related works are [START_REF] Chalons | Well-balanced time implicit formulation of relaxation schemes for the Euler equations[END_REF][START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF]. The approach enables to naturally handle the energy inequality (1.12), and also preserves the positivity of density. Its structure has however to be well-chosen in order to resolve exactly isolated Alfven contacts.

We introduce new variables π,π ⊥ , the relaxed pressures, and c a , c intended to parametrize the speeds. The form of the relaxation system is as follows,

∂ t h + ∂ x (hu) = 0, (2.1)
∂ t (hu) + ∂ x (hu 2 + π) = 0, (2.2) ∂ t (hv) + ∂ x (huv + π ⊥ ) = 0, (2.3)
∂ t (ha) + u∂ x (ha) = 0, (2.4) ∂ t (hb) + ∂ x (hbu -hav) + v∂ x (ha) = 0, (2.5)
∂ t (hπ) + ∂ x (hπu) + c 2 ∂ x u = 0, (2.6)
∂ t (hπ ⊥ ) + ∂ x (hπ ⊥ u) + c 2 a ∂ x v = 0, (2.7)
∂ t c + u∂ x c = 0, (2.8) ∂ t c a + u∂ x c a = 0.
(2.9)

The approximate Riemann solver can be defined as follows, starting from left and right values U l , U r at an interface.

• Solve the Riemann problem for (2.1)-(2.9) with initial data obtained by completing U l , U r by the equilibrium relations

π l = P l ≡ (gh 2 /2 -ha 2 ) l , π r = P r ≡ (gh 2 /2 -ha 2 ) r , (π ⊥ ) l = (P ⊥ ) l ≡ (-hab) l , (π ⊥ ) r = (P ⊥ ) r ≡ (-hab) r , (2.10)
and with suitable positive values of c l , c r , c a,l , c a,r that will be discussed further on, essentially in Section 3.

• Retain in the solution only the variables h, hu, hv, ha, hb. The result is a vector called R(x/t, U l ,U r ).

We can remark that the relaxation system (2.1)-(2.9) is identical to the 5-wave relaxation system in [9, equations (5.5)-(5.7) with b = 0], with the identification B x = ha, B ⊥ = hb. However, the initialization (2.10) differs from that of the MHD equations, [9, equation (2.10)]. Indeed, the homogeneity in magnetic terms is different in the SWMHD and MHD systems. Intuitively, the solver is consistent because of the equations (2.1)-(2.5), that are consistent with (1.6)- (1.10). The specific values used for c, c a do not play any role in this consistency. However, if we require the solver to have the highest accuracy, i.e. to be "tangent" to the original system, one has to take the speeds c > c a > 0 as approximations of h a 2 + gh and h|a| respectively, in the limit when U l , U r are close to a common value U . This is because, as can be checked with straightforward computations, smooth solutions to (1.6)-(1.10) verify

∂ t (hP ) + ∂ x (hP u) + h 2 (a 2 + gh)∂ x u = 0, ∂ t (hP ⊥ ) + ∂ x (hP ⊥ u) + h 2 a 2 ∂ x v = 0, (2.11)
that have to be compared with (2.6), (2.7). The accuracy of the solver on isolated contacts is described by the following lemma.

Lemma 2.1. The approximate Riemann solver R(x/t, U l ,U r ) solves exactly:

(i) material contact discontinuities, (ii) left Alfven contact discontinuities under the condition c a,l = (h|a|) l , (iii) right Alfven contact discontinuities under the condition c a,r = (h|a|) r .

Proof. Material contacts are solutions to the SWMHD system (1.6)-(1.10) with u, v, P , P ⊥ constant. These solutions are obviously solutions to the relaxation system (2.1)-(2.9) with π = P , π ⊥ = P ⊥ . Thus for these data, R coincides with the exact solver, which proves (i).

Alfven contacts are solutions to (1.6)-(1.10) with h, u, a, b sgn av (for a left contact), b sgn a + v (for a right contact) constant. As previously, it is enough to prove that these solutions, completed with π = P , π ⊥ = P ⊥ , are solutions to the relaxation system (2.1)-(2.9). One can see that only (2.7) is not immediately satisfied. Comparing (2.7) to the second line of (2.11), we get the condition c 2 a = h 2 a 2 where v jumps. Note that according to (2.4) and (2.9), ha and c a are both continuous through the Alfven waves (assuming a = 0). This yields (ii) and (iii).

Following the Godunov approach, the numerical scheme can be defined by the approximate Riemann solver as follows. We consider a mesh of cells

(x i-1/2 , x i+1/2 ), i ∈ Z, of length ∆x i = x i+1/2 -x i-1/2 , discrete times t n with t n+1 -t n = ∆t,
and cell values U n i approximating the average of U over the cell i at time t n . We can then define an approximate solution U appr (t,x) for t n ≤ t < t n+1 and x ∈ R by

U appr (t,x) = R( x -x i+1/2 t -t n , U n i , U n i+1 ) for x i < x < x i+1 , (2.12)
where x i = (x i-1/2 + x i+1/2 )/2. This definition is coherent under a half CFL condition, formulated as

x/t < - ∆x i 2∆t ⇒ R(x/t, U i , U i+1 ) = U i , x/t > ∆x i+1 2∆t ⇒ R(x/t, U i , U i+1 ) = U i+1 .
(2.13)

The new values at time t n+1 are defined by

U n+1 i = 1 ∆x i x i+1/2 x i-1/2 U appr (t n+1 -0, x)dx. (2.14)
Notice that it is only in this averaging procedure that the choice of the particular pseudo-conservative variable U as (1.13) is involved. We can follow the computations of [6, Section 2.3], the only difference being that the system is not conservative. We obtain the update formula

U n+1 i = U n i - ∆t ∆x i (F l (U n i , U n i+1 ) -F r (U n i-1 , U n i )), (2.15)
where

F l (U l ,U r ) = F (U l ) - 0 -∞ (R(ξ, U l , U r ) -U l )d ξ, F r (U l ,U r ) = F (U r ) + ∞ 0 (R(ξ, U l , U r ) -U r ) dξ, (2.16)
the variable ξ stands for x/t, and the pseudo-conservative flux is chosen as

F (U ) ≡ (hu, hu 2 + P, huv + P ⊥ , 0, hbu -hav).
(2.17)

In (2.17), the two last components could be chosen differently since the two magnetic equations in our system are not conservative. We can remark that the choice of F has no influence on the update formula (2.15).

Energy inequality

Here we do not use the entropy extension defined in [6, Definition 2.14], because the minimization principle is a bit too restrictive. We instead follow the strategy used in [START_REF] Bouchut | A multiwave approximate Riemann solver for ideal MHD based on relaxation I -Theoretical framework[END_REF][START_REF] Bouchut | A new model for shallow viscoelastic fluids[END_REF]. We define the left and right numerical energy fluxes as

G l (U l ,U r ) = G(U l ) - 0 -∞ E R(ξ, U l , U r ) -E(U l ) dξ, G r (U l ,U r ) = G(U r ) + ∞ 0 E R(ξ, U l , U r ) -E(U r ) dξ, (2.18)
where E and G are respectively the energy and the energy flux from (1.12),

E(U ) = 1 2 h(u 2 + v 2 ) + 1 2 gh 2 + 1 2 h(a 2 + b 2 ), G(U ) = 1 2 h(u 2 + v 2 ) + gh 2 + 1 2 h(a 2 + b 2 ) u -ha(au + bv). (2.19)
Following [6, Section 2.3], a sufficient condition for the scheme to be energy satisfying is that

G r (U l , U r ) -G l (U l ,U r ) ≤ 0. (2.20)
When this is satisfied, because of the convexity of E with respect to U and of the CFL condition (2.13), one has the discrete energy inequality

E(U n+1 i ) -E(U n i ) + ∆t ∆x i G(U n i ,U n i+1 ) -G(U n i-1 ,U n i ) ≤ 0, (2.21)
where the numerical energy flux

G(U l , U r ) is any function satisfying G r (U l ,U r ) ≤ G(U l , U r ) ≤ G l (U l , U r ).
In order to analyze the condition (2.20), let us introduce the sum of the gravitational potential energy and the magnetic energy

e = 1 2 gh + 1 2 (a 2 + b 2 ), (2.22)
that enables to rewrite the energy as

E = 1 2 h(u 2 + v 2 ) + he. (2.23)
Then, while solving the relaxation system (2.1)-(2.9), we solve simultaneously the equation for a new variable ê,

∂ t (ê -π 2 /2c 2 -π 2 ⊥ /2c 2 a ) + u∂ x (ê -π 2 /2c 2 -π 2 ⊥ /2c 2 a ) = 0, (2.24)
where ê has left and right initial data e l = e(U l ), e r = e(U r ). The reason for writing (2.24) is that combined with (2.1)-(2.9) it implies

∂ t 1 2 h(u 2 + v 2 ) + hê + ∂ x 1 2 h(u 2 + v 2 ) + hê u + πu + π ⊥ v = 0. (2.25)
Indeed, (2.25) can be obtained as follows. From (2.2), (2.3), (2.6), (2.7) (combined with (2.1)), we get

∂ t u + u∂ x u + 1 h ∂ x π = 0, ∂ t v + u∂ x v + 1 h ∂ x π ⊥ = 0, ∂ t π + u∂ x π + c 2 h ∂ x u = 0, ∂ t π ⊥ + u∂ x π ⊥ + c 2 a h ∂ x v = 0.
(2.26)

Multiplying these equations respectively by u, v, π, π ⊥ , we obtain

∂ t u 2 /2 + u∂ x u 2 /2 + u h ∂ x π = 0, ∂ t v 2 /2 + u∂ x v 2 /2 + v h ∂ x π ⊥ = 0, ∂ t π 2 /2 + u∂ x π 2 /2 + πc 2 h ∂ x u = 0, ∂ t π 2 ⊥ /2 + u∂ x π 2 ⊥ /2 + π ⊥ c 2 a h ∂ x v = 0.
(2.27)

Using the advection equations (2.8), (2.9) for c and c a , the two last equations of (2.27) give

∂ t π 2 2c 2 + u∂ x π 2 2c 2 + π h ∂ x u = 0, ∂ t π 2 ⊥ 2c 2 a + u∂ x π 2 ⊥ 2c 2 a + π ⊥ h ∂ x v = 0.
(2.28)

Adding up the two first equations of (2.27) and (2.28) yields

∂ t u 2 + v 2 2 + π 2 2c 2 + π 2 ⊥ 2c 2 a + u∂ x u 2 + v 2 2 + π 2 2c 2 + π 2 ⊥ 2c 2 a + 1 h ∂ x (πu + π ⊥ v) = 0.
(2.29) Then, adding this to (2.24) and using (2.1) finally gives our stated identity (2.25).

Using the value of the Riemann solution to the relaxation system at x/t = 0, we define where here U = R(x/t, U l , U r ), e(U ) is defined in (2.22) and ê is defined by

G(U l ,U r ) = 1 2 h(u 2 + v 2 ) + hê u + πu + π ⊥ v x/t=0 . ( 2
(2.24), then G r (U l ,U r ) ≤ G(U l , U r ) ≤ G l (U l , U r ),
G(U l ,U r ) = G(U l ) - 0 -∞ 1 2 h(u 2 + v 2 ) + hê (ξ) -E(U l ) dξ = G(U r ) + ∞ 0 1 2 h(u 2 + v 2 ) + hê (ξ) -E(U r ) dξ.
(2.32) Therefore, comparing to (2.18), we see that in order to get

G r ≤ G ≤ G l it is enough that for all ξ E(R(ξ, U l ,U r )) ≤ 1 2 h(u 2 + v 2 ) + hê (ξ), (2.33)
which is equivalent to (2.31).

Intermediate states

In this subsection we describe the solution to the Riemann problem for (2.1)-(2.9) with initial data completed by the relations (2.10). The analysis is similar to that in [START_REF] Bouchut | A multiwave approximate Riemann solver for ideal MHD based on relaxation I -Theoretical framework[END_REF][START_REF] Bouchut | A multiwave approximate Riemann solver for ideal MHD based on relaxation II -Numerical implementation with 3 and 5 waves[END_REF] for the full MHD system, and to [START_REF] Bouchut | A new model for shallow viscoelastic fluids[END_REF] for shallow elastic fluids.

The quasilinear system (2.1)-(2.9) has the property of having a quasi diagonal form

∂ t (π + cu) + (u + c/h)∂ x (π + cu) - u h c∂ x c = 0, (2.34)
∂ t (π -cu) + (u -c/h)∂ x (π -cu) - u h c∂ x c = 0, (2.35) ∂ t (1/h + π/c 2 ) + u∂ x (1/h + π/c 2 ) = 0, (2.36) ∂ t (π ⊥ + c a v) + (u + c a /h)∂ x (π ⊥ + c a v) - v h c a ∂ x c a = 0, (2.37) ∂ t (π ⊥ -c a v) + (u -c a /h)∂ x (π ⊥ -c a v) - v h c a ∂ x c a = 0, (2.38) ∂ t (b + ha c 2 a π ⊥ ) + u∂ x (b + ha c 2 a π ⊥ ) = 0, (2.39)
∂ t (ha) + u∂ x (ha) = 0, (2.40) ∂ t c + u∂ x c = 0, (2.41)
∂ t c a + u∂ x c a = 0. (2.42)
One deduces its eigenvalues, which are

u - c h , u - c a h , u, u + c a h , u + c h , (2.43)
the central eigenvalue u having multiplicity 5 and the other being simple. From the above form one also checks easily that the system is hyperbolic, with all eigenvalues linearly degenerate. As a consequence, Rankine-Hugoniot conditions are well defined (the weak Riemann invariants do not jump through the associated discontinuity), and are equivalent to any conservative formulation. In the solution to the Riemann problem, the speeds corresponding to the previous eigenvalues will be denoted by

Σ 1 < Σ 2 < Σ 3 < Σ 4 < Σ 5 , (2.44)
the speed Σ 3 corresponding to the eigenvalue u. Thus we get a 5-wave solver with four intermediate states. The variables take the values "l" for x/t < Σ 1 , "l**" for Σ 1 < x/t < Σ 2 , "l*" for Σ 2 < x/t < Σ 3 , "r*" for Σ 3 < x/t < Σ 4 , "r**" for Σ 4 < x/t < Σ 5 , "r" for Σ 5 < x/t, see Figure 1. There are 5 strong Riemann

x U r x/t = Σ 5 U * * r x/t = Σ 4 U * l U * r x/t = Σ 3 U * * l x/t = Σ 2 x/t = Σ 1 U l 0
u -c/h : π + cu, π ⊥ , v, u -c a /h : π ⊥ + c a v, π, u, u + c a /h : π ⊥ -c a v, π, u, u + c/h : π -cu, π ⊥ , v.
(2.46)

We notice that the equations (2.1), (2.2), (2.6), (2.8) form a closed system of equations. Therefore, the variables h, u, π, c can be resolved independently of the knowledge of v, π ⊥ , a, b, c a , and in particular they do not jump through the Σ 2 and Σ 4 waves. This means that for these unknowns h, u, π, c, the "l**" and "l*" states are identical, as well as the "r**" and "r*" states. In particular, u and π take the constant values u * and π * on the whole fan

Σ 1 < x/t < Σ 5 , (u,π) * * l = (u,π) * l = (u,π) * , (u,π) * * r = (u,π) * r = (u,π) * .
(2.47)

The second velocity-pressure set of variables v, π ⊥ jump only though the c a waves Σ 2 , Σ 4 , thus

(v,π ⊥ ) * * l = (v,π ⊥ ) l , (v,π ⊥ ) * * r = (v,π ⊥ ) r , (v,π ⊥ ) * l = (v,π ⊥ ) * r = (v,π ⊥ ) * .
(2.48)

Then, because of (2.45), the variables a, c a do not jump through the Σ 2 and Σ 4 waves, as h,u,π,c. Moreover, b does not jump through the Σ 1 and Σ 5 waves, as v,π ⊥ , see Figure 2. In addition, computations using the Riemann invariants

x v r ,π ⊥r ,b r x/t = Σ 5 v r ,π ⊥r ,b r x/t = Σ 4 v * ,π * ⊥ ,b * l v * ,π * ⊥ ,b * r x/t = Σ 3 v l ,π ⊥l ,b l x/t = Σ 2 x/t = Σ 1 v l ,
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 , v * = c al v l + c ar v r + π ⊥l -π ⊥r c al + c ar , π * ⊥ = c ar π ⊥l + c al π ⊥r -c al c ar (v r -v l ) c al + c ar , (2.49) and 1 h * l = 1 h l + c r (u r -u l ) + π l -π r c l (c l + c r ) , 1 h * r = 1 h r + c l (u r -u l ) + π r -π l c r (c l + c r ) .
(2.50)

Next, using the invariants in (2.45) that involve a, b, we get

a * l = a l h l h * l , a * r = a r h r h * r , (2.51
)

and b * l = b l + h l a l c al (c al + c ar ) π ⊥l -π ⊥r + c ar (v r -v l ) , b * r = b r + h r a r c ar (c al + c ar ) π ⊥r -π ⊥l + c al (v r -v l ) .
(2.52)

Using the previous formulas one can compute the speeds,

Σ 1 = u l -c l /h l , Σ 2 = u * -c al /h * l , Σ 3 = u * , Σ 4 = u * + c ar /h * r , Σ 5 = u r + c r /h r .
(2.53)

Finally, the intermediate values for ê defined by (2.24) are given by

ê * l = (π * ) 2 2c 2 l + (π * ⊥ ) 2 2c 2 al + e l - (π l ) 2 2c 2 l - (π ⊥l ) 2 2c 2 al , ê * * l = (π * ) 2 2c 2 l + (π ⊥l ) 2 2c 2 al + e l - (π l ) 2 2c 2 l - (π ⊥l ) 2 2c 2 al , ê * r = (π * ) 2 2c 2 r + (π * ⊥ ) 2 2c 2 ar + e r - (π r ) 2 2c 2 r - (π ⊥r ) 2 2c 2 ar , ê * * r = (π * ) 2 2c 2 r + (π ⊥r ) 2 2c 2 ar + e r - (π r ) 2 2c 2 r - (π ⊥r ) 2 2c 2 ar , (2.54)
where e l , e r stand for the values of e defined by (2.22) on the left and right states respectively.

3 Analysis of the energy inequality and choice of the speeds

Sufficient stability conditions for a fixed intermediate state

In this subsection, we derive sufficient conditions for the inequality (2.31) in Lemma 2.2 to hold, for a fixed intermediate state

U * = (h * ,h * u * ,h * v * , h * a * , h * b * ), (3.1)
where the star '*' stands for any of "l**", "l*", "r*", "r**". These values are completed with those of the relaxation variables π * , π * ⊥ , c, c a , and of ê * . Note that the notation '*' differs here slightly from the one in the previous paragraph, in particular, v * and π * ⊥ do not coincide with the values in (2.49) in case of the "l**" and "r**" states. It is convenient to denote by the subscript 'l/r' any data evaluated on the initial state on the same side of the central wave as the intermediate state considered. The short notation c, c a will mean that these quantities are evaluated locally, i.e. at the '*' state, or equivalently at the 'l/r' location since these quantities are strong invariants for the central wave.

In order to analyze (2.31), we use the same strategy as in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF]Lemma 2.20], and [START_REF] Bouchut | A multiwave approximate Riemann solver for ideal MHD based on relaxation I -Theoretical framework[END_REF]Lemma 3.2], that consists in using a decomposition of e(U * )ê * in elementary entropy dissipation terms along each waves. This idea was introduced in [START_REF] Bouchut | Entropy satisfying flux vector splittings and kinetic BGK models[END_REF] in the case of constant speeds.

Lemma 3.1. With the preceding notations, we have the identity

e(U * ) -ê * = D --U * , u * - π * c + D -U * , v * - π * ⊥ c a + D 0 U * , U l/r + D + U * , v * + π * ⊥ c a + D ++ U * , u * + π * c , (3.2)
where we set for any state U and any scalar Λ, with P and P ⊥ defined in (1.11),

D --(U,Λ) = 1 4 (u - P c ) 2 - 1 4 Λ 2 (3.3) -(-P,u) • 1 2c (u - P c -Λ), 1 2 (u - P c -Λ) , D ++ (U,Λ) = 1 4 (u + P c ) 2 - 1 4 Λ 2 (3.4) -(-P,u) • - 1 2c (u + P c -Λ), 1 2 (u + P c -Λ) , D -(U,Λ) = 1 4 (v - P ⊥ c a ) 2 - 1 4 Λ 2 (3.5) -(b,v) • ha 2c a (v - P ⊥ c a -Λ), 1 2 (v - P ⊥ c a -Λ) , D + (U,Λ) = 1 4 (v + P ⊥ c a ) 2 - 1 4 Λ 2 (3.6) -(b,v) • - ha 2c a (v + P ⊥ c a -Λ), 1 2 (v + P ⊥ c a -Λ) , D 0 (U 1 ,U 2 ) =e(U 1 ) - P 2 1 2c 2 - P 2 1⊥ 2c 2 a -e(U 2 ) + P 2 2 2c 2 + P 2 2⊥ 2c 2 a -(-P 1 , b 1 ) •   1/h 1 + P 1 /c 2 -1/h 2 -P 2 /c 2 b 1 + h 1 a 1 c 2 a P 1⊥ -b 2 - h 2 a 2 c 2 a P 2⊥   . (3.7)
Proof. We have to sum up the contributions in the right-hand side of (3.2). We look first at terms that are in factor of P * ≡ P (U * ),

1 2c u * - P * c -u * + π * c in D --, (3.8
)

-1 2c u * + P * c -u * - π * c in D ++ , (3.9)
1 h * + P * c 2 - 1 h l/r + P l/r c 2 in D 0 . (3.10)
Then, since 1 h + π c 2 is a strong invariant associated to the eigenvalue u,

1 h l/r + P l/r c 2 = 1 h * + π * c 2 . (3.11)
Thus the sum of (3.8), (3.9), (3.10) equals zero. Then, the terms in factor of u * in D --, D ++ also sum up to zero, as well as the terms in factor of v * in D -, D + . Then we look at the terms in factor of -b * ,

h * a * 2c a v * - P * ⊥ c a -v * + π * ⊥ c a in D -, (3.12
) (u * -

-h * a * 2c a v * + P * ⊥ c a -v * - π * ⊥ c a in D + , ( 3
P * c ) 2 - 1 4 (u * - π * c ) 2 + 1 4 (u * + P * c ) 2 - 1 4 (u * + π * c ) 2 = (P * ) 2 2c 2 - (π * ) 2 2c 2 . (3.16)
Summing up the first lines from D -and D + we get

1 4 (v * - P * ⊥ c a ) 2 - 1 4 (v * - π * ⊥ c a ) 2 + 1 4 (v * + P * ⊥ c a ) 2 - 1 4 (v * + π * ⊥ c a ) 2 = (P * ⊥ ) 2 2c 2 a - (π * ⊥ ) 2 2c 2 a .
(3.17) The last terms are those from the first line of D 0 , e(U * ) -

(P * ) 2 2c 2 - (P * ⊥ ) 2 2c 2 a -e(U ) - P 2 2c 2 - P 2 ⊥ 2c 2 a l/r . (3.18)
Moreover, according to (2.24), ê -π 2 2c 2 -

π 2 ⊥ 2c 2
a is a strong invariant associated to the eigenvalue u, which gives e(U ) -P 2 2c 2 -

P 2 ⊥ 2c 2 a l/r = ê * - (π * ) 2 2c 2 - (π * ⊥ ) 2 2c 2 a . (3.19)
Summing up (3.16), (3.17) and (3.18), we get e(U * )ê * , which proves the lemma.

We notice now that all dissipations excepted the central one in Lemma 3.1 can be written as opposite of squares,

D --(U,Λ) = - 1 4 (u - P c -Λ) 2 , D ++ (U,Λ) = - 1 4 (u + P c -Λ) 2 , (3.20)
D -(U,Λ) = - 1 4 (v - P ⊥ c a -Λ) 2 , D + (U,Λ) = - 1 4 (v + P ⊥ c a -Λ) 2 . (3.21)
Thus they are all nonpositive, and in order to obtain a sufficient stability condition via Lemma 2.2 and Lemma 3.1, we need only to prove that D 0 from (3.7) is nonpositive.

In order to analyze D 0 (U * ,U l/r ), we group the terms in factor of 1/2c 2 , 1/2c 2 a and the terms where b is involved, because they are squares, which gives using the expression of P ⊥ from (1.11) and the relation h * a * = (ha) l/r ,

D 0 (U * , U l/r ) = gh * 2 - gh l/r 2 + (a * ) 2 2 - a 2 l/r 2 + P * ( 1 h * - 1 h l/r ) (3.22) + 1 2c 2 (P * -P l/r ) 2 + 1 2c 2 a (P * ⊥ -P ⊥l/r ) 2 - 1 2 (b * -b l/r ) 2 .
Next, we use the expression of P and P ⊥ from (1.11) to obtain

D 0 (U * , U l/r ) = gh * 2 - gh l/r 2 + g(h * ) 2 2 ( 1 h * - 1 h l/r ) (3.23) + (a * ) 2 2 - a 2 l/r 2 -h * (a * ) 2 ( 1 h * - 1 h l/r ) + 1 2c 2 g(h * ) 2 2 - g(h l/r ) 2 2 2 + (ha) 2 l/r (a * -a l/r ) 2 -2(ha) l/r (a * -a l/r ) g(h * ) 2 2 - g(h l/r ) 2 2 - 1 2 1 - (ha) 2 l/r c 2 a (b * -b l/r ) 2 .
Obviously, under the condition c a ≥ h l/r |a l/r |, the last line is nonpositive. The following result is a particular case of [8, Lemma 3.3].

Lemma 3.2. For any h * > 0, h l/r > 0 one has

gh * 2 - gh l/r 2 + gh 2 * 2 ( 1 h * - 1 h l/r ) + 1 2 1 (gh 3 ) * ,l/r gh 2 * 2 - gh 2 l/r 2 2 ≤ 0, (3.24) with (gh 3 ) * ,l/r = sup h∈[h l/r ,h * ] gh 3 . (3.25)
Using the inequality (3.24), we get an upper bound on the first line of the right-hand side of (3.23),

D 0 (U * , U l/r ) ≤ - 1 2 1 (gh 3 ) * ,l/r g(h * ) 2 2 - g(h l/r ) 2 2 2 (3.26)
- (ha) 2 l/r 2 ( 1 h * - 1 h l/r ) 2 + 1 2c 2 g(h * ) 2 2 - g(h l/r ) 2 2 2 + (ha) 4 l/r ( 1 h * - 1 h l/r ) 2 -2(ha) 2 l/r ( 1 h * - 1 h l/r ) g(h * ) 2 2 - g(h l/r ) 2 2 .
Finally, we can rewrite the right-hand side of (3.26) as a quadratic form, Multiplying by (gh 3 ) * ,l/r we get (3.28), which concludes the proof.

D 0 (U * , U l/r ) ≤ - 1 2(gh 3 ) * ,l/r - 1 2c 2 g(h * ) 2 2 - g(h l/r ) 2 2 2 - 1 2 (ha) 2 l/r 1 - (ha) 2 l/r c 2 ( 1 h * - 1 h l/r ) 2 - (ha) 2 l/r c 2 ( 1 h * - 1 h l/r ) g(h * ) 2 2 - g(h l/r ) 2 2 . ( 3

Choice of signal speeds

In this subsection we derive explicit values for the signal speeds c, c a that are sufficient for having positivity of height and (3.28), (3.29), that yield the energy inequality. Such values have been found in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF]Proposition 2.18] for Euler equations and in [START_REF] Bouchut | A multiwave approximate Riemann solver for ideal MHD based on relaxation II -Numerical implementation with 3 and 5 waves[END_REF] for MHD equations. We use here the approach of [START_REF] Bouchut | A multiwave approximate Riemann solver for ideal MHD based on relaxation II -Numerical implementation with 3 and 5 waves[END_REF], that enables to treat negative pressures π. We make the following a priori choice of the relaxation speeds c l , c r ,

c l = h l s l + 3 2 h l (u l -u r ) + + (π r -π l ) + h l s l + h r s r , c r = h r s r + 3 2 h r (u l -u r ) + + (π l -π r ) + h l s l + h r s r , (3.33)
with s = a 2 + gh. (3.34) This choice for c l and c r implies in particular that c l ≥ h l s l and c r ≥ h r s r . We set

X l = 1 s l (u l -u r ) + + (π r -π l ) + h l s l + h r s r , ( 3
= 1 h l + c r (u r -u l ) + π l -π r c l (c l + c r ) (3.37) ≥ 1 h l - c r (u l -u r ) + c l (c l + c r ) - (π r -π l ) + c l (c l + c r ) (3.38) ≥ 1 h l - (u l -u r ) + c l - (π r -π l ) + c l (h l s l + h r s r ) . ( 3
1 h * l ≥ 1 h l   1 - X l 1 + 3 2 X l    . (3.40)
Thus h * l > 0, and 0 < h * l ≤ h l /x l , (3.41)

with x l = 1 - X l 1 + 3 2 X l ∈ (1/
A(U l , U r ) ≤ C(|u l | + |u r | + s l + s r ) (4.10)
with C an absolute constant, bounding the propagation speed of the approximate Riemann solver by the left and right true speeds. This property is also valid in [START_REF] Bouchut | A new model for shallow viscoelastic fluids[END_REF] for shallow water elastic fluids and is more general than the possibility of treating data with vacuum considered in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF]. Note that for getting (4.10), no restriction on the ratio h l /h r is required.

We are now able to obtain our main result.

= Σ 3 = Σ 4 .
The item (vi) has been proved with Lemma 3.6. The item (vii) is obvious with (3.33), (3.52). The property that h, hu, hv are conservative is also obvious. Thus it remains to prove the statement of (viii) concerning the consistency. Denote as in the introduction b = (a,b), u = (u,v). Then the magnetic equations of the original system (1.9), (1.10) can be written

∂ t (hb) + ∂ x (hbu -hau) + u∂ x (ha) = 0.
(4.12)

The numerical values for Test 2 are given in Table 2. The reference solution is plotted in Figure 5. It consists of, from left to right, a left rarefaction wave, a left Alfven contact, a material contact, a right Alfven contact, and a right shock. We observe in Figure 6 that the second order resolution improves the sharpness of contact discontinuities, but sometimes gives rise to slight instabilities.

The numerical values for Test 3 are given in Table 3. The reference solution is plotted in Figure 7. It consists of, from left to right, a left rarefaction wave, a left Alfven contact, a right rarefaction wave. There is no material contact, nor right Alfven contact. The numerical results are shown in Figure 8. We notice that on the right the height h vanishes, and the variables u, v, a, b take eventually non-physical values, due to the fact that only the conservative variables hu, hv, ha, hb make sense. Taking this into account, we observe that the computed solution achieves a good accuracy.

Figure 1 :

 1 Figure 1: Intermediate states in the Riemann solution

But since b + ha c 2 aπ 2 aP 15 )

 2215 ⊥ is a strong invariant associated to the eigenvalue u, one has b l/r + (ha) l/r c ⊥l/r = b * + h * a * c Thus we get that the sum of (3.12), (3.13), (3.14) equals zero. Now it remains to sum up the first lines from (3.3)-(3.7). Summing up the first lines from D --and D ++ we get 1 4

  the right-hand side and simplifying, we obtain (ha) 2 l/r (gh 3 ) * ,l/r ≤ c 2 (gh 3 ) * ,l/r -1. (3.32)

1 Table 1 : 1 Table 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :

 11123456 Figure 3: Riemann solution for Test 1 at time t = 0.1 computed with the HLL flux or the 5-wave solver with 10000 points

Figure 7 :Figure 8 :

 78 Figure 7: Reference solution for Test 3 at time t = 0.05 computed at first order with 10000 points

  Since (2.25) is a conservative equation, and its conserved quantity and flux reduce to E and G on the left and right states, integrating it over rectangles one gets

and the discrete energy inequality

(2.21)

holds under the CFL condition

(2.13)

.

Proof.

  ) * ,l/r is defined in(3.25), and (ha) l/r , c, c a are evaluated on the same side of the central wave as is U * .Proof. We use Lemma 3.1 and the previous analysis. With the condition (3.29), it yields (3.27). It is now sufficient that the right-hand side in (3.27) is a non positive quadratic form. Thus it is sufficient that

		.27)
	It leads to the following proposition stating that the entropy condition reduces
	to subcharacteristic conditions.
	Proposition 3.3. In order to have e(U * ) -ê * ≤ 0 at the intermediate state U *
	(ensuring the discrete energy inequality according to Lemma 2.2), it is enough
	that h * > 0 and
	(ha) 2 l/r + (gh 3 ) * ,l/r ≤ c 2	(3.28)
	and
	(ha) 2 l/r ≤ c 2 a ,	(3.29)
	where (gh 3 (ha) 2 l/r

  Theorem 4.1. Assume that the initial data U l , U r satisfy h l > 0, h r > 0 and use the choice of the speeds (3.33),(3.52) with π l = P (U l ), π r = P (U r ). Under the CFL condition (4.8), (4.9), the Riemann solver defined by the intermediate states and speeds Σ i computed in Subsection 2.3 and the numerical fluxes F l (U l ,U r ), F r (U l ,U The numerical viscosity is sharp, in the sense that the propagation speeds Σ i of the approximate Riemann solver tend to the exact propagation speeds when the left and right states U l , U r tend to a common value, (viii) The variables h, hu, hv are conservative, and the scheme satisfies an asymptotic consistency with the non-conservative part of the system for smooth data.Proof. The items (i), (ii) are consequences of Proposition 3.5 and Lemma 2.2. The items (iii), (iv) hold according to Lemma 2.1. The proof of (v) is straightforward with the formulas of Section 2.3, since in this case Σ 2

	0,	(4.11)
	are resolved exactly,
	(vi) Data with bounded propagation speeds give finite numerical propagation
	speed, according to (4.10),
	(vii)

r ) defined via (4.1)-(4.7), has the following properties. (i) It keeps the positivity of h, (ii) The approximate Riemann solver satisfies the discrete energy inequality (2.21), with numerical energy flux (2.30), (iii) Isolated material contact discontinuities are resolved exactly, (iv) Isolated Alfven contact discontinuities are resolved exactly, (v) Resonant material-Alfven contact discontinuities defined by h = cst, u = cst, a =

The assumptions concerning the pressure (3.44)-(3.45) are satisfied with p(ρ) = gρ 2 /2 and α = 3/2. Therefore we can apply this result with ρ = h l , X = X l , which gives

(3.47) Thus, using (3.43) and (3.47), for getting (3.28) it is enough that

(3.48)

But using (3.36) and the definition of s l from (3.34),

which yields (3.48). The same analysis is valid on the right, with Proof. We apply Proposition 3.3. The above computations show that h * l ,h * r > 0 and that (3.28) holds. The choice (3.52) gives obviously (3.29). Lemma 3.6 (Bounds on the propagation speeds). The formulas (3.33) ensure the following estimate on propagation speeds:

where C > 0 is an absolute constant.

Proof. We have |P | ≤ hs 2 with P = g h 2 2ha 2 and s 2 = a 2 + gh. Since π l = P l and π r = P r , the result follows obviously.

Properties of the relaxation approximate Riemann solver

Before stating our main result, we have to explicit the numerical fluxes and the CFL condition. The solution to the Riemann problem for the relaxation system (2.1)-(2.9) has five wave speeds Σ 1 < Σ 2 < Σ 3 < Σ 4 < Σ 5 , that can be computed by (2.53). The intermediate states l * * , l * ,r * ,r * * have been determined in Section 2.3. We would like now to compute the left/right numerical fluxes F l , F r that are involved in the update formula (2.15). All components of the system except ha and hb are conservative, thus classical computations give the associated numerical fluxes,

where the conservative part involves the Riemann solution evaluated at x/t=0,

More explicitly (4.2) yields that the quantities between parentheses are evalu-

As usual there is no ambiguity when equality occurs in these conditions. Similarly, the numerical energy flux is computed according to (2.30).

We complete these formulas by computing the left/right numerical fluxes for the variables ha and hb from (2.16),

Using the computation performed in [9, Subsection 5.3], these fluxes can be also written

If Σ 3 ≤ 0 then

where v * is the central value of v defined in (2.49) (and indeed Σ 3 = u * ). The maximal propagation speed is then

The CFL condition (2.13) becomes

At the discrete level, the hb components of the update formula (2.15) involve the numerical fluxes from (4.1)

Using the formulas for the fluxes (4.3) and (4.6), (4.7), we observe that

Thus, asymptotically when U l , U r → U , one has

Since obviously F hb l (U,U ) = F hb r (U,U ) = hbuhau, we conclude the asymptotic consistency with (4.12) for smooth solutions (in the sense proposed in [6, Section 4.2]), which concludes the proof of (viii).

Numerical tests

In this section we perform numerical computations in order to evaluate the properties of the scheme, in relation with Theorem 4.1. First and second order methods in time and space are evaluated, the latter using an ENO reconstruction, as described in [6, section 4.13]. The conservative variable is U as in (1.13), and the slope limitations are performed on the variables h, u, v, ha, b.

We take 200 points, and plot a reference solution obtained by a first order computation with 10000 points. The CFL-number is taken 1/2 in all tests. The space variable x is taken in [0,1], g = 9.81. Neumann boundary conditions are applied. Three test cases are investigated:

• Test 1 is a Riemann problem with (ha) l = (ha) r ,

• Test 2 is a generic Riemann problem with positive height,

• Test 3 is a Riemann problem where the initial heights are taken positive on the left side and vanishing on the right side.

The numerical values for Test 1 are given in Table 1. In this case, ha = 1/2 remains constant, in accordance with (1.9). The first order method in time and space is evaluated, with our numerical fluxes defined by (4.1)-(4.7). Note that ha remains stricly constant at the discrete level, because of (4.3). Our results are compared to those obtained with the HLL flux, see for example [6, Equation (2.111)], applied here to the system (1.6)-(1.8), (1.10) with ha = 1/2, which is conservative. We ommit to plot the ha component, since it is constant in both methods. The reference solution is plotted in Figure 3. It consists of, from left to right, a left rarefaction wave, a left Alfven contact, a right Alfven contact, a right shock. There is no material contact. We observe as expected that the components h, u only vary through fast waves whereas the components v and b only vary through Alfven waves. In Figure 4 we observe that the performances on fast waves are extremely similar for both methods, whereas for Alfven waves, our 5-wave solver shows a much better resolution than the 2-wave HLL solver.


Florence Alberge
email: florence.alberge@u-psud.fr

Min-sum decoding of irregular LDPC codes with adaptive scaling based on mutual information

Keywords: Parity check codes, iterative decoding, linear approximation, mutual information, adaptive LLR scaling

Min-Sum decoding (MS) is an alternative to belief propagation decoding with substantially lower complexity. MS often results in an overestimation of the log likelihood ratio (LLR) in particular in the early stage of the iterative process. A linear post-processing is usually performed as a compensation. With regular low density parity check codes (LDPC), a fixed scaling of the LLR yields sufficiently good results. In contrast, adaptive strategies are mandatory with irregular codes. It is well known that the scaling factor is an increasing function of the reliability of the LLR. In most of the publications, the scaling factor is envisioned as a function of both the iteration number and the signal-to-noise ratio. It is proposed here to use the mutual information between extrinsics as a measure of the reliability of the LLR. A practical implementation is derived with reasonable complexity. Compared to the literature, the proposed method yields slightly better results in terms of BER and a significant reduction in the number of iterations.

I. INTRODUCTION

LDPC were introduced in the literature by Gallager [START_REF] Gallager | Low-density parity-check codes[END_REF] in the early sixties. Limited computing resources prevented Gallager from demonstrating the capabilities of message-passing decoders even for moderate code length. This work was ignored and rediscovered 30 years later by MacKay [START_REF] Mackay | Good error-correcting codes based on very sparse matrices[END_REF]. Combined with belief propagation (BP) decoding, LDPC codes exhibit performance near to the Shannon limit. The computational complexity of BP, especially at check-node, is not suitable for hardware implementation. Complexity reduction can be obtained through an approximation of the check-node calculation leading to min-sum decoding [START_REF] Fossorier | Reduced complexity iterative decoding of low-density parity check codes based on belief propagation[END_REF]. This efficient implementation suffers performance degradation which can be compensated by linear post-processing of the messages computed at the check nodes [START_REF] Chen | Density evolution for two improved BP-Based decoding algorithms of LDPC codes[END_REF]. In this seminal work, a single correcting function is used during the whole iterative process and independently of the check-node. A density evolution analysis is often conducted in order to find the optimal value of the coefficients of the linear function [START_REF] Chen | Density evolution for two improved BP-Based decoding algorithms of LDPC codes[END_REF]. This approach is not efficient enough with many irregular codes. It is observed in [START_REF] Zhang | Two-dimensional correction for min-sum decoding of irregular LDPC codes[END_REF] that considering different post-processing depending on the degree of the nodes can improve the overall performance. Several authors suggest to consider in addition an adaptive rule for updating the linear function within the iterative process. It is proved in [START_REF] Lechner | Improved Sum-Min Decoding for Irregular LDPC Codes[END_REF] that good results could be obtained with irregular codes by scaling the message at the check nodes with a scaling factor which is a function of both the degree of the node and of the mutual information. Since the mutual information is expected to increase during the iterative process, the value of the scaling factor should not be fixed. The practical implementation of the method is not solved in [START_REF] Lechner | Improved Sum-Min Decoding for Irregular LDPC Codes[END_REF], a sub-optimal approach is derived in which the scaling factor is not anymore a function of two variables but is computed as a function of the sole iteration number. Recently, a variable scaling scheme based on the generalized mutual information was proposed in [START_REF] Xu | Variable LLR Scaling in Min-Sum Decoding for Irregular LDPC Codes[END_REF]. The method select the scaling factors as per iteration and as per check node degree independently overcoming the multi-dimensional issue faced by previous methods. We propose in this paper to make a step further. Here, the variable scaling factors are adjusted with respect to both check node degrees and mutual information with a particular focus on implementation aspects that are missing in [START_REF] Lechner | Improved Sum-Min Decoding for Irregular LDPC Codes[END_REF]. It will be proved that the method can be implemented at low cost. Comparison with previously mentioned methods is also provided.

II. ITERATIVE DECODING

The iterative decoding of LDPC codes consists in an exchange of information between the variable and check nodes of the code. Let assume a (N, K) LDPC code over an AWGN channel. Let R i denote the log-likelihood ratio (LLR) of bit i derived from the observations. At iteration k, let

E (k) ji and M (k)
ji be the LLR of bit i which are sent from check node j to variable node i and from variable node i to check node j respectively. The LLR E (k) ji and M (k) ji are defined if and only if the variable node i and the check node j are connected in the Tanner graph of the code. The standard BP algorithm proceeds in two steps [START_REF] Gallager | Low-density parity-check codes[END_REF]. The horizontal step reads:

E (k) ji = 2tanh -1 i =i tanh M (k-1) ji 2 (1)
and is followed by a vertical step as :

M (k) ji = R i + j =j E (k) j i (2)
In min-sum (MS), the horizontal step is replaced by a simpler updating rule by modifying (1) into

E (k) ji = i =i sgn(M (k-1) ji ) min i =i |M (k-1) ji | (3)
This approximation degrades the performance of the iterative process and yields over-estimated LLR especially in early iterations. From a statistical point of view, exchanged messages are usually modelled using Gaussian approximation [START_REF] Brink | Convergence behavior of iteratively decoded parallel concatenated codes[END_REF] that can be seen as noisy estimates of the desired message X. It is assumed here that X ∈ {-1, 1}. Let L denote a LLR. Under Gaussian approximation L reads:

L = σ 2 2 X + σn (4)
where n N (0, 1) and σ ∈ R + \ {0}. In particular, the LLR that proceeds from the received data at the output of an AWGN channel can be modelled by [START_REF] Chen | Density evolution for two improved BP-Based decoding algorithms of LDPC codes[END_REF]. The validity of this assumption has been questioned in [START_REF] Fu | On gaussian approximation for density evolution of low-density parity-check codes[END_REF] for both regular and irregular LDPC codes and under AWGN channel. The main conclusion is that it is inappropriate to use a single parameter σ to model the statistical behaviour of L. Alternatively, the mean and the variance of the density should be considered as independent variables leading to the following model:

L = α σ 2 2 X + σn (5)
where

α ∈ R + \ {0}. Let G α (X) denote the set of LLR 1 following (5), then αL ∈ G 1 (X) where G 1 (X)
is the set of LLR following (4). LLR scaling is thus equivalent to propagating messages which are elements of G 1 (X). From theoretical point of view [START_REF] Zhang | Two-dimensional correction for min-sum decoding of irregular LDPC codes[END_REF], [START_REF] Xu | Variable LLR Scaling in Min-Sum Decoding for Irregular LDPC Codes[END_REF], [START_REF] Lechner | Improved Sum-Min Decoding for Irregular LDPC Codes[END_REF], the scaling factor should depend on the check node degree. The rationale behind that is the relationship between the check node degree and the reliability.

It is well known that the scaling factor should be an increasing function of the reliability of the LLR and that the reliability is larger with low degrees. Putting all these facts together, the vertical step in min-sum decoding is replaced with

M (k) ji = R i + d∈D(i) α (k) (d) j ∈M(d) j =j E (k) j i (6)
where D(i) is the set of degrees corresponding to parity-check equations in which bit i is involved and where M(d) is the set of indexes of parity check equations with degree d. In eq. ( 6),

j ∈M(d) j =j E (k) j i is scaled in order to be in G 1 (X). Since G 1 (X) is stable under addition, M (k)
ji ∈ G 1 (X). As a consequence, the horizontal step remains unchanged since

{E (k)
ji } is already scaled in [START_REF] Lechner | Improved Sum-Min Decoding for Irregular LDPC Codes[END_REF] with {α (k) (d)}. It has been previously mentioned that the numerical value of {α (k) (d)} depends on the reliability of the LLR. A natural measure of the reliability is the mutual information I(L, X). Since X is unknown at the receiver side, an accurate estimation of I(L, X) is unlikely to be obtained at the decoder. It was observed in [START_REF] Alberge | On Some Properties of the Mutual Information between Extrinsics with Application to Iterative Decoding[END_REF] for both turbo-codes and regular LDPC codes that the scaling factor can be obtained as a function of a single parameter: the mutual information between extrinsics. An efficient computation scheme for finite sequences is also provided leading to an estimation of the mutual information between extrinsics at the receiver side. In the next section, the results in [START_REF] Alberge | On Some Properties of the Mutual Information between Extrinsics with Application to Iterative Decoding[END_REF] are briefly reviewed and extended to irregular LDPC codes. The method for computing α (k) (d) from the mutual information between extrinsics is also detailed. The complexity of the proposed method is addressed in section IV in order to be compatible with a low-cost implementation. Numerical results are provided in section V.

III. MUTUAL INFORMATION BETWEEN EXTRINSICS

Let L y , L z ∈ G α (X), the mutual information I(L y , L z ) is defined as:

I(L y , L z ) = E p Ly ,Lz ( y , z ) log 2 p Ly,Lz ( y , z ) p Ly ( y )p Lz ( z ) (7)
The mutual information between extrinsics is related to I(L y , X) and I(L z , X) through [START_REF] Alberge | On Some Properties of the Mutual Information between Extrinsics with Application to Iterative Decoding[END_REF]:

I(L y , L z ) = I(L y , X) + I(L z , X) -I(L y + L z , X) (8)
where

I(L, X) = 1 - +∞ -∞ p L ( |X = 1)log 2 (1 + e -α
)d . and L stands for L y or L z . The evolution of I(L y , L z ) for a given system can be obtained offline via an histogram method applied to each of the three terms involved in [START_REF] Brink | Convergence behavior of iteratively decoded parallel concatenated codes[END_REF] as in EXIT Charts. Let L y resp. L z denote the random variable associated with { j E j,i } resp. { j M j,i }. In min-sum decoding, E j,i is the LLR of the probability that bit i causes parity check j to be satisfied. On the other side, M ji is the LLR of the probability of bit i based on the knowledge of the received sequence and of the probability of the other bits in paritycheck j. Thus, j E j,i is the opinion of the parity-check nodes on bit i whereas j M j,i is the opinion of the variable nodes on bit i. Figure 1 depicts the extension of the singlecurve Exit-Chart presented in [START_REF] Alberge | On Some Properties of the Mutual Information between Extrinsics with Application to Iterative Decoding[END_REF] to irregular LDPC codes. The principle is quite simple. All messages entering a vertical or a horizontal node belong to G 1 (X). The scaling factor α(d) is computed easily from {E ji } j∈M(d) since X is known. Then {E ji } j∈M(d) is scaled and the vertical step is processed with [START_REF] Lechner | Improved Sum-Min Decoding for Irregular LDPC Codes[END_REF]. In parallel, I(L y , L z ) is computed with (8) and an histogram method. From this scheme, the evolution of α(d) can be plotted as a function of I(L y , L z ). We consider as an example, an irregular code of rate 12 and length 10032 with degree distribution ρ(x) = 0.25x 4 + 0.75x 14 and λ(x) = 0.1917x + 0.0125x 2 + 0.0417x 4 + 0.0750x 5 + 0.2040x 6 + 0.0667x 7 + 0.225x 8 + 0.0833x 9 + 0.058x 10 + 0.0541x 12 . This LDPC code was previously considered in [START_REF] Xu | Variable LLR Scaling in Min-Sum Decoding for Irregular LDPC Codes[END_REF] and is used through this paper for comparison purpose. This code has two check node degrees, 5 and 15. The evolution of α(5) and α(15) as a function of I(L y , L z ) is given in Fig. 2. As expected, the scaling factor takes different values depending on the check-node degree. This fact was already established in [START_REF] Xu | Variable LLR Scaling in Min-Sum Decoding for Irregular LDPC Codes[END_REF] where the scaling factor is plotted as a function of

Step Horizontal (Histogram Method)

Mutual Information Evaluate I(L y , L z ) Prior ∈ G 1 (X) σ 2 g Gaussian R i + L y {E ji } Prior ∈ G 1 (X) R i + σ 2 z
Step with Scaling Compute Compute α(1) the iteration number. The shape of curves in Fig. 2 is very similar to the ones in Fig. 5 of [START_REF] Xu | Variable LLR Scaling in Min-Sum Decoding for Irregular LDPC Codes[END_REF]. The main difference is that here the value of the scaling factor is plotted as a function of I(L y , L z ) and is valid for any signal to noise ratio (SNR) and at any iteration since this information is already captured in the mutual information. In contrast, the curves given in [START_REF] Xu | Variable LLR Scaling in Min-Sum Decoding for Irregular LDPC Codes[END_REF] are obtained for a particular value of the SNR and are accurate only locally. Let I = I(L y , L z ). Using polynomial fitting and Fig. 2, we obtain the following rules:

{E ji } j∈M(1) Vertical j E j,i {α(1)E ji } j∈M(1) {α(d)E ji } j∈M(d) α(d) Gaussian j M j,i L z {M ji } α(1) {E ji } j∈M(d) α(d)
α(5) = 0.416 + 1.956I -2.975I 2 + 1.832I 3 (9) α(15) = -0.232 + 2.747I -3.690I 2 + 2.268I 3 (10)
The next section is devoted to the estimation of I(L y , L z ) at the receiver side with finite length sequences.

ONLINE COMPUTATION OF I(L y , L z )

Let consider a message X = (x 1 , x 2 , ..., x K ) with finite length K. Let y,i resp. z,i be a realization of random variable L y resp. L z . Both y,i and z,i are noisy estimates of x i with x i ∈ {-1, +1}. An estimator I yz of I(L y , L z ) is given in [START_REF] Alberge | On Some Properties of the Mutual Information between Extrinsics with Application to Iterative Decoding[END_REF] with expression:

I yz = 1 + 1 K K k=1 log 2 (s k ) (11) with s k = x k e x k ( y,k + z,k ) (1+e x k y,k )(1+e x k z,k )
. An equivalent formulation is given below:

I yz = 1 + 1 K K k=1 log 2 1 + e y,k + z,k (1 + e y,k )(1 + e z,k ) (12)
The accuracy of this estimator was evaluated in [START_REF] Alberge | On Some Properties of the Mutual Information between Extrinsics with Application to Iterative Decoding[END_REF] for a serially concatenated turbo-code. Obviously, the distance between I yz and I(L y , L z ) depends of the system and also of the blocklength K. A numerical evaluation is provided in Fig. 3 for the irregular LDPC code defined in the previous section. The mean square error (MSE) is plotted as a function of I(L y , L z ) where I(L y , L z ) is evaluated with the histogram method whereas I yz is obtained with (12). It can be observed that the histogram method is computationally demanding and necessitates the knowledge of X. This is an offline tool only. In contrast, I yz can be computed offline an online as well. From Fig. 3 we can conclude that, for the irregular code under consideration, the MSE is lower than 6.10 -4 whatever the value of the mutual information is. On average, |I(L y , L z ) -I yz | ≤ 2.5.10 -2 which gives an idea of the accuracy of the proposed estimator.

The implementation of the method in min-sum decoding is considered below. Let f ( ) = log 1 + e then straightforward to prove that

I yz = 1 + 1 Klog(2) K k=1 f ( y,k + z,k ) -f ( y,k ) -f ( z,k ) (13) Using the fact that log(e a +e b ) = max(a, b)+log(1+e -|a-b| ), f ( ) = max(0, ) + log(1 + e -| | ) def = max(0, ) + g( ). It is
max(0, y,k + z,k ) -max(0, y,k ) -max(0, z,k ) = 0 if sign( y,k ) = sign( z,k ) -min(| y,k |, | z,k |) otherwise (14)
To limit the computational complexity of the iterative decoding, it is proposed to use pre-computed values of g(.) for

| q | ∈ Q with Q = {0, q, 2q, ..., L max } and where q = Q( ) is the closest value to in Q. Let ∆ = g( ) -g( q )
denote the corresponding error in the evaluation of g(.). Two different values of ∆ are considered. The first one is ∆ = 10 -2 which is the same order of magnitude than the maximal average error obtained with unquantized LLR in Fig. 3. In that case L max = 4.6 and card(Q) = 115. Alternatively, ∆ = 10 -3 is also considered. This last choice leads to L max = 6.9 and card(Q) = 1725. Let I Q yz denote I yz evaluated with quantized values for g(.). The reference is still the value of the mutual information provided by the histogram method (M SE(I(L y , L z ) -I Q yz )). The MSE of I Q yz is depicted in Fig. 3 for ∆ = 10 -2 and ∆ = 10 -3 . The MSE of I Q yz when ∆ = 10 -3 is almost equal to the MSE of I yz . Surprisingly, the quantization tends to lower the error when ∆ = 10 -2 and I(L y , L z ) < 0.85. In any case, the MSE due to both the estimator and the quantization effect remains lower than 6.10 -4 . In the same way, α(5) and α(15) can be computed in advance for quantized values of I Q yz . The extra cost of the proposed method compared to fixedscaling factor techniques is (i) the storage of the table containing the pre-determined values of α(5), α(15) and g(.), (ii) the arithmetical complexity due to the computation of I Q yz which has to be performed at each iteration. Point (i) is directly connected with the desired precision. With ∆ = 10 -2 , 115 different value of g(.) have to be stored and if an extra quantization error of 10 -2 is allowed on I Q yz then 100 values have to be stored for α(5), α(15). The arithmetical complexity in (ii) for I Q yz is evaluated in table I and compared to the min-sum decoding in eq. [START_REF] Mackay | Good error-correcting codes based on very sparse matrices[END_REF][START_REF] Fossorier | Reduced complexity iterative decoding of low-density parity check codes based on belief propagation[END_REF]. In table I, c(i) denotes the number of ones in column i of the parity-check matrix of size M × K, r(i) denotes the number of ones in raw i of the parity-check matrix and U is the number of bit positions in which sign( y,k ) = sign( z,k ). The comparison is conducted for a single iteration and for the main operations (addition and minimum or maximum selection). For the irregular LDPC Horizontal step (MS) Vertical step (MS)

I Q yz Addition 0 K i=1 c(i) 2 3K -U + 1 Min./Max. M j=1 r(i)(r(i) -2) 0 K -U selection Table I COMPLEXITY OF I Q yz COMPARED TO MS
defined in section III, numerical values can be obtained. We have

3K-U +1 K i=1 c(i) 2 ≤ 3K K i=1 c(i) 2 = 8.44% and K-U M j=1 r(i)(r(i)-2) ≤ K M j=1 r(i)(r(i)-2) = 2%
. This leads to the conclusion that the extra cost for the evaluation of the mutual information is negligible compared to the total cost inherent to MS decoding.

V. SIMULATION RESULTS

In this section, the performance of the proposed method (labelled "MS, α = f (IM )" in the following) is evaluated in terms of bit error rate (BER) and number of iterations. The irregular LDPC code of rate 1 2 and length 10032 with degree distribution ρ(x) = 0.25x 4 + 0.75x 14 and λ(x) = 0.1917x + 0.0125x 2 + 0.0417x 4 + 0.0750x 5 + 0.2040x 6 + 0.0667x 7 +0.225x 8 +0.0833x 9 +0.058x 10 +0.0541x 12 is again considered. Comparison is provided with classical methods and also with the most recent literature. The following schemes are considered: (SP) iterative decoding with sum-product implementation, (MS, α = f (IM )) min-sum implementation with scaling factor as in [START_REF] Lechner | Improved Sum-Min Decoding for Irregular LDPC Codes[END_REF] in which the scaling factors are given by (9-10) and I = I Q yz , (MS, α = f (EbN 0 threshold , it)) min-sum implementation with scaling factor acquired using training LLR and with a value of EbN 0 corresponding to the practical threshold (EbN 0 threshold ) of the code [START_REF] Xu | Variable LLR Scaling in Min-Sum Decoding for Irregular LDPC Codes[END_REF], (DNMS) MS implementation with fixed values of α(5) and α(15) (generalization of the normalized MS to irregular codes [START_REF] Chen | Density evolution for two improved BP-Based decoding algorithms of LDPC codes[END_REF] with best scaling pair is α(5) = 0.88 and α(15) = 0.68 [START_REF] Xu | Variable LLR Scaling in Min-Sum Decoding for Irregular LDPC Codes[END_REF]).

The results are given in Fig. 4 for BPSK modulation and in Fig. 5 for 256-QAM modulation. For each implementation, the same scaling sequence (if apply) is used for both BPSK and 256QAM. We observe the following. Adaptive LLR scaling outperforms fixed LLR scaling in terms of BER. The BER curves obtained with adaptive scaling are very close to the BER curve obtained with SP implementation. The BER with the proposed method (MS, α = f (IM )) is slightly better than with the method in [START_REF] Xu | Variable LLR Scaling in Min-Sum Decoding for Irregular LDPC Codes[END_REF] (MS, α = f (EbN 0 threshold , it)) for both BPSK and 256QAM. The main gain obtained with our method compared to [START_REF] Xu | Variable LLR Scaling in Min-Sum Decoding for Irregular LDPC Codes[END_REF] lies in the reduction in the number of iterations. Indeed, as mentioned above, the scaling sequence α = f (EbN 0 threshold , it) has been acquired for a particular value of SNR EbN 0 threshold . When the true SNR is significantly above EbN 0 threshold , the reliability of the LLR is rapidly increasing and the mutual information between extrinsics as well. In that case, higher values of the scaling factors could be chosen even in first iterations. Using the rule α = f (EbN 0 threshold , it) slows down the convergence since the same rule is used whatever the value of EbN0 is. The same problem arises with DNMS. With BPSK, the iteration gain with the proposed method compared to α = f (EbN 0 threshold , it) is approximately 45% with EbN 0 = 2dB and up to 60% with higher EbN 0.

VI. CONCLUSION

In this paper, the practical implementation of adaptive LLR scaling has been investigated. The originality of the proposed method is the use of a mutual information driven rule for the scaling factor. This is an old idea but, up to our knowledge, this is the first time that an implementation is proposed with acceptable complexity. The efficiency of the proposed method has been demonstrated in terms of BER and convergence rate. It is worth mentioning that the estimate of the mutual information between extrinsics could also be used as an early stopping criterion or sent back to the transmitter via a feedback channel. The last point could be of interest in several applications including secure transmissions.

Figure 1 .Figure 2 .

 12 Figure 1. Single curve EXIT Chart for irregular LDPC codes

Figure 3 .

 3 Figure 3. M SE(I(Ly, Lz) -Iyz) with quantized and unquantized LLR.

Figure 4 .

 4 Figure 4. BER (up) and iterations (down) under BPSK AWGN channel

Figure 5 .

 5 Figure 5. BER (up) and iterations (down) under 256QAM AWGN channel

The term LLR is used indifferently for the log-ratio in (4) and (5) even if[START_REF] Zhang | Two-dimensional correction for min-sum decoding of irregular LDPC codes[END_REF] does not strictly match the definition of a log-likelihood ratio.