Second‑Order Invariant Domain Preserving ALE Approximation of Euler Equations

An invariant domain preserving arbitrary Lagrangian-Eulerian method for solving non-linear hyperbolic systems is developed. The numerical scheme is explicit in time and the approximation in space is done with continuous finite elements. The method is made invariant domain preserving for the Euler equations using convex limiting and is tested on various benchmarks


Introduction
An arbitrary Lagrangian-Eulerian (ALE) technique is developed in this paper in an attempt to combine the advantages of classical kinematic descriptions, and to minimize their drawbacks as much as possible. The computational mesh given by an ALE formulation, allows us to handle greater distortions of the fluid with more precision than that obtained by a purely Eulerian formulation. In our previous work, Guermond et al. [7,9], we investigated the artificial viscosity problem with hyperbolic systems using an ALE formulation, continuous finite elements and explicit time stepping. In this work, we propose a first-order artificial viscosity that does not depend on any ad hoc parameters and that leads to precise 1 3 invariant domain properties and entropy inequalities. Combining a low-order method with an entropy consistent high-order method via a convex limiting process, we describe a highorder method that preserves the invariant domains.
We consider the following hyperbolic system in conservative form: where the dependent variable u is ℝ m -valued, the flux is f ∈ C 1 (ℝ m ;ℝ m×d ) , and u 0 is any admissible initial data. For simplicity in the presentation of the theoretical results, we consider periodic boundary conditions or assume that u 0 is constant outside a compact set. Some examples of well-known equations that can be written in this general form are the Euler equations, shallow water equations, and Burgers' equation.
The main objective of this paper is to present a simplified version of our high-order ALE schemes for Euler equations, mainly focusing on the computational aspects and key features of the ALE motion description and omitting the details of the precise mathematical formulation necessary to rigorously demonstrate the properties of the methods. For this reason, this paper describes in details the algorithm and the implementation of our high-order invariant domain preserving method in a finite element code. One of the key issues in the definition of the ALE motion is the mesh optimization process. We solve it in two stages: (i) high-order reconstruction of the Lagrangian velocity; (ii) smoothing of the mesh. The smoothing stage is driven by information provided by a simulation. We define a blending parameter to smooth the mesh using an adaptation of the method proposed by Loubère et al. [11] which has also been used by other research groups (see, e.g., Boscheri and Balsara [1], Boscheri et al. [2]). There are other ways to approximate the ALE movement, where the ALE velocity is defined not only using the velocity of the fluid, but also more information from the simulation, such as the gradients of the variables or shock positions, see e.g., Dobrev et al. [4], Zhao et al. [15], but we do not explore this direction in this paper. One novelty here, with respect to Guermond et al. [7,9], is that we present new numerical tests in which we compare results obtained using ALE and Eulerian formulations of our high-order scheme.
The remainder of this paper is organized as follows. The ALE motion and its approximation are described in details in Sect. 2. Specifically, we present the main difference between purely Lagrangian and ALE descriptions and describe the finite element spaces used to approximate the ALE motion as well as our method of computing the ALE velocity. To compute the ALE velocity, we reconstruct the approximation of the Lagrangian velocity and then apply a smoothing procedure. The general ALE scheme is introduced in Sect. 3, where the low-order invariant domain preserving method and the high-order (possibly invariant domain violating) method are briefly described. The convex limiting technique, whose purpose is to make the high-order method invariant domain preserving, is described in Sect. 4. The new ALE method is illustrated numerically in Sect. 5 on various benchmark problems and compared with the same method in Eulerian coordinates.

ALE Motion and Its Approximation
We consider an ALE motion X A ∶ ℝ d × ℝ + → ℝ d as a uniformly Lipschitz mapping, and define the spatial description of the associated velocity as In the notation used herein, the subindex A stands for arbitrary. (Fig. 1). To illustrate the differences between ALE and Lagrangian descriptions, we consider the transport equation

Remark 1 Arbitrary Lagrangian Eulerian vs. Lagrangian motion
where is the velocity of the fluid.
By computing the time derivative of u along the characteristic curves of the Lagrangian motion, we get If we discretize (4) using the explicit Euler method, we obtain where x i = X( , t i ), i = n, n + 1 ; that means that the solution is translated along the characteristic curves.
By writing the material derivative along the characteristic curves of the ALE motion, we get Therefore, using the explicit Euler method, we obtain Fig. 1 Description of an ALE motion and its relation with the Lagrangian motion of a body where X n+1 A ( ) = x n+1 . The new solution is the solution at the previous time step moved along the characteristics plus a correction term due to the difference between the ALE velocity and the fluid velocity.

Geometric Finite Elements and Mesh
We now present the discrete setting of the ALE motion. Let T n h , h > 0 , n = 0, 1, ⋯ , N be a sequence of shape-regular matching meshes, where T 0 h denotes the initial mesh and T n h is the deformed mesh by means of the user-defined ALE motion at time t n . Given mesh T n h , we denote by D n the computational domain generated by T n h . To simplify the presentation, we assume that our meshes are formed by triangles in the case of two space dimensions and tetrahedra in the case of three space dimensions. In Guermond et al. [7], a more general setting for the configuration of the meshes that can be used was described.
We introduce the reference Lagrangian finite element (K,P geo ,Σ geo ) , where K is the reference element, P geo is the reference polynomial space, and Σ geo is the set of linear forms that define the degrees of freedom. We denote by {â i } i∈{1∶N geo } and {̂ geo i } i∈{1∶N geo } the Lagrange nodes of K and the associated Lagrange shape functions, respectively, where N geo ∶= dimP geo . Let {a i } i∈{1∶V geo } be the collection of all Lagrange nodes of mesh T n h and let geo ∶ T n h ×{1∶N geo } ⟶ {1∶V geo } be the geometric connectivity array (assumed to be independent of the time). Then the geometric transformation T n K ∶K → K ∈ T n h is defined by We now define the finite-dimensional spaces based on the geometric Lagrange finite elements and the following vector-valued spaces which we use to represent the ALE motion: By considering the approximate ALE motion mapping X n,n+1 That is, the approximate ALE motion is fully described by the position of geometric Lagrange nodes at time t n+1 : a n+1 i ∶= X n,n+1 A (a n i ). In the numerical simulations presented at the end of this paper, we use triangles and P geo (T n h ) is one of the following polynomial spaces: ℙ i , i = 1, 2, 4.

Approximation of the ALE Motion
Although the ALE motion of the mesh, or equivalently the ALE velocity field, is userdefined, we give in this section some details on how we compute the ALE velocity which (10) X n,n+1

3
may be useful to the reader. All the numerical results presented in this paper have been obtained using the techniques described below. We compute the new position of the geometric Lagrange nodes of the mesh using an explicit Euler step, that is where v n A is the approximate ALE velocity at time t n and Δt ∶= t n+1 − t n is the time step. We implement the calculation of the ALE velocity v n A in two stages: reconstruction and smoothing of the ALE velocity. A brief description of the two stages is as follows.

Reconstruction of the ALE Velocity
One of the motivations for using ALE methods in compressible hydrodynamics is that by making the mesh motion as close as possible to the actual fluid motion, one can significantly reduce the effects of the artificial viscosity. Hence, at each time t n , the approximate solution of the hyperbolic system should be used to reconstruct the fluid velocity, the ALE velocity should be defined from that. Let v n be the approximated velocity of the fluid obtained by solving our hyperbolic system. The easiest way to define the ALE velocity is v n A = v n which make the scheme Lagrangian. However, if we use low-order finite elements to solve our system, v n A , and therefore X n,n+1 A , would be also in a low-order space and would not represent well smooth mesh motions with vorticity. It has been reported in the literature and is also our experience that using higher-order polynomials to represent the mesh motion limits the risks of mesh entangling and, in the case of vortical motions, postpones the time when the mesh eventually entangles. To illustrate this observation, we solve an isentropic vortex problem 1 where the analytical velocity is known, and we move the mesh using two different approximations of this velocity: piecewise linear and piecewise quadratic. In Fig. 2, we plot the mesh at the instant (11) X n,n+1 (a n i ) = a n+1 i = a n i + Δtv n A (a n i ), before collapsing, which occurs at t = 1.81 when using linear polynomial approximation of the velocity and at t = 3.48 when using high-order (quadratic) polynomial approximation of the velocity. Reconstructing a high-order field from a low-order representation is a problem that has been thoroughly investigated in the computer graphics literature. We propose to reconstruct the ALE velocity from v n using an algorithm called the butterfly subdivision algorithm, initially described in Dyn et al. [5] for surface reconstruction. We refer the reader to our previous work [9] for more details on how this algorithm is used to reconstruct the ALE velocity.
Once the velocity is reconstructed, we compute the new position of the mesh nodes using where v n Lag ∈ P geo (T n h ) is the reconstructed velocity. The notation Lag is used to indicate that nodes are moved following the fluid motion as no mesh smoothing is applied yet.

Smoothing of the ALE Velocity (Maintaining the Connectivity)
When the fluid undergoes large deformations, for instance, when a vortex or a shock wave is formed, a purely Lagrangian algorithm can experience a loss of accuracy or the mesh will eventually tangle. Therefore, it is not possible to describe the computational domain motion without doing remeshing operations. One way to avoid this breakdown is to use in the algorithm an ALE velocity which is a smoothed version of the Lagrangian velocity.
We propose a method that involves blending the Lagrangian velocity using an averaging technique. The blending is done through a parameter that controls the local deformation of the cells; the definition of this parameter is inspired by a technique proposed by Loubére et al. [11]. From the Lagrangian motion given by (12), we compute an averaged version a n+1 i,Sm as follows: where i ∈ V geo , and I(i) denotes the collection of indices of the neighboring nodes of i. This smoothing technique is not sufficient to avoid mesh tangling when the fluid motion is complex; a more sophisticated curvilinear mesh smoothing is a topic of our ongoing research. Finally, the actual position of the geometric nodes is defined by To define the blending parameter, we first construct first the Jacobian matrix of the map- . Then, we compute the right Cauchy-Green strain tensor | K |K and the two eigenvalues 1,K (a n i ) , 2,K (a n i ) at all the geometric Lagrange nodes a n i ∈ K , with the convention 1,K (a n i ) ⩽ 2,K (a n i ) . We define T i ∶= {K ∈ T n h | a n i ∈ K} for all i ∈ V geo , and compute the blending parameter as follows: i,Lag = a n i + Δtv n Lag (a n i ), In each application, we will specify this natural number, which can be larger if we want our ALE motion to be closer to the Lagrangian motion. Finally, after the two stages, the ALE motion is given by and the ALE velocity of the geometric nodes is given by

ALE Scheme
We describe in this section two ALE approximations of (1). One is invariant domain preserving and entropy satisfying, while the other is at least second-order accurate in space but may violate the invariant domain property. We use continuous finite elements and explicit time stepping.
If is a weak solution of the hyperbolic system (1), then the following result that we demonstrated in Guermond et al. [7] is the main motivation for the ALE formulation.

Lemma 1 The following identity holds in the distribution sense
where (x, t) ∶= (X −1 A,t (x)).

Approximating Finite Elements
The construction of the approximate solution of (1) is based on a reference finite element {(K,P,Σ)} . It is important to note that P geo and P are different objects. (K,P geo ,Σ geo ) is a Lagrange element but (K,P,Σ) may not be; it may for instance be a Bernstein-Bezier finite element, see for example Lai and Schumaker [10, Chap. 2] or another modal finite element.
In our code, we use Lagrange finite elements; specifically we always take P = ℙ 1 but we use The shape functions of the reference element are denoted by

and introduce the vector-valued spaces
The global shape functions in P(T n h ) are denoted by { n i } i∈V . Recall that these functions form a basis of P(T n h ) . We denote by ∶ N×T n h ⟶ V the connectivity array associated with the global shape functions { n i } i∈V . This array, which we assume to be independent of n, is defined such that n (i,K) (x) =̂ i ((T n K ) −1 (x)) , for all i ∈ N and for all K ∈ T n h . For any i ∈ V , I(i) is the collection of indices of the shape functions whose support has a nontrivial intersection with the shape function i ; that is, we set where for any measurable set E ⊂ D n , |E| denotes the measure of E. Henceforth we call the connectivity graph of P(T n h ) , the graph (V, E) where the vertices are all members of V , and the edges are pairs (i, j) in V 2 such that (i, j) is in E iff j ∈ I(i) and i ∈ I(j) . Notice that, actually, j ∈ I(i) iff i ∈ I(j) . The connectivity graph does not depend on n, since we assumed that the connectivity array ∶ N×T n h ⟶ V does not depend on n. The solution of (1) is approximated in P m (T n h ) , so we can write

Generic Algorithm
We introduce in this section a generic algorithm to obtain u n+1 n , the approximate solution of (1) at time t n+1 . This algorithm was introduced in Guermond et al. [7] and was sightly modified in Guermond et al. [9] to allow the ALE velocity to be represented with a polynomial degree much higher than that used to approximate the solution of (1). A detailed description of this method and its conservation properties can be found in those references; here we describe only the main steps of the method.
Let ( 0 i ) i∈V be the approximation of the mass of the shape functions at time t 0 defined by 0 be a reasonable approximation of the initial data u 0 , and let ( n i ) i∈V be the approximations of the mass of the shape functions at time t n . If N is the number of time steps, for each n ∈ {0, 1, ⋯ , N} , we do the following steps.
(i) Compute the ALE velocity v n A . This can be done following the stages described in Sect. 2.2.
(ii) Move the mesh: a n+1 i = a n i + Δtv n A (a n i ).
and d n ij is an artificial viscosity for the pair of degrees of freedom (i, j), that is clearly defined in Sects. 3.3 and 3.4. We call d n ij the graph viscosity since this coefficient only involves the connectivity graph of P(T n h ) . We henceforth assume that d n ij = 0 if j ∉ I(i) and Notice that d n ii does not really need to be defined for (24) to make sense; this quantity is nevertheless introduced to shorten the definition of the CFL number.

Remark 2 Recall that the initialization is done by setting
At this point we assume again that the user-defined ALE velocity field w n is reasonable or that Δt is small enough that n+1 i is positive.
Δt is the forward Euler approximation of the left-hand side in (18). Notice that we have replaced the consistent mass matrix by the approximate lumped mass matrix to approximate the time derivative. The expression ∑ j∈I(i) ( n j ⊗ n j − ( n j )) ⋅ c n ij is just the Galerkin approximation of the right-hand side of (18).

First-Order Viscosity: GMS-GV Method
We define in this section the artificial graph viscosity that makes the method described in the previous section to be invariant domain preserving and entropy satisfying. The method is called GMS-GV because it is based on the guaranteed maximum speed (GMS) and firstorder graph viscosity (GV); see Guermond and Popov [6], Guermond et al. [7,8].
We consider the 1D Riemann problem: (26) d L,n ij = max( max (g n j , n n ij , n i , n j )‖c n ij ‖ 2 , max (g n i , n n ji , n j , n i )‖c n ji ‖ 2 ).

Second-Order Viscosity: GMS-EV Method
We now describe a technique that is formally high-order accurate in space but may be invariant domain violating. Let ( ( ), ( )) be an entropy pair for the hyperbolic system.

Convex Limiting Technique: a Second-Order Invariant Domain Preserving Scheme for Euler Equations
Some limitations must be applied to the high-order update, H,n+1 i , given by (21) and (22) with the viscosity given by (30) to guarantee the preservation of the physical bounds and to be invariant domain preserving.
Given a convex invariant set B ⊂ A , and assuming that n i ∈ B, ∀i ∈ V and some n ⩾ 0 , we explain in this section how to push back H,n+1 i in the invariant domain using a convex limiting technique. For simplicity, we will give a brief overview of the limiting strategy for Euler equations; the precise description of the method for a general hyperbolic system and the rigorous mathematical framework can be seen in Guermond et al. [9]. We note that when the equation of state (EOS) implies finite compressibility of the fluid (like with the co-volume EOS), we can automatically enforce that with our method and this seems to be difficult with other formulations. We start by estimating the difference H,n i − L,n i . This is done by subtracting (22), written with the high-order viscosity d H,n ij , from (22), written with the low-order viscosity d L,n ij . We obtain Now we apply this result to limit the high-order scheme to solve Euler equations with an ideal gas EOS given by p = ( − 1)(E − 1 2 ‖u‖ 2 2 ) where > 1 . Euler equations can be written as with = ( , , E) T ∈ ℝ d+2 and ( ) ∶= ( , ⊗ + p , (E + p)) T . The specific energy is That means that, for any smooth ALE velocity, the scheme is not only invariant domain preserving ( n+1 i ∈ A ), but also preserves the local domain ( n+1 i ∈ B i ).

Remark 4
In the numerical tests presented in Sect. 5, the quantity s n,min i is relaxed as explained in Guermond et al. [8,Sect. 4.7.2], to maintain the second-order accuracy of the method.

Remark 5
Imposing local lower and upper bounds on the density, with the bounds coming from the bar sates, automatically enforces positivity and in case of the co-volume EOS, it also enforces the maximum compressibility of the method.

Numerical Results
In this section, we consider a set of test problems for the Euler equations, which are designed to verify the theoretical conservation and invariant domain properties of our highorder ALE finite element scheme. We first solve problems with known solutions to compare the errors obtained with our low-and high-order methods (with and without limiting). Then, we present a more challenging set of problems, where the mesh motion is more complex and we focus the presentation on the robustness of our ALE method. All simulations use RK3 strong stability preserving time integration derived from the explicit Euler scheme given in (22). The EOS for all test cases is the ideal gas law, p = ( − 1) e , where the value of is given for each test case. When the exact solution of a specific test case is known, we compute the following relative error indicator: The above norm is estimated using local Gaussian quadrature rules of order 8 with 16 points.

Isentropic Vortex
The first test case we consider is the so-called isentropic vortex problem. The exact solution is isentropic and is given by The free-steam conditions we use are ∞ = p ∞ = T ∞ = 1 and u ∞ = (2, 0) T . The perturbations are where r = ‖ − c (t)‖ is the Euclidean distance from the vortex center c (t) ∶= (x 0 1 + 2t, x 0 2 ) T , = 5 is a constant defining the vortex strength, and = 7 5 . We set the initial computational domain to be D 0 = (−5, 5) 2 . The first mesh consists of 20×20 squares divided into two triangles, then the mesh is refined uniformly five times. We set the density, momentum and total energy to the free-stream values on the boundary of the computational domain at all times. We use the GMS-EV method with the entropy commutator computed with the generalized entropy (u) = p 1 to set d H,n ij , see Sect. 3.4. Table 1 displays the errors obtained with our ALE schemes at t = 2 , where the CFL number is 0.25. We observe the first-order accuracy with the GMS-GV method and second-order accuracy with the GMS-EV method with and without limiting. The solution of the isentropic vortex is smooth so the limiting technique does not offer an additional advantage over the initial GMS-EV method; however, we have verified that it does not significantly increase the computational cost or increase the error of the method.
The CPU time used for the three methods is provided in Table 2. It can be observed that the limiting technique increases the CPU time by about 3% when the GMS-EV is used. (46) Notice that the time step changes in each simulation, because it depends on the viscosity and the fluid velocity. Therefore, a better indicator of the computational cost of each method would be the ratio CPU time∕(N#dof ) where N is the number of time steps done.
From the values of this parameter displayed in Table 2, we can deduce that the computational cost of the low-and high-order methods is similar. Finally, to illustrate how the accuracy of the methods increases when the ALE description is used, and Table 3 displayed the error indicator obtained with the GMS-EV scheme with limiting using the Eulerian (setting v A = 0 ) and ALE approaches. When the Eulerian approach is used, as the mesh is fixed and the vortex is moving to x > 0 direction, we have to increase the computational domain to ensure that the vortex remains in the domain. We set D 0 = (−5, 10) × (−5, 5) . As in the previous simulations, we consider meshes formed by squares divided into two triangles, and we use the same meshes for both techniques. It can be observed that the errors obtained with the ALE approach are smaller than those obtained with the Eulerian approach, whereas the rate of convergence is the same.

Noh Problem
Noh problem is a test problem with a known exact solution; a shock wave propagating radially outwards at a constant speed given by  The initial computational domain is D 0 = (−1, 1) × (−1, 1) and we do the computations up to t f = 0.6 , with CFL = 0.4 and = 5 3 . The calculations were performed in a sequence of uniform meshes generated from N×N, N ∈ {30, 60, 120, 240} squares divided into two triangles. A summary of the data used in our simulations is provided in Table 4.
Convergence rates and L 1 relative errors are presented in Table 5. With both schemes, GMS-GV and GMS-EV limited, the convergence rates are close to first order which is consistent with the literature on this problem. As expected, the errors obtained using the highorder scheme are smaller. However, when using the GMS-EV high-order scheme without applying the limiting technique, the convergence rates are worse. This is consistent with the results shown in Fig. 3. On the left panel of Fig. 3, a scatter plot of the density field is plotted for three ALE methods: GMS-GV, GMS-EV and GMS-EV with convex limiting. To obtain these plots, we compute the distance to the center of each mesh vertex and plot it versus the value of the density at this point. We compare our results with the exact solution. If we look at the results obtained with the GMS-EV with and without limiting, plotted with green and red solid lines, respectively, we observe that the non-limited solution is more precise, but it exhibits undershoots and overshoots, as expected.
To examine the behavior of our ALE technique, Fig. 4 compares the solution obtained with our best scheme (GMS-EV + limiting) with the same method using an Eulerian description of the motion. It can be concluded that the solution with the Eulerian description is much more diffusive.

3
To illustrate the robustness of our ALE method with respect to the mesh regularity, we solve the problem on a non-uniform mesh. The initial mesh is divided into four quadrants: the bottom left quadrant is composed of 32×32 squares; the top left is composed of 32×64 squares; the top right is composed of 64×64 ; and the bottom right is composed of 64×32 squares. Each square cell is divided into two triangles. Figure 5 indicates that our ALE high-order method preserves the radial symmetry of the solution and does not develop any hour-glass-like instability with a non-uniform mesh.

Radial Sod Problem
We run the well-known Sod shock tube problem using a radial configuration, so the initial condition is ij . This high-order viscosity, together with the limiting technique, allows us to use a purely Lagrangian motion for the mesh without smoothing the velocity. Thus, for this test problem, in the definition of the ALE velocity, we only do the reconstruction stage using ℙ 2 finite elements.
The results for the ALE (Lagrangian in this case) and Eulerian GMS-EV schemes are presented in Fig. 6. The density fields for both methods are plotted together with the final mesh obtained by the ALE method. The Eulerian versions of the scheme suffer from a strong numerical diffusion, whereas with the ALE method, the mesh moves to better capture the shock wave and the contact discontinuity. The same results can be observed from the surface plots of the density in Fig. 7.

Double Mach Reflection Problem
The double Mach reflection problem involves a strong shock (Mach 10) moving in air ( = 1.4 ) that impinges a wall with 60 degree angle. When the shock runs up the wall, a self-similar structure (with two triple points) emerges. This is a standard benchmark test popularized by Woodward and Colella [14] and then by many others for evaluating the resolution of Euler codes. Different setups are proposed in the literature for this problem to prevent the formation of undesirable numerical artifacts (see Vevek et al. [13]). The setup in the case of an ALE code is even more problematic, we use a setting similar to that of Zhao et al. [15] but allow the non-wall boundaries to move with the flow.
We solve this problem in the initial computational domain D 0 = (−1.5, 4) × (0, 1) , and with the shock initially positioned at x 1 = 1∕6 where the horizontal wall begins. A free-slip boundary condition is enforced on x 2 = 0, x 1 ⩾ 1∕6, at the top boundary, the flow values are set to describe the exact motion of the Mach 10 shock, and inflow and outflow boundary conditions are enforced on the left and right boundaries, respectively. The initial condition is 3(x 1 − 1 6 ) − 20t is the position of the shock at instant t. The flow is computed at time t = 0.2 with CFL = 0.5 . The initial mesh is a uniform mesh formed by 65 142 triangles and 33 072 ℙ 1 nodes. Regarding the mesh motion, we allow the inflow and the top boundary nodes to move right with the flow. Because part of the bottom boundary is a wall, the second component of the velocity in the wall nodes is 0. Therefore, in order to maintain the computational domain as a square, we set on all boundaries ( A ) x 2 = 0. Figure 8 presents the density field obtained with the GMS-EV scheme using ALE and Eulerian approaches. In the works of Zhao et al. [15] and Boscheri and Balsara [1], the double Mach reflection problem is also solved using ALE techniques. The mesh regularization stage in our ALE movement is much simpler than that used in these two references where a much more sophisticated mesh motion is defined. However, we believe that our results are in good agreement with the results of Boscheri and Balsara [15], Zhao et al. [1] as the final results for the density are similar. The result from the GMS-EV scheme using ALE mesh motion appears to be superior to the Eulerian one. The triple points are resolved better and the jet is well defined (Fig. 9).

Rayleigh-Taylor Instability
This instability occurs in the interface between two fluids of different densities when the more dense fluid is on top. Under the force of gravity, the initial equilibrium is unstable to any perturbation of the interface and the heavier fluid goes into the lighter one developing the so-called Rayleigh-Taylor instability. The computational domain is a fixed rectangular box D 0 = (0, d∕2) × (0, 3d) . In the simulations, we set d = 1∕3 . The gravitational force is = (0, −0.1x 2 ) T . The  heavy fluid has density max = 2 and the light fluid has density min = 1 . The interface between the two fluids at t = 0 is x 2 = (x) = 0.5 − 0.1d cos(2πx 1 ∕d) . The growth of the perturbation at the interface is exponential and the rate is exp( t) where ∝ A , and A = ( max − min )∕( max + min ) is the Atwood number; here, we take A = 1∕3 . The initial density field is slightly regularized by setting The initial pressure is hydrostatic. The slip boundary condition is enforced on the four walls of D 0 . In Fig. 10, we can see the computational domain and the initial conditions, as well as the details of the initial mesh around the interface. The mesh is composed of 50×200 squares divided into two triangles. The ALE velocity is computed with P geo = ℙ 4,2 Lagrange finite elements. The simulations are run until t = 7.5 with CFL = 0.6 . A summary of the data used in the simulations is given in Table 6.
The mesh undergoes very large deformations in the time interval [0, t] and is smoothed using s = 3 to define the blending parameter given by (15). Figure 11 displays the density field obtained with the GMS-EV scheme with limiting using ALE and Eulerian approaches. Finally, to illustrate that the method is robust when using fourth-order polynomials to describe the mesh motion, that is P geo = ℙ 4,2 , Fig. 12 displays the meshes and the density field at six time: t = 2, 3, 4, 5, 6, 7.5 . These pictures show that the mesh undergoes very large deformations in this test.   We do not compare here the solutions obtained by the different schemes proposed, as it was done in the previous sections. This is because when we try to solve this problem with the low-order method, with the mesh described above, the interface between the two fluids is completely diffused. On the other hand, when we use the GMS-EV scheme (not limited), the solution is not between the physical bounds. As a consequence the mesh collapses at a certain time (before the final time of t = 7.5 s ), even when doing the smoothing, because the ALE velocity is not smooth enough to keep the mesh untangled.

Conclusion
We have extended in this paper our work from Guermond et al. [7,9]. The main novelty is that we give complete details on the implementation of the method for the Euler equations and present new numerical tests. In particular, we compare the performance of the ALE and Eulerian formulations of our high-order scheme for the radial Sod and the double Mach reflection problems. The results clearly show the advantage of using the ALE method and confirm the robustness of our methodology. More sophisticated mesh motions and handling of multiple materials will be addressed in future works.