IMPLICIT DISCRETIZATION OF LAGRANGIAN GAS DYNAMICS

. We construct an original framework based on convex analysis to prove the existence and uniqueness of a solution to a class of implicit numerical schemes. We propose an application of this general framework in the case of a new non linear implicit scheme for the 1D Lagrangian gas dynamics equations. We provide numerical illustrations that corroborate our proof of unconditional stability for this non linear implicit scheme.


Introduction
To approach equations traducing the movements of compressible fluids, explicit schemes are traditionally used, see [17,30]. Explicit schemes need to satisfy a stability CFL condition such as ∆ ≤ ∆ , where is the speed of sound, ∆ is the time step, and ∆ is the size of the discretization in space of the mesh. In some cases, this CFL constraint contributes to have such small time steps that it becomes unfavourable to use explicit methods.
An alternative is to use implicit in time schemes which have always aroused interest in the literature. In particular, they are much less sensitive to the CFL number: for example a finite difference algorithm is proposed in [2], implicit and semi-implicit schemes are investigated in [37] and an experimentation on implicit upwind methods for Euler equations is done in [26]. Some other implicit algorithms are explained in the following articles, see for example [10,25,29]. A large part of them use the method of predictor-corrector scheme like in [21,27,[39][40][41]. More recent works can be found in [7,28]. A major reference in the context of our work is [15] where an implicit Lagrangian scheme for non viscous compressible gas dynamics is studied for astrophysical purposes, but only by means of numerical experiments and without further theoretical foundation.
Nonetheless, major technical difficulties appear for the numerical resolution of fully implicit non linear schemes. At a theoretical level, it is difficult to prove the existence and uniqueness of a solution to implicit schemes. Currently, a powerful theory is the one developed in [16], for Navier-Stokes equations, using the topological degree. The existence of a solution is proved in details, but the uniqueness requires restrictive hypothesis.
Another strategy is explained in [4] for piecewise linear functions in the case of a symmetrical structure of the linear part of the system. The existence of a solution is proved, and the uniqueness is also studied in the case of special hypothesis. The non-linear implicit-explicit strategy in [15] is second order in both space and time, but there are no proof of existence or uniqueness even if the numerical results indicate good robustness. In some sense, our work answers positively to the theoretical issue raised in [15] about the construction of fully justified implicit solvers.
Our original contributions to this field in this work are, firstly the elaboration of a general framework which allows to prove existence and uniqueness of implicit solution of some numerical schemes, and secondly the application of the general framework in the case of a non linear implicit scheme for the 1D model problem (1) written in semi-Lagrangian coordinates (semi-Lagrangian coordinates means Lagrange+update) (1) One has = 1 > 0 the mass density, is the pressure, is the velocity and is the total energy density. The variables and are taken positive to be physically admissible, see [34] or [13] for more details. The material derivative used in (1) and afterwards is = + . The following set of equations is the isentropic Euler equations that is approximated by the prediction step of our implicit scheme (2) The first two equations are identical to (1), only the last one is different, with denoting the physical entropy.
To simplify, these equations are equipped with a perfect gas equation of state where is the thermal capacity at constant volume, > 1 is the adiabatic index, = − 1 2 2 is the internal energy density, is the temperature and is the speed of sound given by 2 = .
The justification of using a prediction step based on the discretization of (2) comes from ideas in the work of Chalons, Coquel and Marmignon [5]. It will be explained in details in this article.
Consider a mesh ℳ composed of cells noted ∈ {1, . . . , }. The time is discretized with a time step ∆ that corresponds to one iteration. The mass of the cell is . The boundary conditions are supposed periodic, and the fluxes are defined with the acoustic impedance + 1 2 > 0. The scheme has a predictor-corrector structure. The prediction step is written as The correction step is given by where only the total energy is modified. The correction step is explicit so the main difficulty is in the prediction step.
In (4) and (5), the fluxes are defined by where the coefficient + 1 2 > 0 is for simplicity equal to a mean value of the acoustic impedance: In Lagrangian formalism, the mesh moves according to the velocity of the fluid: +1 . The predictor-corrector scheme (4-5) is naturally conservative since it is expressed in terms of fluxes. Such a scheme can be proved to be weakly consistent as in Després [11,12].
The predictor-corrector scheme is an adaptation of ideas from the article [5] which is dedicated to solve the Euler equation (6) in Eulerian coordinates The authors explain that the difficulties of solving this system come from the flux terms in the second and third equations. Indeed, there is a strong non linearity due in particular to the pressure. To overcome this complexity, the authors propose a predictor-corrector strategy that we use also in this work. Firstly is to solve the isentropic Euler equation (7) during the prediction step Secondly, the classical Euler equation (6) are solved in order to restore the conservation of the total energy. At a discrete level, the fluxes are expressed thanks to an isentropic scheme, and then inserted in the scheme associated to (6). For the prediction step, the authors use a relaxation scheme on the pressure. They prove the existence of a solution to the relaxation implicit scheme. Nonetheless, the robustness of the scheme depends on an extra equation as mentionned in a report, see [33].
Let us now describe our results. For physically admissible data, the system (4) will be said unconditionally stable if there exists a unique solution to the implicit non-linear scheme. Our proof of the unconditional stability of (4) comes from a rewriting of the prediction step (4) under the form where is a vector of real unknowns, is a functional defined on a domain and is a matrix of real coefficients. In our case, is strictly convex over its domain and is skew-symmetric, so the transformation ↦ → ∇ ( ) − is a monotone operator, see Brézis [3]. The proof that (8) has a unique solution relies on Theorem 1.4 that seems to be new considering the classical literature of convex analysis, see [1,18,20]. The ingredients to establish Theorem 1.4 are the following.
We made a slight abuse of notations by using the same letter in the Hypothesis 1.1 as the iteration index in the scheme (4)(5). We believe this does not interfere with the readability.
Hypothesis 1.2. The function : ∈ → ( ) ∈ R is 3 , strictly convex and coercive in the sense that Moreover for all ∈ there exists a unit direction ∈ R + which is outward from such that Also for all ∈ , one has ||∇ ( )|| The verification of (9), (10) and (11) will be obtained directly from the perfect gas law equations (3). For an isentropic gas, it can be simplified. Applying this Theorem, we show that (4) is well defined for all ∆ > 0.
Corollary 1.5. Considering physically admissible data ( > 0 for all ), the prediction scheme (4) can be written under the form (8). Therefore, it is unconditionally stable.
Moreover, the predictor-corrector scheme (4-5) satisfies two entropy inequalities. Theorem 1.6. For all ∆ > 0, the solution of the prediction step (4) satisfies The solution of the correction step (5) verifies the entropy inequality 5 The case = 0 corresponds to one unknown systems. For example the Traffic flow equations where the only unknown is the density. Otherwise, > 0 and > 0.
The organization of this article is as follows. In Section 2 we write the scheme (4) under the form (8). In Section 3, we prove Theorem 1.4. The proof is split in several parts. On the one side we rapidly prove the uniqueness of a solution, and on the other side we decompose the proof of the existence in different steps. In Section 4, we apply Theorem 1.4 for the isentropic Euler equations, and prove Corollary 1.5. In Section 5, the correction step is introduced and the complete scheme is fully justified. The proof of Theorem 1.6 concerning entropy inequalities for both steps is detailed. In Section 6, we provide a few numerical illustrations. The Appendix contains a brief description of the modification to treat an isothermal equation of state.

Formulation under the form (8)
The objective of this Section is to provide the details of the transformation from the implicit scheme (4) to the form (8). The verification of Hypothesis 1.1, 1.2 and 1.3 will be performed in Section 4.
We consider Euler isentropic equations in one dimension (2) for compressible perfect gas with periodic boundary conditions. As the physical entropy is constant during this prediction step, its equation is not necessary for the proof. Replacing the fluxes and rearranging the terms in (4), one obtains Let us define a vector of unknowns where the domain is defined as One then has = = using notations of Hypothesis 1.1. The matrix A is defined by The functional : → R is defined as a sum of elementary functions where the elementary functions are Replacing the equation of state (3) by another one leads to a new definition of the functions 1 and 2 . The functions 1 and 2 depend only on the scheme and remain the same. Proof. For a perfect gas law, the correspondence between , and can be written as = ( −1) exp ( ) Therefore the equivalence between a solution of (14) and a solution (15) of ∇ ( ) = is explicited by To finish the proof it is sufficient to calculate explicitly ∇ ( ). The derivatives of 1 , 2 , 1 and 2 are By (18), one calculates all the components of the vector ∇ ( ) ∈ R 2 . With the definition (17) of the matrix , one obtains immediately that the first equations in the vectorial identity ∇ ( ) = are while the last equations are With the correspondence (19), one finds that (20 -21) is equal to (14).

Proof of Theorem 1.4
In this Section, we prove Theorem 1.4 stated in the introduction under the Hypothesis 1.
By construction, is lower semi-continuous because is continuous over . For a function which is coercive on its domain like in (9), the closure is also coercive in the sense of Hirriart-Urruty [18], [Chapter 2, p 41]

Uniqueness
It relies on elementary considerations which are classical for monotone operators [3].
Lemma 3.1. Assuming that the problem (8) admits a solution in , then it is unique.
Proof. Let 1 ∈ and 2 ∈ be two solutions of the problem (8) One has Since is a skew-symmetric matrix therefore ( Since is strictly convex, this is only satisfied if 1 = 2 .

Existence
The existence of a solution relies on a few intermediate results which are convenient for our model problem.

Existence of a minimum for
The first result to prove is the existence of a minimum point for the function , using a classical result from convex analysis. Proof. We use [Definition 3.2.6, p. 180] of Hirriart-Urruty and Lemarechal [19] to the function = . So there exists ∈ R + such that ( ) ≤ ( ) for all ∈ R + . Necessarily ( ) < ∞ is finite so ∈ . It remains to show that / ∈ . Let us assume on the contrary that ∈ . Thanks to the convexity of and inequality (10) in Hypothesis 1.2, one can write It is a contradiction. Therefore ∈ is a minimum and is unique thanks to the strict convexity of on .

A continuation method
We prove in this Section that the problem (23) admits a solution in the domain for all 0 ≤ ≤ 1 For = 0, the problem is treated in Section 3.2.1. This allows to write (23) with a continuation method under the form of an initial value problem (24) A rearrangement yields (∇ 2 ( ) − ) = . For ∈ , the matrix ∇ 2 ( ) − is invertible thanks to the following result.

Lemma 3.3. Let and
be two matrices of ℳ , (R), > 0 ∈ N, such that is skew symmetric and is positive definite. Then the matrix = + ∈ ℳ , (R) is invertible.
The Lemma 3.3 is applied with − ∈ ℳ + , + (R) as the skew symmetric matrix, and ∇ 2 ( ) ∈ ℳ + , + (R) as the positive definite matrix. Thus (∇ 2 ( ) − ) −1 exists. The initial value problem can be rewritten as Let us define ℐ = [0, +∞[ and the function : To obtain the existence of a maximal solution to (25), one can apply standard results from the theory of ODEs that we recall now.
Definition 3.4 (see Coddington and Levinson [6], Chap. 1). Let ∈ N and : ℐ × R → R , let 0 ∈ ℐ, and ∈ R where ℐ is a non empty interval of R. A solution of the differential equation is given by a non empty interval ⊂ ℐ and a differentiable function : → R satisfying (26) for all ∈ . A solution of the initial value problem (or Cauchy problem) associated to (26) is a solution of (26) such that 0 ∈ and ( 0 ) = 0 .
Theorem 3.5 (Cauchy Lipschitz Theorem for locally Lipschitz functions [6], Th. 3.1, p 12). Let the function be a 1 function, then, for all initial data ( 0 , 0 ) ∈ ℐ × R , there exists an interval ∈ ℐ containing 0 such that there exists in a unique solution to the associated initial value problem.
In particular for all such data, there exists a unique maximal solution and all other solution verifying the condition of Cauchy is a restriction of the maximal solution. Proof. One applies the Cauchy Lipschitz Theorem. The function is well defined, and differentiable in terms of . The derivative is continuous because is a matrix of scalar coefficients, and is a function of class 2 thanks to Hypothesis 1.2. In terms of the second variable, as ∇ 2 is locally Lipschitz, and the other terms are locally bounded, is then of class 1 . Thanks to Theorem 3.5, the problem (25) admits a unique maximal solution. To prove the last part of the Lemma, one notes that Since ∇ ( 0 ) = 0 it shows that ∇ ( ) − = 0 on the maximal interval, (23) is satisfied.
In the rest of this Section, we prove that max > 1.

Upper bound of ( )
In this Section, the objective is to prove that ( ) is bounded, which is necessary to conclude that the solution of the problem (23) stays in the domain .
Lemma 3.7. There exists ∈ such that the following inequality is satisfied on the maximal interval Proof. Let us take ∈ ker( ) ∩ that is non empty by Hypothesis 1.3. The convexity of implies that The matrix is skew symmetric, hence ( , ) = 0. So Using again the property of skew symmetry, one has ( In addition, since is a coercive function by Hypothesis 1.2, there exists < +∞ such that for all in the maximal interval.

End of the proof of Theorem 1.4
The end of the proof of Theorem 1.4 is based on the following standard result. For our problem, one has the Property.
Proof. Thanks to (27), It remains to prove that is a compact of . Let us take a sequence ∈ for ∈ N. Since ( ) is bounded, there exists ∈ (0, ) and a subsequence still denoted such that → . Necessarily ∈ , so either ∈ or ∈ . Let us assume that ∈ . Thanks to Hypothesis 1.2, inequality (11), one has ||∇ ( )|| → +∞. It is a contradiction with the definition of . Therefore ∈ . Since is 2 , ∇ is a continuous function and ||∇ ( )|| ≤ 2|| || . So ∈ , which shows that is a compact of .
Proof of Theorem 1.4. Thanks to Proposition 3.9 and Theorem 3.8, one has that max > 1. Therefore, one takes = 1 which concludes the proof of existence of the solution of (8).

Proof of the unconditional stability of Corollary 1.5
We prove the unconditional stability for the prediction step corresponding to the implicit discretization of the isentropic Euler equations. In this purpose we apply Theorem 1.4 after a precise check of the hypothesis. The definitions of , and are given in Section 2.

Properties of J
We verify the properties of Hypothesis 1.2, that is the strict convexity of , its coercivity and its limits at the boundary of the domain .   (18) is strictly convex on .
The proof is easily verified, but we detail the calculations.
Proof. The second derivatives, for all and ∈ {1, . . . , } are Hence, for all ∈ , and for all ∈ R 2 , one has For ̸ = 0, one has (∇ 2 ( ) , ) > 0. Therefore the function is strictly convex on . Proof. The function (18) is the sum of elementary functions 1 , 2 , 1 and 2 . The quadratic functions 1 and 2 are clearly bounded from below, as well as 2 . The function 1 is also bounded from below because > 0, > 0, > 1 and 1− 1 is dominated by for → +∞. So 1 is coercive. It is evident that 2 is coercive. Since is defined by a sum over all indices 1 ≤ ≤ , then is coercive.  (10) and (11).
Proof. The first derivative of with respect to − is Let ∈ . It means that there exists a non empty subset ⊂ {1, . . . , } such that for all ∈ , = 0. One takes as the unit outward direction ∈ R + , such that ∀ ∈ , > 0 and for all ̸ ∈ , = 0. The limit of the first derivative of when → 0 + is ∀ ∈ Indeed, lim

Properties of matrix A
Let us prove Hypothesis 1.3 holds.
Property 4.6. The kernel of and the domain intersect.
Proof. We first evaluate the kernel of the matrix . Let ∈ R 2 satisfying = 0, so One obtains  (17) and (18). All the hypothesis of Theorem 1.4 are satisfied. The existence and uniqueness of a solution to the implicit isentropic Euler scheme for all time step is proved.

Proof of the entropy inequalities (Thm. 1.6)
In this Section, we study the two steps of the predictor-corrector scheme, and prove the stability inequality in each step.
We recall the complete scheme Each inequality is studied separately.

Proof of stability inequality (12)
Actually, inequality (12) is a mathematical entropy inequality. Usually, entropy stability is defined for continuous in time problems, or for explicit in time schemes. A continuous definition is found in [30], [Thm. 3.3, p 27], and a discrete version of entropy stability is written in [12], [Prop. 2.25, p. 66].
Definition 5.1. The implicit conservative scheme (4) is said entropic if there exists a mathematical entropy pair ( , ) and a numerical entropy flux Φ( , ) (such that Φ( , ) = ( )) such that the following inequality holds for all We use this definition onto the scheme of the prediction step (4) where the vector of unknowns is = ( , , ) .
The entropy corresponding to the isentropic system is ( ) = . As is convex, one gets ( . In Lagrangian formalism, the calculus can be performed easily. Thanks to Gibbs formula, see [30], [Chapter 4, p. 318], one has = − + .
For the flux, the one state solver is End of the proof of inequality (12). First is the evaluation of ∇ ( ) · + 1 2 , and second ∇ ( ) · − 1 2 . Then the difference between the two expressions is performed. In the rest of the calculations, the time dependence will be omitted to simplify the writing. One has and with the same method of calculus The difference is One has the equalities 1 2 ( + ) These results are injected in the previous inequality One then denotes also written differently as This corresponds exactly to the entropy inequality because Φ is the entropic flux. As a matter of fact, one remarks Therefore as = and Φ = , one rewrites (29) as

Proof of inequality (13). Let us denote =
During the correction step, the discretization of the total energy is ) .
The variation of total energy is evaluated between the correction and the prediction step as To conclude, it is important to have in mind the Gibbs formula = + , where the variable is fixed because +1 = = 1 . One has During the prediction step = , so Thus, thanks to (30), ( , ) is a growing function of for fixed. One concludes that +1 − Δ ≥ 0.

Numerical illustrations
In this Section, we provide numerical illustrations which show that the theoretical properties of the numerical methods are transferred to real calculations. The implicit scheme is solved using a Newton algorithm and the final update of the solution is performed in a conservative way, so the scheme is implemented in a perfectly conservative fashion, as for finite volume methods [14,23,36]. For all our test problems, we have observed that the Newton algorithm converges without any difficulties in only few iterations (approximately 5) and we do not comment this issue further.
The numerical illustrations can be divided between those related to robustness issues and those related to accuracy issues. Robustness issues are illustrated in all numerical simulations, from reasonable CFL numbers (CFL = 0.4) to huge ones (CFL = 537). Accuracy issues are illustrated in Section 6.1.1 (gas dynamics with CFL = 0.4 to CFL = 537). The position of the contact discontinuity is discussed in Section 6.1.2. The authors also provide a simulations performed on non uniform meshes, with stiffened gas in Section 6.2.1, and a perturbed Sod shock tube in Section 6.2.2.

Fully implicit treatment
For each example, the implicit scheme is compared to the explicit acoustic scheme (31) which provides a reference solution Explicit scheme where the fluxes are defined by ).

Sod shock tube
For this problem, the initial conditions are The boundary conditions are = ℎ = 0. The adiabatic index is = 1.4. The equation of state for the gas is given by (3). The final time of the simulation is = 0.2. Several CFL (0.4, 40, 80, 537) are taken to evaluate the robustness of the scheme. The number of cells is 100. To give an idea of the computational time, the explicit solver is six times faster than the implicit solver for this test case.
For a CFL equal to 0.4 the curves are quasi identical in Figure 1. There are oscillations on the density and the entropy curves near contact discontinuities. This is a classical phenomenon for explicit Lagrangian schemes named wall heating. The implicit scheme does not correct these oscillations. Numerous papers have been dedicated to the correction of these oscillations, the interested reader can see [31] or [38]. This phenomenon is not discussed further in the rest of this paper.
For an implicit CFL 100 times larger, the numerical smearing is visible in Figure 2 for the rarefaction waves as well as the shock. On the contrary, the contact discontinuity is still at the correct location.
As it can be seen in Figure 3, one observes that when the CFL tends to be very large, there is more numerical dissipation on shocks and rarefaction waves but the contact discontinuity still seems to be at the right position.
The solution of Figure 4 shows the unconditional stability of the implicit scheme.  A remark can be done on the position of the contact discontinuity that is slightly shifted. To explain the origin of this misplacement, we can evoke the fact that the interaction with the boundaries of the domain is very important. To validate this hypothesis, we performed another calculation (same initial conditions and final time) on a domain 9 times larger. The results are visible in Figure 5, and the contact discontinuity is once again well positioned.

Position of the contact discontinuity
We develop hereafter a possible explanation for the precision of the position for the contact discontinuity in Figure 5. It deals with the integration of the Riemann problem.
Indeed, consider the following initial conditions The equations (32) are discretized in time but the space part is left continuous, ∆ > 0, ∈ R. This mimics the implicit scheme, indeed ∆ can be taken extremely big, but the space step ∆ is very small. This corresponds of having a discrete ∆ compared to a small continuous ∆ . One obtains Using the boundary conditions (33), the integration constant is = .
One finally obtains For > 0, the solution verifies the equations We apply the same method than in the case < 0, and check that the solution is of the same kind. One finally finds an expression for ′ ′ ( ) = ± √︀ −2 ( , ) + + .
At the interface, when = 0, the continuity conditions are with ⋆ ∈ R. One gets Using (35) and (36), one finds the scalar equation where ⋆ is the unknown.
In the numerical examples, we took initial conditions of a Sod shock tube, that are recalled hereafter.

As
Using the perfect gas law, one can rewrite the equation in terms of and .
One obtains  Numerically, we calculated with a Newton method that the solution to (37) is approximately equal to ⋆ = 0.2559. It corresponds to a velocity of ⋆ = 0.8789. The exact value of the velocity for the Riemann problem at the contact discontinuity is = 0.9275. This value is found in Toro [36], [ Table 4.3, p 131]. The difference between ⋆ and is equal to 5.2%, which is a satisfying accuracy considering that the implicit simulation performs with only one time step. This small relative error of 5.2% is the reason why the contact discontinuity of the implicit solver is approximately superimposed with the reference one in Figure 5.We observed a similar behavior for all the other test problems and we believe it is a strong asset of this family of implicit Lagrangian schemes.

Blast wave problem
To give the reader an idea of how the implicit scheme behaves in other configurations, we present here a more severe test case called the blast wave problem. It is also known as the Woodward and Colella problem, see [42]. This test case involves multiple interactions. Indeed, two strong blast waves develop right after the beginning of the simulation, then collide at = 0.028 s, inducing a decrease of the time step for explicit schemes, and at the final time = 0.038 s, a new contact discontinuity can be observed.
For this example, the initial conditions consist in three constant states of perfect gas at very different pressures. The simulation is performed on a mesh of 400 cells between [0, 1], following the article [43], and there are reflexive boundary conditions on both sides of the domain. As the shocks are very strong, it would not be relevant to take large time steps. In Figure 6, the precision of the implicit scheme is visible for a same CFL than the explicit scheme. In Figure 7, we increase the CFL of the implicit scheme and we observe some oscillations in the density curve for = 5. The wall heating is visible in the density curve for both the explicit and the implicit scheme at the location of the contact discontinuity. For strong shocks, it seems that the numerical dispersion is sensible to the CFL number. We do not comment this result further as we do not know if it comes from a numerical artifact.

Implicit-Explicit coupling
When more than one fluid is represented, it is interesting to have a cell size appropriated to each fluid. It can lead to great disparities in the dimension of the cells and hence an implicit-explicit coupling is a solution. The mesh is separated into subdomains. Each one is treated using either the acoustic explicit solver, or the implicit prediction-correction scheme developed in this work. The mesh is considered as described in Figure 8, namely the implicit part of the mesh contains the smaller cells.  At each interface between an explicit and an implicit cell, the values of the fluxes must be the same. Consider there is an interface between cell and cell + 1 of a mesh ℳ, the problem at the interface is the following where and are respectively the values of the pressure and the velocity in cell at time . The values obtained by the prediction step of the implicit scheme in cell + 1 are +1 for the pressure +1 for the velocity. The fluxes at the interface are * For the computation, the implicit part of the mesh ℳ imp is treated first using a Newton algorithm to obtain ( ) ∈ℳimp and ( ) ∈ℳimp . At the end of this step, the values of the fluxes are evaluated on the implicit and explicit part of the mesh then used to update , and at time +1 in each cell of ℳ.
Let us emphasize that the implicit-explicit coupling described in this Section is not an IMEX method. Indeed, the IMEX strategy consists in a numerical methods which contain an implicit part, but as far as we know, the non linear global part is always treated in an explicit manner as in [35] and references therein. In our case, the strategy differs in the sense that the scheme used is totally implicit or totally explicit depending on the subdomains treated.

Water-gas simulation
To validate this coupling, we present a water-gas simulation and compare the results to the solution obtained with a total explicit solving. For this example, one considers the case of a stiffened gas provided with the The coefficient describes the attractive effects that lead to a cohesion in the matter, it is also called the reference pressure and must satisfies > 0. A modified expression of is thus evaluated and implemented = (︂ ( − 1) exp( ) Theorem 1.4 still applies on the variable = (−( + ) ∈ℳimp , ( ) ∈ℳimp ). The two-phase shock test case presented originates from [32]. It considers having two materials with all the variables strongly discontinuous.
On the left part of the tube there is water (high pressure) and on the right part air (low pressure). The initial conditions are The variable is set to 0 ( ) = 6 × 10 8 for < 0.7. The simulation is performed on a mesh of 1000 cells. There are 950 cells between [0,0.7] and 50 cells between [0.7, 1]. The smaller cells lie in the left region that is thus solved using the implicit scheme, and the right part is solved with the explicit acoustic solver. The explicit solver reaches the final time = 240 × 10 −6 in 2160 iterations and a time step of = 1.11 × 10 −7 s. The implicit-explicit solver runs during 585 iterations and a time step of = 3.9 × 10 −7 s. The computational time is approximately the same. In Figure 9 one notices that the rarefaction wave is more dissipated with the implicit-explicit treatment. The contact discontinuity and the shock are well placed. This validates the implicit-explicit coupling.

Sod shock tube perturbed with a water drop
We present here an original test case where the implicit-explicit coupling algorithm is more efficient than the explicit acoustic solver of reference. It consists of a Sod shock tube perturbed by a water drop. In Figure 10, the explicit curve and the implicit-explicit one are similar in shape with the reference solution. The gas hits the water drop from the left side, creating an important reflexive pressure wave as can be seen in Figure 10. It is well modeled by both of the methods. The explicit solution is evaluated in 65161 iterations in time, corresponding to a = 2.45 × 10 −9 . The time of computation is around 76.5 s. The implicit-explicit solution is obtained in 158 iterations in time, with a time step dictated by the size of the bigger explicit cells, corresponding to = 2.58 × 10 −7 , and a computational time of about 4.5 s. For this type of test case, the implicit-explicit coupling performs well.

Conclusions
We have used a strategy of predictor-corrector scheme, based on the previous work [5] in Eulerian coordinates, to solve numerically the Euler equations. We have defined an abstract frame in order to analyze a family of implicit schemes written under the peculiar form (8). We have proved the existence and uniqueness of a solution to the prediction step of our implicit scheme. We provided examples using this result and led numerical tests that have indeed corroborated the theoretical statements of stability. The numerical illustrations compare the implicit scheme to an explicit scheme of reference, and show the precision of this new algorithm in these cases. It also provides examples in the case of an implicit-explicit coupling for stiffened gas.
In a future work, it would be interesting to generalize this method to the case of thin elasto-plastic structures using Kluth and Després [22] or Maire et al. [24]. We could also try to improve this work by using a more elaborate flux or increase the scheme order at the order 2. It will probably ameliorate the precision, but it stays to evaluate the cost of simulation that it would generate. The multi-dimensional version would need a more advanced management for the displacement of the mesh, but the principal ingredients of Theorem 1.4 remain similar.
Other numerical examples that are more realistic have to be performed to evaluate the pertinence of this algorithm. Theoretically as well, it would be great to have an explanation on the rapid convergence of the Newton algorithm for the prediction step, and to have a more elegant proof of the entropy inequalities using the frame (8).
Finally our approach can be the basis of a fully implicit Lagrange+remap strategy for the development of implicit solvers for the non viscous Euler system, where it is sufficient to treat the linear remap stage in an implicit fashion to obtain a fully implicit Eulerian numerical scheme. The evaluation of such approaches in particular in the context of low-Mach flows will be the topic of further examination.