Emanuela Abbate
email: eabbate@uninsubria.it

Angelo Iollo
email: angelo.iollo@inria.fr

AND Gabriella Puppo
email: gabriella.puppo@uniroma1.it

AN ASYMPTOTIC-PRESERVING ALL-SPEED SCHEME FOR FLUID DYNAMICS AND NONLINEAR ELASTICITY *

Keywords: low Mach limit, all-speed schemes, asymptotic-preserving property, relaxation, non-linear elasticity, entropy production AMS subject classifications. 65M08, 65M99, 76T15

An implicit relaxation scheme is derived for the simulation of multi-dimensional flows at all Mach numbers, ranging from very small to order unity. An analytical proof of the asymptotic preserving property is proposed and the divergence-free condition on the velocity in the incompressible regime is respected. The scheme possesses a general structure, which is independent of the considered state law and thus can be adopted to solve gas and fluid flows, but also deformations of elastic solids. This is achieved by adopting the Jin-Xin relaxation technique in order to get a linear transport operator. The spatial derivatives are thus independent of the EOS and an easy implementation of fully implicit time discretizations is possible. Several validations on multi-dimensional tests are presented, showing that the correct numerical viscosity is recovered in both the fully compressible and the low Mach regimes. An algorithm to perform grid adaptivity is also proposed, via the computation of the entropy residual of the scheme.

1. Introduction. In the propagation of flows of gases and of other compressible materials, complex and non stationary phenomena are generated. The complexity is often related to the local stiffness of the involved media, to the geometry of the physical problem or to boundary and initial conditions. All these aspects can produce waves propagating at very different speeds inside the considered materials, giving rise to specific numerical problems.

The accurate simulation of these phenomena requires the construction of numerical schemes that are able to deal with different regimes. The Mach number M is given by the ratio between the flow velocity and the speed of sound, namely it measures how much the acoustic waves are faster with respect to the flow velocity. In gas and fluid dynamics, M is used to identify the regime of the considered flow, ranging from fully compressible (local Mach number of order unity or higher) to incompressible (very small local Mach number).

Standard explicit-upwind codes developed for the simulation of compressible flows usually fail in approximating fluid flows or solid deformations at low speed. This is mainly due to the excessive numerical viscosity introduced on the slow waves by an upwind Godunov-like spatial discretization when the Mach number of the flow becomes small, as it has been proved in [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF][START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit: II. Godunov type schemes[END_REF][START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF]. Specifically, an upwind scheme on a Cartesian grid leads to pressure fluctuations of order O (M ), while in the continuous case the pressure fluctuations are of order O M 2 . Moreover, when explicit methods are used to solve low Mach number flows, the time step ∆t becomes extremely small due to the enforcement of the CFL stability condition. By imposing this constraint, ∆t is roughly proportional to the Mach number M :

(1.1) ∆t ≤ ∆x λ max = ∆x max|u ± c| = M ∆x max|u (M ± 1) | ,
∆x being the space step and λ max the fastest characteristic speed. For the Euler system we have that λ max = u ± c, c being the sound speed and u the velocity.

Consequently, compressible codes require an increasingly large computational time as the incompressible regime is reached.

The derivation of all-speed solvers is motivated by all the above mentioned reasons.

In general, the purpose of an all-speed scheme is to handle both the compressible regime and the incompressible one. To this end, the derivation of preconditioning methodologies has been triggered by Chorin [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF] and by Turkel [START_REF] Turkel | Preconditioned methods for solving the incompressible and low speed compressible equations[END_REF][START_REF] Turkel | Local preconditioners for steady state and dual time-stepping, Mathematical Modeling and Numerical Analysis[END_REF], proposing a modification of the Roe matrix that moderates the numerical diffusion of upwind schemes, inside a fully implicit time discretization [START_REF] Viozat | Implicit upwind schemes for low Mach number compressible flows[END_REF]. More recent methods falling in this category have been proposed in [START_REF] Li | An all-speed Roe-type scheme and its asymptotic analysis of low Mach number behaviour[END_REF][START_REF] Barsukow | A numerical scheme for the compressible low-Mach number regime of ideal fluid dynamics[END_REF][START_REF] Van Leer | Characteristic time-stepping or local preconditioning of the Euler equations[END_REF]. However, the main problem of these techniques is related to the difficulty in handling the non-linearities of classical upwind discretizations (e.g. approximate Riemann solvers) within the fully implicit time integration.

Numerical schemes based on low Mach asymptotics have been proposed by Klainerman and Majda [START_REF] Klainerman | Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF][START_REF] Klainerman | Compressible and incompressible fluids[END_REF]. Klein adapted this technique to derive an operator splitting into convection and pressure waves [START_REF] Klein | Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics i: One-dimensional flow[END_REF]. These methods have been derived for the low Mach number regime, but they also have been the starting point for the development of different all-speeds schemes. These are mainly based on the separation of the fast and slow scales, that are respectively integrated in time explicitely and implicitely in an IMEX logic, producing stability conditions on ∆t that are independent of the Mach number. Some examples may be found in [START_REF] Boscarino | All Mach number second order semi-implicit scheme for the Euler equations of gas dynamics[END_REF][START_REF] Chalons | An all-regime Lagrange-projection like scheme for the gas dynamics equations on unstructured meshes[END_REF][START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equation[END_REF][START_REF] Cordier | An asymptotic-preserving all-speed scheme for the Euler and Navier-Stokes equations[END_REF][START_REF] Noelle | An asymptotic preserving all Mach number scheme for the Euler equations of gas dynamics[END_REF][START_REF] Dimarco | Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit[END_REF].

Klainerman and Majda have also shown that solutions of the compressible Euler equations converge to the solutions of the incompressible Euler equations as the Mach number tends to zero [START_REF] Klainerman | Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF]. The asymptotic preserving (AP) property is a consistency criterion for numerical schemes that have to deal with these two regimes. This property was firstly introduced by Jin in [START_REF] Jin | Efficient asymptotic-preserving (ap) schemes for some multiscale kinetic equations[END_REF]: a scheme is AP, if its lowest order multiscale expansion is a consistent discretization of the incompressible limit. Thus, an AP scheme for the Euler equations should provide a consistent and stable discretization independently of the Mach number. In the fully compressible regime, it should possess the desirable features of a compressible solver (e. g. good resolution of shocks and discontinuities and non-oscillatory solution profiles). Moreover, it should provide a consistent discretization of the incompressible equations when M → 0.

In the present work, our interest is the numerical simulation of flows with Mach numbers ranging from the incompressible limit to compressible regimes with M O(1). To this end, we have derived in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF] an implicit all-speed scheme, which is here extended for the accurate solution of multi-dimensional flows. Moreover, its AP property is proved analytically and numerically. The scheme has a general formulation that can be adopted without any structural modification for the simulation of waves propagating inside materials with different behavior, such as fluids and elastic solids. The ability of the scheme in dealing with very general equations of state (EOS) is achieved by adopting the Jin-Xin relaxation technique [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF]. With this method, the fluxes are relaxed at the continuous level and a linear transport operator is obtained:

this avoids a direct dependence of the spatial derivatives on the specific EOS. The second advantage consists in the fact that Riemann solvers are not necessary and fully implicit time integrators are easily implemented. This allows to get rid of the This manuscript is for review purposes only.

acoustic stability constraint (1.1), avoiding the requirement of a time stepping decreasing as M → 0, which would enormously increase the computational effort. The spatial discretization is built by a suitable combination of upwind and centered fluxes, in order to recover the correct viscosity at all speeds. The ability of the scheme in moderating the numerical viscosity with respect to standard explicit-upwind methods is here investigated by solving the Gresho vortex test and also different 2D Riemann problems involving material waves. Moreover, we propose an adaptive mesh refinement (AMR) algorithm, specifically designed for this all-speed scheme. We build a numerical estimate of the scheme entropy production and use it as an error indicator to drive the mesh refinement and coarsening. This estimate is also useful to assess the fact that the all-speed discretization is superior in approximating the different waves with respect to standard upwind-like schemes.

The structure of the paper is the following: in Section 2 we briefly revise the asymptotic analysis of the Euler equations for fluid dynamics. The 2D all-speed relaxation scheme is presented in Section 3. Then, the AP property is analytically proved in Section 4. The scheme is validated on fluid dynamics tests in Section 5,

where we also describe the adaptive mesh refinement algorithm based on the entropy residual for the all-speed scheme. Then, the scheme is extended to solve non-linear elasticity problems in Section 6, where two different limits occurring in elastic solids are analyzed. Conclusions are drawn in Section 7.

2. Asymptotic analysis of the continuous Euler equations. Let us consider the two dimensional Euler system:

(2.1)

     ∂ t ρ + ∇ x • (ρu) = 0 ∂ t (ρu) + ∇ x • (ρu ⊗ u + pI) = 0 ∂ t (ρe) + ∇ x • ((ρe + p) u) = 0,
ρ being the density, u the velocity field, p the pressure, I the identity matrix and e the total energy per unit mass, which is given by the sum of the kinetic energy and the internal energy per unit mass :

(2.2) e = 1 2 |u| 2 + .
For fluid dynamics problems, we consider the following state law, which extends the applicability of the ideal gas EOS:

(2.3) (ρ, s) = κ (s) γ -1 ρ γ-1 + p ∞ ρ = p + γp ∞ ρ (γ -1)
.

The internal energy is written as a function of the density and the entropy s and we have used the definition of the pressure for the second equality (see [START_REF] Whitham | Linear and nonlinear waves[END_REF][START_REF] Leveque | Numerical methods for conservation laws[END_REF]). Here,

γ = c p /c v is the polytropic gas constant, κ (s) = exp ((s -s 0 ) /c v ) with s 0 reference
entropy and p ∞ is a constant describing the intermolecular interaction, typical of stiffened gases, liquids and also some solids. The perfect gas EOS is recovered by setting p ∞ = 0. The general formulation for the speed of sound reads as follows

(2.4) c = ∂p ∂ρ s=const = γ (p + p ∞ ) ρ .
We now briefly revise the analysis of the low Mach number regime in fluid dynamics. For simplicity of notation, we consider a perfect gas, i. e. p ∞ = 0.

2.1. Non-dimensional Euler system. To write the non-dimensional Euler system, we follow Klainerman and Majda's works [START_REF] Klainerman | Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF][START_REF] Klainerman | Compressible and incompressible fluids[END_REF], by decomposing each quantity into a product of a reference value (denoted with the subscript * ) and a dimensionless number (denoted with • ), e.g. ρ = ρ * ρ for the density. The reference value is a "scaling" factor and it should be chosen in such a way that the dimensionless value is of order one. In this perspective, the reference pressure p * is defined with the reference sound speed c * = p * /ρ * (see definition (2.4)), which is also used in the scaling of the energy.

After some simple algebraic manipulations, most of the reference quantities cancel out, leading to the non-dimensional Euler system (we omit the hat notation • for the sake of simplicity):

       ∂ t ρ + div (ρu) = 0 (2.5) ∂ t (ρu) + div (ρu ⊗ u) + 1 M 2 ∇p = 0 (2.6) ∂ t (ρe) + div ((ρe + p) u) = 0. (2.7)
The non-dimensional Euler system depends only on a single non-dimensional reference quantity. This is the reference Mach number, defined as the ratio between the reference flow velocity and the reference sound speed as follows:

(2.8) M = u * c * .
The non-dimensional state law takes the following formulation:

(2.9)

p = (γ -1) ρe - M 2 2 ρ|u| 2 .

Low Mach number asymptotics.

In order to analyze the asymptotic behaviour of the Euler equations in the zero Mach number limit, we perform an expansion of the scaled variables in terms of the Mach number [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF]. The density is expanded as follows:

(2.10)

ρ = ρ 0 + M ρ 1 + M 2 ρ 2 + O M 3 ,
and all other variables are developed in the same way. The terms of zeroth order (subscript • 0 ) represent the zero Mach number limit. By inserting this expansion in the scaled system (2.5)-(2.6)-(2.7) and collecting terms with equal powers of M , one obtains:

• order O 1/M 2 :
(2.11) ∇p 0 = 0

• order O (1/M ):
(2.12)

∇p 1 = 0 • order O (1):      ∂ t ρ 0 + ∇ • (ρ 0 u 0 ) = 0 (2.13) ∂ t (ρ 0 u 0 ) + ∇ • (ρ 0 u 0 ⊗ u 0 ) + ∇p 2 = 0 (2.14) ∂ t (ρe) 0 + ∇ • ((ρe) 0 u 0 + p 0 u 0 ) = 0, (2.15)
with the first order of the state law (2.16)

p 0 = (γ -1) (ρe) 0 .
This manuscript is for review purposes only.

Relations (2.11) and (2.12) prove that the pressure is constant in space up to fluctuations of order M 2 . Hence, we can write the following pressure asymptotic:

(2.17) p (x, t) = P 0 (t) + M 2 p 2 (x, t) ,

where P 0 (t) is a thermodynamic pressure constant in space. In presence of open boundaries, the thermodynamic pressure P 0 is imposed to be equal to the exterior pressure P ext . For the sake of simplicity, we assume that the exterior pressure does not change in time, i. e. P 0 is constant in both space and time. From (2.16), we deduce that also the energy is constant in space and time and the incompressibility constraint ∇ • u 0 = 0 is easily derived.

Introducing this constraint into the continuity equation (2.13), the material derivative of the density is zero Dρ Dt = 0. This means that the density is constant along a trajectory of any fluid element. Therefore, when the incompressibility constraint is respected, the density of the fluid is constant in time and space, i. e. ρ 0 = const, in the case where the initial density of the fluid is constant in space.

With the study of the asymptotics carried out above, system (2.13)-(2.14)- (2.15) reduces to the incompressible Euler system in its non-dimensional form. This system is the zero Mach number limit of the compressible Euler system and reads as follows

(2.18)      ρ 0 = const ρ 0 (∂ t u 0 + (u 0 • ∇) u 0 ) + ∇p 2 = 0 ∇ • u 0 = 0.
3. Implicit relaxation all-speed scheme. To allow for an efficient and robust numerical procedure, we adopt the multi-dimensional Jin-Xin relaxation approach [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF]: a linear hyperbolic relaxation system is built to approximate the original system (here the Euler system (2.1)) with a small dissipative correction. Thanks to the linearity of the advection terms, the spatial derivatives lose their dependence on the state law. Numerical schemes are then derived in a general formulation that is not related to the EOS of a given material. Nonlinear terms appear only in the right hand side of the relaxation system, thus only diagonal terms are interested by the need of a linearization in the discretization process.

A fully implicit time integration of the relaxation system is proposed. Thus, acoustic CFL constraints are not required and a centered spatial discretization of the stiff parts when the M → 0 can be adopted, without stability issues. This is crucial to get the correct numerical viscosity in the low Mach number regime.

3.1. The Jin-Xin relaxation. Letting x = (x 1 , x 2 ) be the coordinates in the canonical basis of R 2 , u = (u 1 , u 2 ) the velocity components, system (2.1) may be rewritten in a general compact form in the following way:

(3.1) ∂ t ψ + ∂ x1 F (ψ) + ∂ x2 G (ψ) = 0,
where we have adopted the directional splitting of the flux function, with ψ ∈ R n and

F (ψ) , G (ψ) ∈ R n
. This is a system of n equations, with (x, t) ∈ R d , R + , with d = 2. In our case, we have n = 4 and the conservation variables and the fluxes along the two directions read as follows:

(3.2) ψ =     ρ ρu 1 ρu 2 ρe     , F (ψ) =     ρu 1 ρu 2 1 + p ρu 1 u 2 (ρe + p) u 1     , G (ψ) =     ρu 2 ρu 1 u 2 ρu 2 2 + p (ρe + p) u 2     .
This manuscript is for review purposes only.

The corresponding Jin-Xin relaxation system is constructed by introducing two vectors containing the relaxation variables v ∈ R n and w ∈ R n in the two directions, approximating the fluxes F (ψ) and G (ψ) respectively. The relaxation system of dimension n × (d + 1) takes the following formulation:

(3.3)          ∂ t ψ + ∂ x1 v + ∂ x2 w = 0 ∂ t v + A 1 ∂ x1 ψ = 1 η (F (ψ) -v) ∂ t w + A 2 ∂ x2 ψ = 1 η (G (ψ) -w) .
The small positive parameter η is called relaxation rate. The right hand sides of the second and third equations are stiff source terms and are the only non-linear parts of the system. The state law dependence is given by the fluxes in the sources, thus the EOS expression does not enter in the resulting linear transport operator.

The relaxation matrices A 1 and A 2 are positive diagonal matrices and are chosen by enforcing the subcharacteristic condition [START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF][START_REF] Whitham | Linear and nonlinear waves[END_REF]:

(3.4) A 1 -F (ψ) 2 ≥ 0 and A 2 -G (ψ) 2 ≥ 0 ∀ψ.
Specifically, we construct these matrices by a-priori estimating the wave speeds λ i of the original system (2.1). In the simplified case of A = A 1 = A 2 = diag{a i }, the eigenvalues of the relaxation system are computed as follows

(3.5) µ j = ± √ a i , i = 1, ..n, j = 1, ..2n,
thus having µ 1 ≥ max|u -c|, µ 2 = µ 3 ≥ max|u| and µ 4 ≥ max|u + c| for the Euler system (for more details on the construction of the relaxation matrices see [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF]). We remark that the relaxation is performed direction by direction, namely the equations on the relaxation variables are one dimensional problems.

The relaxation system (3.3) approximates the original system (3.1). This can be shown by applying the Chapman-Enskog expansion of the variables for small η [10]:

(3.6) v = v 0 + ηv 1 + O η 2 w = w 0 + ηw 1 + O η 2
At leading order O η (1) (zero relaxation limit η → 0), the original system is easily recovered, with the relaxation variables identically equal to the fluxes:

(3.7)      v 0 = F (ψ) w 0 = G (ψ) ∂ t ψ + ∂ x1 F (ψ) + ∂ x2 G (ψ) = 0.
The state satisfying (3.7) is called local equilibrium. Then, after some manipulations, the following first order approximation in expansion (3.6) is obtained

(3.8)          v 1 = -A 1 -F (ψ) 2 ∂ x1 ψ w 1 = -A 2 -G (ψ) 2 ∂ x2 ψ ∂ t ψ + ∂ x1 F (ψ) + ∂ x2 G (ψ) = η ∂ x1 A 1 -F (ψ) 2 ∂ x1 ψ + ∂ x2 A 2 -G (ψ) 2 ∂ x2 ψ ,
where F (ψ) and G (ψ) are the Jacobian matrices of the flux functions. It is evident that satisfying the subcharacteristic condition (3.4) amounts to ensuring the dissipative nature of system (3.8).

This manuscript is for review purposes only.

3.2. Numerical discretization. The derivation of the scheme is here extended to solve two-dimensional problems. System (3.3) is discretized with finite volumes on a Cartesian mesh. The first reason of this choice is due to the fact that the parallelization is very efficient on this kind of grids. Moreover, on a Cartesian grid, the direction-bydirection relaxation approximation (3.3) is easily discretized. As pointed out above, this produces one dimensional problems on the relaxation variables.

3.2.1. Implicit time discretization. We adopt a fully implicit integration in time with the aim of avoiding the acoustic CFL constraint (1.1). The resulting scheme possesses a general structure, with a time integrator that is not based on the specific EOS. Letting ∆t = t n+1 -t n be the time stepping, the implicit discretization at first order is a simple backward Euler scheme and reads:

(3.9)

               ψ n+1 -ψ n ∆t + ∂ x1 v n+1 + ∂ x2 w n+1 = 0 v n+1 -v n ∆t + A 1 ∂ x1 ψ n+1 = 1 η F ψ n+1 -v n+1 w n+1 -w n ∆t + A 2 ∂ x2 ψ n+1 = 1 η G ψ n+1 -w n+1
The non-linear fluxes in the right hand sides are solved with one iteration of the Newton's method, which is enough to obtain convergence for the proposed first order scheme. The resulting approximation consists in a Taylor expansion:

(3.10) F ψ n+1 F (ψ n ) + F (ψ n ) ψ n+1 -ψ n
and the same holds for G (ψ). F (ψ n ) and G (ψ n ) are the Jacobians of the fluxes in the two directions and can be computed analytically. Linearization (3.10) introduces a coupling among all the equations and the following linear system is obtained:

(3.11)      LΨ n+1 + MV n+1 + NW n+1 = r PΨ n+1 + QV n+1 = s 1 TΨ n+1 + UW n+1 = s 2
where Ψ n+1 , V n+1 and W n+1 are the vectors containing the grid point values of the conservative and relaxation variables and the matrices structure is given by the spatial discretization (see Sec. 3.2.2). Relation (3.10) implies that the Jacobians enter only on the diagonal terms of the sub-matrices of P and T. The full linear system is solved with the GMRES iterative solver implemented in the PETSc library [START_REF] Balay | PETSc users manual[END_REF].

Scheme (3.9) is fully implicit, thus unconditionally stable (see [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF]). Nevertheless, ∆t has to be chosen accordingly to the desired accuracy. If the aim is the accurate resolution of material waves, the enforcement a material CFL condition is recommended. The material Courant number is defined on the speed of the material wave as ν mat = µ mat ∆t/∆x, where we use the "material" eigenvalue of the relaxation system µ mat = √ a mat ≥ max|u|, with µ j defined in (3.5). A material CFL condition is enforced by setting ν mat ≤ 1. This way the problem of a Mach-dependent ∆t shown by relation (1.1) is avoided. On the other hand, since the material CFL does not depend on the speed of the fast waves, these waves are not captured in the low Mach limit. If an accurate resolution also of the acoustic waves is required, a standard acoustic CFL has to be enforced, adapting relation (1.1) to the relaxation system.

This gives ν ac = µ max ∆t/∆x, where µ max = √ a max ≥ max|u + c|.

3.2.2.

All-speed spatial discretization. Let ∆x l be the grid spacing in the

x l direction and Ω ij the control volume centered in the node (i∆x 1 , j∆x 2 ). For a generic variable h, h ij denotes the approximate cell average of the variable in the cell

C ij = x i-1/2,j , x i+1/2,j × x i,j-1/2 , x i,j+1/2
and h i+1/2,j denotes the approximate point value at the cell interface x i+1/2,j . The spatial discretization of (3.3) reads (3.12)

             ∂ t ψ ij + v i+1/2,j -v i-1/2,j ∆x 1 + w i,j+1/2 -w i,j-1/2 ∆x 2 = 0 ∂ t v ij + A 1 ψ i+1/2,j -ψ i-1/2,j ∆x 1 = 1 η F ψ ij -v ij ∂ t w ij + A 2 ψ i,j+1/2 -ψ i,j-1/2 ∆x 2 = 1 η G ψ ij -w ij .
The numerical fluxes in (3.12) are computed by constructing a convex combination of upwind and centered schemes. This is based on the local Mach number of the flow M loc . For a generic variable h this yields

(3.13) h i+1/2,j = f (M loc ) h i+1/2,j upw + (1 -f (M loc )) h i+1/2,j cent .
The same holds for direction x 2 . The function f (M loc ) has to satisfy 0 ≤ f (M loc ) ≤ 1 (in the numerical experiments we choose f (M loc ) = min{1, M loc }) and M loc is computed at the previous time step at the numerical interface x i+1/2,j . This allows to recover the correct numerical viscosity for each considered regime [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF][START_REF] Abbate | A relaxation scheme for the simulation of low Mach number flows[END_REF].

The numerical fluxes with the centered scheme are computed as follows:

(3.14) h i+1/2,j cent = 1 2 (h i+1,j + h ij ) .
This discretization is stable inside the adopted fully implicit framework. Since we are centering all conservative and relaxation variables, the pressure gradients result to be discretized with central differencing, thus providing the correct Mach scaling [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF]. The upwind scheme is built along the characteristic variables of the relaxation system, as in [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF], getting the following numerical fluxes along direction x 1 :

(3.15)

   ψ i+1/2,j upw = 1 2 ψ i+1,j + ψ ij -1 2 A -1/2 1 (v i+1,j -v ij ) v i+1/2,j upw = 1 2 (v i+1,j + v ij ) -1 2 A 1/2 1 ψ i+1,j -ψ ij .
The full 2D implicit all-speed scheme then reads:

(3.16)                                          ψ n+1 ij -ψ n ij ∆t + 1 2∆x 1 v n+1 i+1,j -v n+1 i-1,j - f (M loc ) A 1/2 1 2∆x 1 ψ n+1 i+1,j -2ψ n+1 ij + ψ n+1 i-1,j + 1 2∆x 2 w n+1 i,j+1 -w n+1 i,j-1 - f (M loc ) A 1/2 2 2∆x 2 ψ n+1 i,j+1 -2ψ n+1 ij + ψ n+1 i,j-1 = 0 v n+1 ij -v n ij ∆t + A 1 2∆x 1 ψ n+1 i+1,j -ψ n+1 i-1,j - f (M loc ) A 1 1/2 2∆x 1 v n+1 i+1,j -2v n+1 ij + v n+1 i-1,j = 1 η F ψ n+1 ij -v n+1 ij w n+1 ij -w n ij ∆t + A 2 2∆x 2 ψ n+1 i,j+1 -ψ n+1 j,i-1 - f (M loc ) A 1/2 2 2∆x 2 w n+1 i,j+1 -2w n+1 ij + w n+1 i,j-1 = 1 η G ψ n+1 ij -w n+1 ij .
This manuscript is for review purposes only.

The asymptotic preserving property. The Asymptotic-Preserving (AP)

property is defined considering a continuous physical model S M (in our case the compressible Euler system (2.1)) that involves a perturbation parameter M (the acoustic Mach number). The perturbation parameter can range from M 1 to M 1 values.

In the Euler case, there exists a reduced system S 0 , which is the limit system of S M as M → 0, i. e. S 0 is the incompressible system (2.18). Then, let S ∆ M be a numerical scheme providing a consistent discretization of S M , with discrete time and space steps ∆ = (∆t, ∆x). The scheme S ∆ M is said to be AP if the two following properties are verified [START_REF] Jin | Efficient asymptotic-preserving (ap) schemes for some multiscale kinetic equations[END_REF]:

1. its stability condition is independent of M , namely the time step ∆t does not depend on the Mach number of the flow;

2. as M goes to zero, there exists the limit discrete S ∆ 0 , which provides a consistent discretization of the continuous limit system S 0 .

Here we show that scheme (3.16) is AP. Property 1 is satisfied, since the scheme is fully implicit, thus unconditionally stable. In order to prove that property 2 is respected, we write the limit discrete scheme S ∆ 0 as M → 0 of the implicit relaxation scheme and show that it is consistent with the continuous limit model S 0 .

In what follows, we begin with the non-dimensionalization of the scheme and then we carry out the analysis of its asymptotics. We expose the reasoning on the time semi-discrete scheme (3.9) for readability. The extension to the full time and space discretization is straightforward.

4.1. Non-dimensional implicit relaxation scheme. In deriving the nondimensional implicit relaxation scheme, we adopt the same notation of Section 2.1.

Here, we also have to scale the relaxation variables v and w: due to the relaxation leading order (3.7), these variables have the same physical dimensions of the fluxes F (ψ) and G (ψ) defined in (3.2). By considering the pressure as predominant, we choose to scale v 2 and w 3 using the speed of sound. Thus, we obtain the following non-dimensional formulation of the semi-discrete scheme (3.9) (we omit from now on the subscript • for simplicity of notation):

1. the non-dimensional conservation of mass, setting

z 1 = [v 1 , w 1 ] T , reads (4.1)
     ρ n+1 -ρ n + ∆t∇ • z n+1 1 = 0 z n+1 1 -z n 1 + a 1 ∆t∇ρ n+1 = ∆t η (ρu) n+1 -z n+1 1 .
2. the non-dimensional conservation of momentum is given by the two parts:

(4.2)                (ρu 1 ) n+1 -(ρu 1 ) n + ∆t ∂ x1 v n+1 2 M 2 + ∂ x2 w n+1 2 = 0 v n+1 2 -v n 2 M 2 + a 2 ∆t∂ x1 (ρu 1 ) n+1 = ∆t η ρu 2 1 n+1 + p n+1 M 2 - v n+1 2 M 2 w n+1 2 -w n 2 + a 2 ∆t∂ x2 (ρu 1 ) n+1 = ∆t η (ρu 1 u 2 ) n+1 -w n+1 2
This manuscript is for review purposes only.

(4.3)                (ρu 2 ) n+1 -(ρu 2 ) n + ∆t ∂ x1 v n+1 3 + ∂ x2 w n+1 3 M 2 = 0 v n+1 3 -v n 3 + a 3 ∆t∂ x1 (ρu 2 ) n+1 = ∆t η (ρu 1 u 2 ) n+1 -v n+1 3 w n+1 3 -w n 3 M 2 + a 3 ∆t∂ x2 (ρu 2 ) n+1 = ∆t η ρu 2 2 n+1 + p n+1 M 2 - w n+1 3 M 2 . 3. setting z 4 = [v 4 , w 4 ]
T , the non-dimensional conservation of energy is given by

(4.4)      (ρe) n+1 -(ρe) n + ∆t∇ • z n+1 4 = 0 z n+1 4 -z n 4 + a 4 ∆t∇ (ρe) n+1 = ∆t η (ρe) n+1 + p n+1 u n+1 -z n+1 4 ,
with the scaled state law (4.5)

p n+1 = (γ -1) (ρe) n+1 - M 2 2 ρ n+1 |u| n+1 .
4.2. Asymptotics of the implicit relaxation scheme. In the spirit of studying the low Mach number asymptotics, we develop all scaled variables, i. e. both conservative and relaxation variables, in powers of the Mach number, as done in Eq.

(2.10). The expansion of the scaled relaxation variable v 1 reads as follows

(4.6) (v 1 ) n+1 = v 0,η 1 n+1 + M v 1,η 1 n+1 + M 2 v 2,η 1 n+1 + O M 3 ,
where we have introduced the notation with two superscripts, the first one indicating the order in the power of M and the second one indicating the order in the power of η.

Here, we are keeping η = η fixed. Terms of zero-th order (superscript • 0,η=0 ) represent the zero Mach number limit in the zero relaxation limit. Relaxation variables may be expanded also in powers of η, as in a Chapman-Enskog expansion (3.6). One could then combine expansions (4.6) and (3.6) and write a full expansion for v and w in powers of both M and η. The expansion for v 1 reads

(v 1 ) n+1 = v 0,0 1 n+1 + η v 0,1 1 n+1 + M v 1,0 1 n+1 + η v 1,1 1 n+1 + M 2 v 2,0 1 n+1 + η v 2,1 1 n+1 + O M 3 + O η 2 .
(4.7)

By setting M = 0 and η = 0, the incompressible Euler system can be recovered.

Since we need to preserve the low Mach number behaviour in the expansions introduced above, we require that η M . This is necessary to recover the correct This manuscript is for review purposes only.

1. order O M 1/M 2 : (4.8)      ∂ x1 v 0,η 2 n+1 = 0 v 0,η 2 n+1 -v 0,η 2 n = ∆t η p n+1 0 -v 0,η 2
n+1
and the same holds for w 0,η 3 . At this point, we substitute expansion (3.6) in powers of η in (4.8) (once again, the same holds for variable w 3 ):

       ∂ x v 0,0 2 n+1 + η v 0,1 2 n+1 = 0 v 0,0 2 n+1 + η v 0,1 2 n+1 -v 0,0 2 n -η v 0,1 2 n = ∆t η p n+1 0 -v 0,0 2 n+1 -η v 0,1 2
n+1
.

We are interested in the zero relaxation limit η → 0, hence we collect the terms

O η (1/η) in the last equation, obtaining (v 0,0 2 ) n+1 = p n+1 0 and (w 0,0 3 ) n+1 = p n+1 0
. This is plugged in (4.8), obtaining:

   v 0,0 2 n+1 = p n+1 0 ∂ x1 p n+1 0 = 0    w 0,0 3 n+1 = p n+1 0 ∂ x2 p n+1 0 = 0. It is clear then that ∇p n+1 0 = 0 is respected. 2. order O M (1/M ): (4.9)      ∂ x1 v 1,η 2 n+1 = 0 v 1,η 2 n+1 -v 1,η 2 n = ∆t η p n+1 1 -v 1,η 2
n+1
and the same holds for w 1,η 3 . After inserting the expansion in powers of η, taking the zero relaxation limit η → 0 and collecting the terms O η (1/η), we obtain the following relations:

   v 1,0 2 n+1 = p n+1 1 ∂ x1 p n+1 1 = 0    w 1,0 3 n+1 = p n+1 1 ∂ x2 p n+1 1 = 0.
This means that also ∇p n+1 1 = 0 is respected.

order O M (1):

• for the conservation of mass we have:

(4.10)        ρ n+1 0 -ρ n 0 + ∆t∇ • z 0,η 1 n+1 = 0 z 0,η 1 n+1 -z 0,η 1 n + a 1 ∆t∇ρ n+1 0 = ∆t η (ρu) n+1 0 -z 0,η 1 n+1
.

Once again, we expand in powers of η and take the zero relaxation limit by collecting terms of order O η (1/η), obtaining (4.11)

   z 0,0 1 n+1 = (ρu) n+1 0 ρ n+1 0 -ρ n 0 + ∆t∇ • (ρu) n+1 0 = 0.
This manuscript is for review purposes only.

• for the conservation of momentum we have:

(4.12)                (ρu 1 ) n+1 0 -(ρu 1 ) n 0 + ∆t ∂ x1 v 2,η 2 n+1 + ∂ x2 w 0,η 2 n+1 = 0 v 2,η 2 n+1 -v 2,η 2 n + a 2 ∆t∂ x1 (ρu 1 ) n+1 0 = ∆t η ρu 2 1 n+1 0 + p n+1 2 -v 2,η 2 n+1 w 0,η 2 n+1 -w 0,η 2 n + a 2 ∆t∂ x2 (ρu 1 ) n+1 0 = ∆t η (ρu 1 u 2 ) n+1 0 -w 0,η 2 n+1 , (4.13)                (ρu 2 ) n+1 0 -(ρu 2 ) n 0 + ∆t ∂ x1 v 0,η 3 n+1 + ∂ x2 w 2,η 3 n+1 = 0 v 0,η 3 n+1 -v 0,η 3 n + a 3 ∆t∂ x1 (ρu 2 ) n+1 0 = ∆t η (ρu 1 u 2 ) n+1 0 -v 0,η 3 n+1 w 2,η 3 n+1 + w 2,η 3 n + a 3 ∆t∂ x2 (ρu 2 ) n+1 0 = ∆t η ρu 2 2 n+1 0 + p n+1 2 -w 2,η 3 n+1
In the expansion in powers of η, terms of order O η (1/η) give expressions for (v 2,0 2 ) n+1 , (v 0,0 3 ) n+1 , (w 0,0 2 ) n+1 and (w 2,0 3 ) n+1 . This yields

(4.14)                      v 2,0 2 n+1 = ρu 2 1 n+1 0 + p n+1 2 w 0,0 2 n+1 = (ρu 1 u 2 ) n+1 0 v 0,0 3 n+1 = (ρu 1 u 2 ) n+1 0 , w 2,0 3 n+1 = ρu 2 2 n+1 0 + p n+1 2 (ρu) n+1 0 -(ρu) n+1 0 + ∆t∇ • (ρ 0 u 0 ⊗ u 0 ) n+1 + ∆t∇p n+1 2 = 0.
• for the conservation of energy we adopt the same procedure of the conservation of mass, getting

+ p n+1 0 ∇ • u n+1 0 = 0.
This manuscript is for review purposes only.

We now assume that the boundary conditions are such that p n+1 the limit M → 0. This implies that a fully centered discretization is adopted in the incompressible limit, i. e. the numerical viscosity is consistent with the incompressible regime.

5. Numerical validations on fluid dynamics problems. At first, the Gresho vortex test case is analyzed, in order to verify that the scheme possesses the correct numerical viscosity in the low Mach number regime. This test verifies that the proposed scheme is able to preserve the incompressible regime at all times, after setting an incompressible initial flow. The vortex is solved for both perfect and stiffened gases. Then, we introduce an estimate of the numerical entropy production of the scheme. We use this indicator to show that the all-speed discretization in more precise with respect to an upwind one and also to perform grid adaptivity. Due to the AP estimates derived in Sec. 4.2, we set η = 10 -8 in all the proposed computational experiments.

5.1. The Gresho vortex. We test the relaxation all-speed scheme (3.16) on the classical Gresho vortex test case [START_REF] Gresho | On the theory of semi-implicit projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix. part 1: Theory[END_REF][START_REF] Liska | Comparison of several difference schemes on 1D and 2D test problems for the Euler equations[END_REF]. This vortex is a stationary solution of the incompressible Euler equations, where centrifugal forces are exactly balanced by pressure gradients. A rotating vortex is positioned at the center (0.5, 0.5) of the computational domain [0, 1] × [0, 1]. The initial conditions are specified in terms of the radial distance r = (x -0.5) 2 + (y -0.5) 2 in the form

     ρ (x, y, 0) = 1 u 1 (x, y, 0) = -u φ (r) sinφ u 2 (x, y, 0) = u φ (r) cosφ
The rotation is initiated by imposing a simple angular velocity distribution of

u φ (x, y, 0) =      5r 0 ≤ r ≤ 0.2 2 -5r 0.2 ≤ r ≤ 0.4 0 r ≥ 0.4, p (x, y, 0) =      p 0 + 25 2 r 2 0 ≤ r ≤ 0.2 p 0 + 25 2 r 2 + 4 (1 -5r -ln 0.2 + ln r) 0.2 ≤ r ≤ 0.4 p 0 -2 + 4 ln 2 r ≥ 0.4.
and the background pressure is adjusted such that it matches the maximum Mach number M max :

(5.1)

p 0 = ρ (u φ ) 2 max (γM 2 max ) = ρ (γM 2 max )
since the maximum velocity (u φ ) max = 1, which is reached in r = 0.2

We follow the flow over one full rotation of the vortex with different M max in a perfect gas, which is completed at time t = 1. We compare the results obtained with the relaxation all-speed scheme (3.16) with the results of a standard explicit-upwind discretization of the relaxation system. All the Gresho vortex tests are performed on a uniform grid of 128×128 cells and a material CFL constraint of ν mat = 0.2. In test 1

we set M max = 10 -1 : the initial condition in t = 0 is presented in Fig. 1(a). In Figs.

This manuscript is for review purposes only. 1(b) and 1(c), we compare the results of the relaxation all-speed scheme (3.16) and of the explicit-upwind scheme of [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF]. A comparison with the initial Mach number distribution shows that the all-speed spatial discretization accurately preserves the shape of the vortex, due to the convex combination of upwind and centered fluxes.

Instead, an upwind scheme is too diffusive for the targeted regime and thus the shape of the vortex is lost. We also plot the results obtained after 2 and 3 rotations of the vortex in Fig. 2, confirming that the all-speed scheme is able to recover the correct weakly incompressible solution also after long times.

For tests 2 and 3 we set M max = 10 -2 and M max = 10 -3 respectively. The results after one full rotation are reported in Fig. 3, for the relaxation all-speed scheme and the explicit-upwind scheme [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF]. The shape of the two vortices is completely diffused when adopting an upwind flux discretization, whereas the all-speed convex combination accurately preserves the initial vortex shape for both cases, besides a small noise probably due to directional splitting.

In Table 1, we report the total kinetic energy in the simulation domain at time t = 1 relative to the total kinetic energy at time t = 0 for the three tests. With the all-speed scheme, the kinetic energy reduces by about 1.5 per cent over one rotation of the vortex. However, this loss is clearly independent of the Mach number of the flow.

On the contrary, when adopting conventional upwind discretizations the dissipation rate of kinetic energy consistently increases as the Mach number decreases.

We also perform a study on the pressure fluctuations p f l = (p max -p min ) /p max , computed on the same grid at time t = 1 for the three considered Mach numbers.

This manuscript is for review purposes only. These results are reported in Table 2 and they show that pressure fluctuations scale exactly with M 2 when adopting the all-speed convex combination. We can thus infer that the simulated flow is kept in the incompressible regime by the proposed numerical scheme. This is clearly not the case with an upwind spatial discretization. The scaling of p f l as M , theoretically expected for this latter scheme, can be asymptotically reached for small Mach numbers (it can be seen by comparing the simulations with M max = 10 -2 and M max = 10 -3 ). We here adapt the standard Gresho vortex test case to a water flow. It suffices This manuscript is for review purposes only. to adjust the background pressure (5.1) with the stiffened gas state law as follows:

M max = 10 -1 M max = 10 -2 M max = 10
p 0 = ρ (u φ ) 2 max (γM 2 max ) -p ∞ .
The initial density is set to ρ = 1000Kg/m 3 and the water parameters γ = 4.4 and p ∞ = 6.8 • 10 8 are given with state law (2.3). The results for M max = 10 -2 obtained with the two schemes are reported in Figs. 4 at time t = 1. The all-speed scheme is preserving the vortex shape also in the case of water.

Numerical entropy production.

The behaviour of the all-speed discretization on the different waves can be assessed by studying an a-posteriori error indicator. We adopt the numerical entropy production introduced in [START_REF] Puppo | Numerical entropy production for central schemes[END_REF][START_REF] Puppo | Numerical entropy and adaptivity for finite volume schemes[END_REF], since the entropy is naturally available for any system of conservation laws with an entropy inequality and has a well defined physical meaning. Moreover, the entropy production of a scheme scales as the truncation error in the regular regions and its behavior allows to distinguish between contact discontinuities and shocks. In [START_REF] Puppo | Numerical entropy and adaptivity for finite volume schemes[END_REF] Let us consider the entropy pair (η, ζ). The entropy inequality ∂ t η +∇ x •ζ ≤ 0 has to be integrated with the finite volume scheme used to integrate the hyperbolic system that we are interested in solving. Coherently with discretizations (3.9)-(3.12), we employ a first order implicit scheme for time integration in a finite volume framework.

We get the following numerical entropy production in every cell:

(5.2)

S n+1 ij = η ψ n+1 ij -η ψ n ij + ∆t ∆x 1 ζ n+1 1;i+1/2,j -ζ n+1 1;i-1/2,j + ∆t ∆x 2 ζ n+1 2;i,j+1/2 -ζ n+1 2;i,j-1/2 ,
where we are considering the spatial discretization of the two components of the nu-

merical entropy flux ζ = (ζ 1 , ζ 2 )
. The interface values have to be computed according to the all-speed convex combination (3.13). This is equivalent to a Lax-Friedrichs scheme with a numerical viscosity modulated by the local Mach number, namely:

(5.3)

ζ n+1 1;i+1/2,j = 1 2 ζ n+1 1;i+1,j + ζ n+1 1;i,j - √ a max f (M loc ) 2 η n+1 i+1,j -η n+1 i,j
.

In the numerical experiments, we build the entropy pair on the physical entropy of This manuscript is for review purposes only.

the Euler system, as follows:

η (ψ) = -ρ log p + p ∞ ρ γ , ζ (ψ) = η (ψ) u.
The proposed estimate can also be adopted as a criterion that is able to pivot mesh adaptivity. An effective adaptive algorithm has to be driven by an indicator able to provide a robust a-posteriori measure of the local error and also to recognize the qualitative structure of the flow. This will be the focus of Section 5.2.2.

2D Riemann problem.

We use the entropy production to compare the behaviour of the all-speed and of the upwind discretizations on contact waves. We study a Riemann problem in the squared domain [0, 1] × [0, 1], with the following initial conditions [START_REF] Kurganov | Solution of two-dimensional riemann problems for gas dynamics without riemann problem solvers[END_REF][START_REF] Schulz-Rinne | Numerical solution of the riemann problem for two-dimensional gas dynamics[END_REF]: We adopt a grid of 256 × 256 cells and we enforce an acoustic constraint on the time step ν ac = 0.9, in order to have stability for the explicit-upwind scheme and in order to have a good resolution of all the propagating waves. In Fig. 5 we compare the density contours and the numerical entropy production for the two schemes. It is evident from the density contours that the all-speed property helps in keeping sharp the two contact waves and the shock. The contact wave at the bottom of the domain is in the low Mach number regime, since M 8 • 10 -3 . Also the small vortex in the center is more accurately captured by the all-speed scheme. The entropy plots confirm that the all-speed discretization is superior in capturing the solution structure. For both schemes, the entropy production is higher in the shock wave region, as expected by the theory of [START_REF] Puppo | Numerical entropy and adaptivity for finite volume schemes[END_REF]. The all-speed discretization produces more entropy than the explicit-upwind one, confirming that this region is more accurately resolved by the proposed scheme.

(p, ρ, u 1 , u 2 ) (x 1 , x 2 , t = 0) =          (1, 1, 0, -0.4) , Ω 1 (1, 2, 0, -0.3) , Ω 2 (0.

Adaptive mesh refinement.

As anticipated above, the numerical entropy production can be exploited to decide where to locally refine or coarsen the mesh. In order to reduce the computational costs, the discretization of the solution can be done by using non-conforming hierarchical meshes. Specifically, we use the adpative mesh refinement (AMR) techinque [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF].

Fig. 6. Decomposition of a squared domain and corresponding quadtree [START_REF] Raeli | A finite-difference method for the variable coefficient poisson equation on hierarchical cartesian meshes[END_REF].

Octree-meshes allow for a strong reduction of the number of degrees of freedom where the problem exhibits smooth behavior and also a strongly localized increase of information in areas needing more accuracy. Octrees are a hierarchical data structure based on the principle of recursive decomposition of space. Each internal node has exactly four children (quadtree) for 2D problems, and eight children (octree) for 3D

problems. Here we focus on quadtree meshes, that are non-conforming hierarchical meshes defined in a square as shown in Fig. 6. We use the Bitpit library [START_REF] Engineering | Bitpit Web page[END_REF] for the efficient implementation of our computational grid. The hierarchical nature of the grid makes mesh generation, adaptivity and partitioning very efficient and with a low-memory footprint. The data structure is based on a linear quadtree [START_REF] Frisken | Simple and efficient traversal methods for quadtrees and octrees[END_REF], namely only the leafs of the tree structure are stored. This structure is easily dispatched to a distributed memory architecture and parallel communications are limited to only the first layer of neighboring cells. This constraint is perfectly in line with the stencil of the proposed numerical scheme (3.16). A Z-order index is assigned to every cell [START_REF] Morton | A computer oriented geodetic data base and a new technique in file sequencing[END_REF].

More details on the library can be found in [START_REF] Engineering | Bitpit Web page[END_REF][START_REF] Raeli | A finite-difference method for the variable coefficient poisson equation on hierarchical cartesian meshes[END_REF].

We design an AMR algorithm based on the entropy production (5.2). We start This manuscript is for review purposes only. from a uniform coarse grid of 2 d × 2 d grid points. This grid is associated to the minimum level L min = d of the quadtree data structure. At every refinement, each cell of the grid may be replaced by four children. Let L max be the maximum refinement level allowed for a grid. At the end of every time step, the following procedure is implemented:

1. the quantity S n ij is computed in every cell with (5.2); This manuscript is for review purposes only. We test the AMR algorithm on the previous 2D Riemann problem (Section 5.2.1).

2. if |S n ij | > S ref
We set L min = 5, S ref = 0.002 and S coa = 0.0001 and we enforce a material CFL constraint ν mat = 0.3. L max is varied from 6 to 9 and these grids at time t = 0.3 are reported in Fig. 7. We observe that the entropy criterion is pivoting the AMR in the correct way, since the refinement levels are introduced where the 4 waves occur.

By comparing the grids of Fig. 7 and the solution structure in Fig. 8, we can see that the postshock region is refined also after the solution has become smooth. For this reason, the refinement needs to be further optimized. In Fig. 8 we compare the density profile obtained with a uniform grid built with L min = L max = 10, which gives 1,048,576 cells and the density obtained with the AMR by setting L min = 5 and L max = 10, which gives a total number of cells of 146,578 at the end of the simulation. The solution structure, the different waves and the small vortex in the center are accurately approximated for both grids, gaining a very similar precision (the small oscillations that can be seen on the contours of Fig. 8(b) are due to the visualization software that is not able to fully handle quadtree grids). However, it is evident that with the AMR technique the computational effort is reduced, since the number of degrees of freedom is consistently reduced (of about 10 times). The CPU times required with the uniform grid and with the AMR grid, using 32 cores, are respectively 9.25 • 10 3 and 1.61 • 10 3 seconds.

6. An extension to non-linear elasticity. As anticipated above, the relaxation all-speed scheme (3.16) can be adopted also for the simulation of flows with more complex EOS. Specifically, the accurate approximation of the deformation of elastic solids can be addressed. To this end, we adopt a monolithic Eulerian model describing different materials with the same system of conservation laws. The model derivation has been exstensively discussed in several previous works [START_REF] Godunov | Elements of continuum mechanics and conservation laws[END_REF][START_REF] Plohr | A conservative Eulerian formulation of the equations for elastic flow[END_REF][START_REF] Gorsse | A simple Cartesian scheme for compressible multimaterials[END_REF][START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF].

Here we briefly report the 2D system and the chosen hyperelastic state law.

Let Ω 0 ∈ R 2 the initial (or reference) configuration of a continuum and Ω t ∈ R 2 the deformed configuration at time t. The considered system of conservation laws models fluids and hyperelastic solids in the Eulerian framework, thus it is written in the deformed configuration Ω t . As for fluids, the equations describing the evolution of elastic solids are the conservation of mass, of momentum and energy. In addition, a law for the description of the deformation is needed: this is done with the introduction This manuscript is for review purposes only.

of equations of transport of the backward characteristics. These functions describe the continuum in the Eulerian framework: for a time t and a point x in the deformed configuration, the corresponding initial point ξ in the initial configuration is given, i. e. Y : Ω t × [0, T ] → Ω 0 , (x, t) → Y (x, t). One can also introduce the forward characteristics X (ξ, t), defined as the image at time t in the deformed configuration of a material point ξ belonging to the initial configuration, i. e. X : Ω 0 × [0, T ] → Ω t , (ξ, t) → X (ξ, t). The corresponding Eulerian velocity field u is defined as u :

Ω t × [0, T ] → R 2 , (x, t) → u (x, t).
∂ t Y + u • ∇Y = 0 Y (x, 0) = x, x ∈ Ω t .
Since the stress tensors have a direct dependence on [∇Y ], the gradient of (6.1) is taken as a governing equation (see for details [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF][START_REF] Brauer | A Cartesian scheme for compressible multimaterial models in 3d[END_REF]). The full monolithic Eulerian model thus is written as follows:

(6.2)

         ∂ t ρ + ∇ x • (ρu) = 0 ∂ t (ρu) + ∇ x • (ρu ⊗ u -σ) = 0 ∂ t ([∇ x Y ]) + ∇ x (u • [∇ x Y ]) = 0 ∂ t (ρe) + ∇ x • ρeu -σ T u = 0.
Here σ is the Cauchy stress tensor, which is derived through the chosen EOS. The total energy e is still given by expression (2.2), but now the internal energy has to account for the behaviour of different materials, including gases, liquids and solids.

Therefore, we adopt the following general constitutive law [START_REF] Holzapfel | Nonlinear solid mechanics: a continuum approach for engineering science[END_REF][START_REF] Gorsse | A simple Cartesian scheme for compressible multimaterials[END_REF]:

(6.3) (ρ, s, [∇ x Y ]) = κ (s) γ -1 ρ γ-1 + p ∞ ρ general gas + χ ρ trB -2 neohookean solid
, where B is the normalized Cauchy stress tensor:

B = B J = [∇ x Y ] -1 [∇ x Y ] -T J , J = det ([∇ x Y ]) -1 .
The energy function (6.3) The general expression of the Cauchy stress tensor σ is easily obtained [START_REF] Brauer | A Cartesian scheme for compressible multimaterial models in 3d[END_REF][START_REF] Abbate | Numerical methods for all-speed flows in fluid-dynamics and non-linear elasticity[END_REF] (6.4)

   σ (ρ, s, [∇ x Y ]) = -p (ρ, s) I + 2χJ -1 B - trB 2 I p (ρ, s) = -p ∞ + k (s) ρ γ .
We rewrite the Eulerian model (6.2) in the compact formulation with the directional splitting (3.1), having ψ, F (ψ) , G (ψ) ∈ R n with n = 8: (6.5)

ψ =             ρ ρu 1 ρu 2 Y 1 ,1 Y 2 ,1 Y 1 ,2 Y 2 ,2 ρe             , F (ψ) =             ρu 1 ρu 2 1 -σ 11 ρu 1 u 2 -σ 21 u 1 Y 1 ,1 + u 2 Y 1 ,2 u 1 Y 2 ,1 + u 2 Y 2 ,2 0 0 ρe -σ 11 u 1 -σ 21 u 2             , G (ψ) =             ρu 2 ρu 1 u 2 -σ 12 ρu 2 2 -σ 22 0 0 u 1 Y 1 ,1 + u 2 Y 1 ,2 u 1 Y 2 ,1 + u 2 Y 2 ,2 ρe -σ 22 u 2 -σ 12 u 2             .
In this notation, the superscript i indicates the component of Y and the subscript ,j stands for the direction along which the derivative is calculated. With this formulation, the Jin-Xin relaxation system (3.3) is easily constructed, by introducing the relaxation variables vectors. Then, the implicit time integration and the all-speed spatial discretization are adopted just as in the case of fluid dynamics, with the analytical derivation of the jacobians of the fluxes F (ψ) and G (ψ) to be used inside linearization (3.10) (for the jacobians derivation see Appendix A). Thus, the scheme has exactly the same formulation derived for the Euler equation case: this means that the relaxation all-speed scheme introduced in Section 3 can solve flows with more complex EOS, including elastic deformations, without any substantial modification.

6.1. Low Mach number limits in elastic solids. We now analyze the low Mach number regime in compressible elastic solids. This regime can be seen as a small deformation or a deformation that is slow with respect to other propagating waves. Alongside with the sound speed, which is computed as for fluid dynamics with expression (2.4), we define an "elastic speed " in the following way: (6.6)

u el = 2χ ρ .
We can thus define an "isochoric Mach number" on this elastic speed, as the ratio between the isochoric speed and the advective velocity [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF]:

(6.7) M χ = u u el = ρu 2 2χ .
Therefore, two different scales can be distinguished. They can be identified when the Eulerian system is non-dimensionalized. As done for the standard case of the Euler system in Section 2.1, every variable is decomposed into a product of a reference value and a dimensionless number. By using the definitions of the two Mach numbers (2.8) and (6.7) and by using the speed of sound to scale pressure and total energy, we get the following formulation of the non-dimensional Eulerian system (for simplicity of This manuscript is for review purposes only.

notation we loose the hat • on the non-dimensional variables):

(6.8)                  ∂ t ρ + ∇ x • (ρu) = 0 ∂ t (ρu) + ∇ x • (ρu ⊗ u) + ∇ x p M 2 - 2χ M 2 χ ∇ x J -1 B - trB 2 I = 0 ∂ t ([∇ x Y ]) + ∇ x (u • [∇ x Y ]) = 0 ∂ t ρe M 2 + ∇ x • ρe M 2 u + p M 2 - 2χ M 2 χ ∇ x J -1 B - trB 2 I u = 0.
Here we have substituted the expression for the stress tensor (6. number consistently smaller with respect to the elastic Mach number. This limit can be observed in rubber-like materials. In this regime, the gradient of stress tensor σ can be split into two different scales: the acoustic scale, which is represented by the pressure gradient and the elastic one, which is represented by the gradient of the elastic deformation. The pressure gradient is the only responsible for the stiffness of system (6.8) and in the EOS only the acoustic part is consistently large. Thus, only the longitudinal acoustic waves are consistently faster with respect to the deformation, whereas the shear elastic waves have a speed which is similar to the material deformation.

Physically, M ≤ M χ is always verified. Therefore, with the aim of recovering the correct numerical viscosity in both low Mach regimes, the proposed convex combination of upwind and centered fluxes (3.13) still holds for the all-speed scheme. For a detailed study of the low Mach regime in elastic solids and of the behaviour of the different waves, we refer the reader to [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF][START_REF] Abbate | Numerical methods for all-speed flows in fluid-dynamics and non-linear elasticity[END_REF]. In the following section we will solve two different Riemann problems representing the two limits inside two elastic solids. condition is the following

(p, ρ, u 1 , u 2 ) (x 1 , x 2 , t = 0) =          10 9
, 8900, 0, 0 , Ω 1 10 9 , 8900, 0, 0 , Ω 2 10 9 , 8900, 0, 0 , Ω 3 10 5 , 8900, 0, 0 , Ω 4 , where we have introduced the subdomains Ω It is evident that the all-speed spatial discretization accurately approximates the two material waves, which are instead diffused by an upwind-like discretization. Then, in Figs. 9(c)-9(d), we solve the same problem with the all-speed relaxation scheme by enforcing a material CFL constraint ν mat = 0.2, which gives ∆t = 7.8 • 10 -6 , and ν mat = 0.3, which gives ∆t = 1.17 • 10 -5 : the accuracy on the material waves is maintained also with large material time steps. Of course the fast waves are diffused, since their speed is too high to be followed by a material ∆t, as explained in Sec.

1 = [1, 2] × [1, 2], Ω 2 = [0, 1] × [1, 2], Ω 3 = [0, 1] × [0, 1] and Ω 4 = [1, 2] × [0,
3.2.1. In Fig. 10 we report the components σ 11 and σ 22 of the stress tensor, the pressure and the velocity field u, which is of course continuous on the two material waves. We report these profiles obtained with the acoustic CFL condition, in order to show all the propagating waves with a good resolution. We impose a small initial pressure discontinuity and also an initial velocity field in Ω 4 . The velocity field is imposed in order to analyze the propagation of slow This manuscript is for review purposes only. shear waves. In this framework, we get M 2.5 • 10 -3 and M χ 0.125 on the two material waves at the edges of Ω 4 . Thus, the longitudinal acoustic waves are consistently faster with respect to all other waves, including the material waves and the shear elastic waves. In fact, these latter waves propagate with a speed similar to the deformation velocity (this wave pattern is briefly presented in the previous section and analysed in detail in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF]). In Figs. 11(a)-11(b) we show the density contours at final time t = 3 • 10 -4 s: we compare the results obtained with the all-speed relaxation scheme and with the explicit-upwind relaxation scheme. For both schemes we enforce ν ac = 0.45, which gives ∆t = 5.6•10 -6 , in order to have a good resolution of all waves.

It is evident that the all-speed scheme is providing a more accurate approximation of the deformation waves also for this specific low Mach number limit. In Figs. 11(c)-11(d) the stress tensor component σ 21 is plotted and we can focus on the shear elastic waves approximation. At time t = 3 • 10 -4 s these slow waves are still very close to the material waves. We can observe that the sharpness and the correct magnitude of the shear waves are accurately approximated by the all-speed discretization, whereas the upwind scheme is once again diffusive. These results are consistent with the results of a similar 1D test proposed in [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF].

7. Conclusions and future developments. In this work we have proposed an asymptotic preserving scheme to solve flows of compressible materials at all speeds.

The scheme exploits the relaxation method proposed by Jin and Xin, in order to deal with an advective operator that is independent of the EOS. Thanks to this, we have been able to simulate fluid flows and elastic deformations without any modification

This manuscript is for review purposes only.

of the scheme structure.

A proof of the asymptotic preserving property has been proposed, showing that the incompressibility condition is respected by the scheme. The Gresho vortex tests have shown that, by setting an incompressible initial flow, the scheme is able to preserve the incompressible regime at all times, with pressure fluctuations of order M 2 .

The scheme has also been tested on the solution of Riemann problems, accurately approximating material waves propagating in both fluids and elastic materials. The material waves are kept sharp also with large material time steps. Moreover, an adaptive mesh refinement algorithm has been designed for the proposed scheme, providing a consistent reduction of the computational effort. The algorithm is based on the estimate of the scheme entropy production, which has also been used to study the accuracy of the proposed spatial discretization.

Different improvements will be proposed in the future, including the extension of the scheme to higher orders and the derivation of suitable preconditioners that can provide a more efficient solution of the linear system. The scheme will also be applied to the simulation of multi-material problems such as low speed-impacts. This can be done with the introduction of a level set function to track the physical interfaces and by extending the implicit multi-material scheme proposed in [START_REF] Abbate | An implicit scheme for moving walls and multi-material interfaces in weakly compressible materials[END_REF] to multi-dimensional problems. Moreover, the AMR technique will be extensively analyzed and improved.

here reported: u1

F (ψ) =                          0 1 0 0 -u 2 1 -σ 11
                   
and a similar expression is derived for G (ψ) in the same way.

  scheme (4.1)-(4.2)-(4.3)-(4.4) and we collect terms of equal power of M :

= 0 is respected. Relations ∇p n+1 0 = 0 and ∇p n+1 1 = 0 -

 010 Scheme (4.11)-(4.14)-(4.15) is clearly a consistent discretization of the Euler system in its incompressible limit, derived in (2.13)-(2.14)-(2.15). This means that the scheme is asymptotic preserving. Nevertheless, in what follows, we show that the incompressibility constraint ∇ • u n+1 0 = 0 is respected.

0= p 0 0 and the same for p n+1 1 .

 1 is independent of n too, i. e. p n+1 0 = p n 0 = ... = p 1 0 Of course, this means that also the energy is independent of n, namely (ρe) 17), one obtains directly the incompressibility constraint ∇ • u n+1 0 = 0. The proof of the AP property for the fully discrete scheme (3.16) easily follows by introducing discretization (3.13) into the spatial derivatives of (4.11)-(4.14)-(4.15), in

Fig. 1 . 1 :

 11 Fig. 1. Test 1: Gresho vortex with Mmax = 10 -1 (initial condition and results at time t = 1 of the relaxation all-speed scheme and with an explicit-upwind scheme).

(a) t = 2 (b) t = 3 Fig. 2 . 1 :

 2321 Fig. 2. Test 1: Gresho vortex with Mmax = 10 -1 , results at time t = 2 and t = 3 obtained with the relaxation all-speed scheme.

Fig. 3 .

 3 Fig. 3. Tests 2 and 3: Gresho vortex with Mmax = 10 -2 (first row) and with Mmax = 10 -3(second row). Initial condition and results at time t = 1 with the relaxation all-speed scheme and with an explicit-upwind scheme.

Fig. 4 .

 4 Fig. 4. Water Gresho vortex with Mmax = 10 -2 (initial condition and results at time t = 1 of the relaxation all-speed scheme and with the explicit-upwind relaxation scheme).

  it has been shown that, if the solution is locally smooth, S n+1 ij = O (h r ) with r equal to the order of the scheme. On the other hand, S n+1 ij = O (1) if there is a contact discontinuity and S n+1 ij = O (1/h) if there is a shock in the considered cell.

Fig. 5 .

 5 Fig. 5. 2D Riemann problem at time t = 0.3: comparison between the explicit-upwind relaxation scheme and the relaxation all-speed scheme. Density plots (30 contours: from 0.53 to 1.98) and numerical entropy production. In panel (a) R indicates a rarefaction, S a shock and C a contact wave.

  4, 1.0625, 0, 0.2145) , Ω 3 (0.4, 0.5197, 0, -1.1259) , Ω 4 , This manuscript is for review purposes only. where we have introduced the subdomains Ω 1 = [0.5, 1]×[0.5, 1], Ω 2 = [0, 0.5]×[0.5, 1], Ω 3 = [0, 0.5] × [0, 0.5] and Ω 4 = [0.5, 1] × [0, 0.5]. This test is characterized by the following wave pattern: two contact waves, one left moving shock and one right moving rarefaction: C 21 , ← -S 32 , C 34 , -→ R 41 .

= 9 Fig. 7 .

 97 Fig. 7. Grids for the 2D Riemann problem at time t = 0.3, obtained with a L min = 5 and with different Lmax (Lmax = 6: 3916 cells, Lmax = 7: 10030 cells, Lmax = 8: 22918 cells, Lmax = 9: 68251 cells).

  and if the level of refinement of cell C ij does not equal L max , then the cell is marked for refinement; this cell is thus split into four children and cell averages in the newly created cells are set by taking the cell average of the "ancestor"; 3. if |S n ij | < S coa for all four children and if the level of refinement of cell C ij does not equal L min then the cell is marked for coarsening; the four children are replaced by the ancestor cell and the cell average in the ancestor is set by taking the mean of the cell averages of the four children; 4. the time step ∆t is computed with the chosen CFL constraint using the smallest cell size of the grid.

Fig. 8 .

 8 Fig. 8. Density for the 2D Riemann problem at time t = 0.3. Comparison between uniform (number of cells: 1,048,576) and AMR grids (number of cells: 146,578).

  Forward and backward characteristics are related as follows: Y (X (ξ, t)) = ξ. By differentiating this latter relation, one gets(6.1)

Fig. 9 .

 9 Fig. 9. Copper Riemann problem: density profiles at time t = 10 -4 s (20 contours: from 8882 to 8935). Panel (a): explicit-upwind scheme with ∆t = 8.8 • 10 -7 , given by νac = 0.45. Panels (b), (c) and (d): all-speed relaxation scheme with ∆t = 8.8 • 10 -7 , ∆t = 7.8 • 10 -6 (given by νmat = 0.2) and ∆t = 1.17 • 10 -5 (given by νmat = 0.3) respectively.

5

 5 scheme, by enforcing the same acoustic CFL ν ac = 0.45, which gives ∆t = 8.8 • 10 -7 .

(a) σ 11 (b) σ 22

 1122

Fig. 10 . 8 ,

 108 Fig. 10. Copper Riemann problem: normal components of the stress tensor σ, pressure p and velocity field u at time t = 10 -4 s (νac = 0.45). The second test is instead representative of the acoustic low Mach regime. In order to simulate a rubber-like material, we adopt the EOS parameters γ = 4.4, χ = 8 • 10 5 Pa and p ∞ = 6.8 • 10 8 Pa and the following initial conditions:

Fig. 11 .

 11 Fig. 11. Rubber-like Riemann problem: density (30 contours: from 994 to 1005) and tangential stress σ 21 at time t = 3 • 10 -4 s (νac = 0.45).

,ψ 1 2u1 -σ 11 ,ψ 2 -σ 11 ,ψ 3 -σ 11 ,ψ 4 -u1u2 -σ 21 ,ψ 1 u2 -σ 21 ,ψ 2 u1 -σ 21 ,ψ 3 -σ 21 ,ψ 4 - 5 -u2σ 21 ,ψ 5 -u1σ 11 ,ψ 6 -u2σ 21 ,ψ 6 -u1σ 11 ,ψ 7 -u2σ 21 ,ψ 7 1 -σ 11 ,ψ 8

 11121131142112122132145215116216117217118

Table 1

 1 Total kinetic energy E kin,tot (t = 1) after one full rotation of the Gresho vortex relative to its initial value E kin,tot (t = 0) for different maximum Mach numbers.

Table 2

 2 Global pressure fluctuations p f l after one full rotation of the Gresho vortex for different maximum Mach numbers.

  includes different physical behaviors. The first part is exactly the internal energy (2.3) for gases and liquids. The last term describes the variation of energy in a neohookean elastic solid due to elastic deformations (χ is the shear elastic modulus). As shown in Table3, classical models are obtained by specific choices of the coefficients.

	Material	γ	p ∞	χ
			[Pa]	[Pa]
	Perfect biatomic gas	1.4	0	0
	Stiffened gas (water) 4.4	6.8 • 10 8	0
	Elastic solid (copper) 4.22 3.42 • 10 10 5 • 10 10

Table 3

 3 Typical parameters for different materials.

  4) inside the conservation equations, in order to separate terms scaling with M and with M χ . The scaled

	state law reads as follows
	(6.9)	ρe M 2 =	1 2	ρ|u| 2 +	p + γp ∞ M 2 (γ -1)	+	χ trB -2 2M 2

χ , where again both the acoustic and the elastic Mach numbers are present. Two different low Mach number regimes can be observed: 1. acoustic and elastic low Mach regime, namely M 1 and M χ 1. This case can be verified only if the parameters of the considered material are such that O(p ∞ ) O(χ). We thus get that the two Mach numbers are of the same order, i. e. O(M ) O(M χ ). Copper can be representative of this limit, since p ∞ = 3.42 • 10 10 Pa and χ = 5 • 10 10 Pa. In this regime, the gradient of stress tensor σ fully scales as O(M 2 ), i. e. the entire stress tensor gradient introduces stiffness in system (6.8). All propagating waves are consistently faster with respect to the material waves, since both the acoustic and elastic terms of the EOS are large and of the same order of magnitude. 2. acoustic low Mach regime, namely M 1 and M M χ . This is verified if the material parameters are such that p ∞ χ, giving an acoustic Mach

This manuscript is for review purposes only.

Acknowledgments. E. Abbate acknowledges the support of the Università Italo-Francese/Université Franco-Italienne ("Bando Vinci 2016"). G. Puppo and E. Abbate were partially supported by the "National Group for Scientific Computation (GNCS-INDAM), 2018". The numerical experiments were carried out using the PLAFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from Bordeaux INP, LABRI and IMB and other entities: Conseil Régional d'Aquitaine, Université de Bordeaux and CNRS (and ANR in accordance to the Programme d'Investissements d'Avenir http://www.plafrim.fr/).

Appendix A. Jacobian of the fluxes in elasticity. The jacobians of the fluxes F (ψ) and G (ψ) inside the Eulerian model (6.2) can be computed analytically.

We introduce the notation σ jk ,ψi , which stands for the derivative of the jk , j, k = 1, 2 component of the tensor σ with respect to the conservative variable ψ i , i = 1, ..8.

These derivatives are computed using definition (6.4) of σ, where it can be useful to rewrite the pressure as function of the energy as follows

We remind that tensor σ is symmetric, i.e. σ 12 = σ 21 . The expression for F (ψ) is

This manuscript is for review purposes only.