Thomas O Gallou

Quentin M Érigot
email: quentin.merigot@math.u-psud.fr

A LAGRANGIAN SCHEME À LA BRENIER FOR THE INCOMPRESSIBLE EULER EQUATIONS

Keywords: 1991 Mathematics Subject Classification. 35Q31, 65M12, 65M50, 65Z05 Incompressible Euler equations, Optimal Transport, Lagrangian numerical scheme, Hamiltonian

We approximate the regular solutions of the incompressible Euler equations by the solution of ODEs on finite-dimensional spaces. Our approach combines Arnold's interpretation of the solution of the Euler equations for incompressible and inviscid fluids as geodesics in the space of measure-preserving diffeomorphisms, and an extrinsic approximation of the equations of geodesics due to Brenier. Using recently developed semi-discrete optimal transport solvers, this approach yields a numerical scheme which is able to handle problems of realistic size in 2D. Our purpose in this article is to establish the convergence of this scheme towards regular solutions of the incompressible Euler equations, and to provide numerical experiments on a few simple test cases in 2D.

Introduction

The purpose of this article is to investigate a discretization of Euler's equation for incompressible and inviscid fluids in a domain Ω ⊆ R d with Neumann boundary conditions:

             ∂ t v(t, x) + (v(t, x) • ∇) v(t, x) = -∇p(t, x), for t ∈ [0, T ], x ∈ Ω , div (v(t, x)) = 0 for t ∈ [0, T ], x ∈ Ω , v(t, x) • n = 0 for t ∈ [0, T ], x ∈ ∂Ω , v(0, x) = v 0 .
(1.1)

As noticed by Arnold [START_REF] Arnold | Sur la géométrie différentielle des groupes de lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits[END_REF], when expressed in Lagrangian coordinates, Euler's equations can be interpreted as the equation of geodesics in the infinite-dimensional group of measurepreserving diffeomorphisms of Ω. To see this, consider the flow map φ : [0, T ] × Ω → Ω induced by the vector field v, that is:

       d dt φ(t, x) = v (t, φ(t, x)) for t ∈ [0, T ], x ∈ Ω , φ(0, •) = id, ∂ t φ(0, •) = v 0 .
(1.2)

Using the formula d dt det Dφ(t, x) = div (v(t, x)) det Dφ(t, x), the incompressibility constraint div (v(t, x)) = 0 and the initial condition φ(0) = id, one can check that φ(t, •) belongs to the set of volume preserving maps S, defined by

S = s ∈ L 2 (Ω, R d ) | s # Leb = Leb ,
where Leb is the restriction of the Lebesgue measure to the domain Ω and where the pushforward measure s # Leb is defined by the formula s # Leb(A) = Leb(s -1 (A)) for every measurable subset A of Ω. Euler's equations (1.1) can therefore be reformulated as

             d 2 dt 2 φ(t) = -∇p(t, φ(t, x)) for t ∈ [0, T ], x ∈ Ω , φ(t, •) ∈ S for t ∈ [0, T ], φ(0, •) = id, ∂ t φ(0, •) = v 0 . (1.3)
To obtain (1.3) one simply needs to derive (1.2). This equation can be formally interpreted as the equation of geodesics in S. In particular the pressure term in the evolution equation in (1.3) expresses that the acceleration of φ should be orthogonal to the tangent plane to S at φ. Indeed, note that the condition φ(t, •) ∈ S in (1.1) encodes the infinitesimal conditions div v(t, •) = 0 and v(t, x) • n(x) = 0 in (1.3). This suggests that the tangent plane to S at a point φ ∈ S should be the set {v • φ | v ∈ H div (Ω)}, where H div (Ω) denotes the set of divergence-free vector fields

H div (Ω) = v ∈ L 2 (Ω, R d ) | Ω v • ∇ϕ = 0, ∀ϕ ∈ C ∞ c (Ω) .
In addition, by the Helmoltz-Hodge decomposition, the orthogonal subspace to H div (Ω) in L 2 (Ω, R d ) is the space of gradients of functions. Therefore the evolution equation in (1.3) expresses that d 2 dt 2 φ(t) ⊥ T φ(t) S , in other words that t → φ(t, •) is a geodesic of S. Note however that a solution to (1.3) does not need to be a minimizing geodesic between φ(0, •) and φ(T, •). The problem of finding a minimizing geodesic on S between two measure preserving maps amounts to solving equations (1.3), where the initial condition ∂ t φ(0, •) = v 0 is replaced by a prescribed coupling between the position of particles at initial and final times. It leads to generalized and non-deterministic solutions introduced by Brenier [START_REF] Brenier | The least action principle and the related concept of generalized flows for incompressible perfect fluids[END_REF], where particles are allowed to split and cross. Shnirelman showed that this phenomenon can happen even when the measure-preserving maps φ(0, •) and φ(T, •) are diffeomorphisms of Ω [START_REF] Shnirelman | Generalized fluid flows, their approximation and applications[END_REF].

Previous work: discretization of geodesics in S. The first numerical experiments to recover generalized minimizing geodesics have been performed by Brenier in 1D [START_REF] Brenier | Generalized solutions and hydrostatic approximation of the Euler equations[END_REF]. He also proposed a scheme to compute the solutions of the Cauchy problem (1.3) in [START_REF] Brenier | A combinatorial algorithm for the Euler equations of incompressible flows[END_REF]. In Brenier's discretization, the measure-preserving maps are approximated by permutations of a decomposition of the domain into cubes. The numerical implementation of this idea relies on the resolution of a linear assignment problem at every timestep, whose cost is unfortunately prohibitive for domains in dimension higher than one.

The discretization we consider in this article is a variant of this approach which is more tractable computationally and leads to slightly better convergence estimates. As in [START_REF] Brenier | Derivation of the Euler equations from a caricature of Coulomb interaction[END_REF], the measure-preserving property (or incompressibility) is enforced through a penalization term involving the squared distance to the set of measure-preserving maps S. This squared distance can be computed efficiently thanks to recently developed numerical solvers for optimal transport problems between probability densities and finitely-supported probability measures [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF][START_REF] Mérigot | A multiscale approach to optimal transport[END_REF][START_REF] De Goes | Blue noise through optimal transport[END_REF][START_REF] Lévy | A numerical algorithm for L 2 semi-discrete optimal transport in 3d[END_REF]. This alternative discretization has already been used successfully to compute minimizing geodesics between measure-preserving maps in [START_REF] Mérigot | Minimal geodesics along volume preserving maps, through semidiscrete optimal transport[END_REF], allowing the recovery of non-deterministic solutions to Euler's equations predicted by Shnirelman and Brenier in dimension two. The object of this article is to study whether this strategy can be used to construct Lagrangian schemes for the more classical Cauchy problem for the Euler's equations (1.1), able to cope with problems of realistic size in dimension two.

Discretization in space: approximate geodesics. The construction of approximate geodesics presented here is strongly inspired by a particle scheme introduced by Brenier [START_REF] Brenier | Derivation of the Euler equations from a caricature of Coulomb interaction[END_REF]. We first approximate the Hilbert space M = L 2 (Ω, R d ) by finite dimensional subspaces. Let N be an integer and let P N be a tessellation of Ω into N subsets (ω i ) 1≤i≤N satisfying

     ∀i ∈ {1, . . . , N }, Leb(ω i ) = 1 N Leb(Ω) h N := max 1≤i≤N diam(ω i ) ≤ C N 1/d
where C > 0 is independent of N . We consider M N ⊆ M the space of functions from Ω to R d which are constant on each of the subdomains (ω i ). To construct our approximate geodesics, we consider the squared distance to the set S ⊆ M of measure-preserving maps:

d 2 S : m ∈ M n → min s∈S m -s 2 M .
The approximate geodesic model is described by the differential equation in the finitedimensional space

M N :    m(t) + ∇d 2 S (m(t)) 2 2 = 0, for t ∈ [0, T ] , (m(0), ṁ(0)) ∈ M 2 N (1.4)
Note that the squared distance d 2 S is semi-concave, so that its restriction to the finitedimensional space M N is differentiable at almost every point. This differential system is induced by the Hamiltonian H :

M N × M N → R H(m, ṁ) = 1 2 ṁ 2 M + d 2 S (m) 2 2 .
(1.5)

We now rewrite the differential system (1.4) in terms of projection on the sets S and M N . Since the space of measure-preserving maps S is closed but not convex, any point in M admits a projection on S, but this projection is usually not uniquely defined. To simplify the exposition we will nonetheless associate to any point m ∈ M one of its projection P S (m), i.e. any point in S such that P S (m) -m M = d S (m). We also denote P M N : M → M N the orthogonal projection on the linear subspace M N ⊆ M, which is a linear map. We can rewrite Eq. (1.4) in terms of these two projection operators:

   m(t) + m(t)-P M N •P S (m(t)) 2 = 0, for t > 0 , (m(0), ṁ(0)) ∈ M 2 N (1.6)
From Proposition 5.2, the double projection P M N • P S (m) is uniquely defined for almost every m ∈ M N .

Remark 1.1. Equation (1.6) can be rewritten as a system of N particles in interaction, whose positions are denoted M 1 (t), . . . , M N (t) ∈ R d . Denoting 1 ω i the indicator function of the set ω i ⊆ Ω, we introduce

W : (M 1 , . . . , M N ) ∈ (R d ) N → d 2 S ( i M i 1 ω i ),
and we denote

B i (M 1 , . . . , M N ) = ∇ M i W (M 1 , . . . , M N ).
As explained in Proposition 5.2, the points (B i (M 1 , . . . , M N )) i are barycenters of a decomposition of Ω into N cells which depend on the solution to the optimal transport problem between Leb and the empirical measure

1 N 1≤i≤N δ M i . With these notations, Equation (1.6) is then equivalent to Mi (t) + 1 2 (M i (t) -B i (M 1 (t), . . . , M N (t))) = 0, for t > 0 and i ∈ {1, . . . , N }, (M (0), Ṁ (0)) ∈ (R d ) N × (R d ) N (1.7)
Loosely speaking, equations (1.4)-(1.6) describe a physical system where each particle M i (t) is subject to the force of a spring with stiffness 1 attached to the point B i (M 1 (t), . . . , M N (t)) which varies in time and depends on the position of all the particles. Equation (1.7) is also the Hamiltonian system associated to H :

(R d ) N × (R d ) N → R H(M, Ṁ ) = 1 2 N i=1 | Ṁi | 2 + W(M ), (1.8)
In the case of an non-homogeneous fluid with varying volume masse, such as a mixture of oil and water, an analogue discretization would involve a system of particles with different masses ρ i . This corresponds to replacing the Hamiltonian by

H(M, Ṁ ) = 1 2 N i=1 ρ i | Ṁi | 2 + W(M ).
(1.9)

In this last formulation, it is also possible to add potential terms, such as gravitation. This will be the case for the simulation of the Rayleigh-Taylor instability in subsection 5.4.

We first prove that the system of equations (1.4) can be used to approximate regular solutions to Euler's equations (1.1). Our proof of convergence uses a modulated energy technique which is similar to that used in [START_REF] Brenier | Derivation of the Euler equations from a caricature of Coulomb interaction[END_REF] and requires C 1,1 regularity assumptions on the solution to Euler's equations. See also [START_REF] Brenier | A geometric approximation to the euler equations: the vlasov-mongeampere system[END_REF][START_REF] Cullen | The semigeostrophic equations discretized in reference and dual variables[END_REF] for related works.

m(0) = P M N (id), ṁ(0) = P M N (v(0, •)),
which is twice differentiable and satisfies the second-order equation (1.4) for all times in [0, T ], possibly up to a countable number of exceptions. Then,

max t∈[0,T ] ṁ -v(t, φ(t, •)) 2 M ≤ C 1 h 2 N ε 2 + C 2 ε 2 + C 3 h N (1.10)
where the constants C 1 , C 2 and C 3 only depend on Ω, on the L ∞ norm (in space) of the velocity v(t, •) and on the Lipschitz norms (in space) of the velocity and its first derivatives ∇v(t, •), ∂ t v(t, •) and of the pressure and its derivatives p(t, •), ∇p(t, •), ∂ t p(t, •).

The values of C 1 , C 2 and C 3 are given more precisely at the end of Section 3. Note that the hypothesis on the solution m to the differential equation (1.4) is introduced here mainly for technical reasons. Removing it is not of our main concern in this paper since we also give a proof of convergence of the fully discrete numerical scheme regardless of this assumption. It is likely that solutions to (1.4) satisfying this hypothesis can be constructed through di Perna-Lions or Bouchut-Ambrosio theory [START_REF] Ambrosio | Transport equation and cauchy problem for BV vector fields[END_REF][START_REF] Bouchut | Renormalized solutions to the vlasov equation with coefficients of bounded variation[END_REF][START_REF] Lions | Sur les équations différentielles ordinaires et les équations de transport[END_REF], see also [START_REF] Brenier | A geometric approximation to the euler equations: the vlasov-mongeampere system[END_REF]Appendix].

Remark 1.3. Remark that (1.10) implies the convergence of the associated flows. In particular integrating inequality (1.10) one can show that max

t∈[0,T ] m(t) -φ(t) 2 M ≤ 2h 2 N + 2T C 1 h 2 N ε 2 + C 2 ε 2 + C 3 h N .
Discretization in space and time. To obtain a numerical scheme we also need to discretize in time the Hamiltonian system (1.6). For simplicity of the analysis, we consider a simple first-order scheme called symplectic Euler scheme with timestep τ > 0. The double projection P M N • P S (m) is defined as above. The discrete solution consists of two sequences M n , V n in the finite-dimensional space M N , given by:

     (M 0 , V 0 ) ∈ M N V n+1 = V n -τ 2 (M n -P M N • P S (M n )) M n+1 = M n + τ V n+1 (1.11)
Note that numerically, the piecewise-constant map M n : Ω → R d (resp. the piecewiseconstant vector field V n : Ω → R d ) is simply encoded by an ordered list of N points (resp. N vectors), so that this scheme can be considered as describing a dynamical system involving N particles. We have the following theorem, where we denote t n = nτ . Theorem 1.4. Let Ω be a bounded domain of R d with Lipschitz boundary, let and τ be positive numbers and let N ∈ N. Let v, p be a strong solution of (1.1), let φ be the flow map induced by v (see (1.2)) and assume that v, p, ∂ t v, ∂ t p, ∇v and ∇p are Lipschitz on Ω, uniformly on [0, T ]. Let (M n , V n ) n≥0 be a sequence generated by (1.11) from

M 0 = P M N (id), V 0 = P M N (v(0, •)). Assuming τ ≤ ε and h N ≤ ε, we have max n∈N∩[0,T /τ ] V n -v(t n , φ(t n , •)) 2 M ≤ C 2 + h N + h 2 N 2 + τ 2 ,
where the constant C only depends on Ω, on the L ∞ norm (in space) of the velocity v(t, •) and on the Lipschitz norms (in space) of the velocity and its first derivatives ∇v(t, •), ∂ t v(t, •) and of the pressure and its derivatives p(t, •), ∇p(t, •), ∂ t p(t, •).

In order to use the numerical scheme (1.11), one needs to be able to compute the double projection operator P M N • P S or equivalently the gradient of the squared distance d 2 S for (almost every) m in M N . Brenier's polar factorization problem [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF] implies that the squared distance between a map m : Ω → R and the set S of measure-preserving maps is equal to the squared Wasserstein distance [START_REF] Villani | Optimal transport: old and new[END_REF] between the restriction of the Lebesgue measure to Ω, denoted Leb, and its pushforward m # Leb under the map m:

d 2 S (m) = min s∈S m -s 2 = W 2 2 (m # Leb, Leb).
Moreover, since m is piecewise-constant over the partition (ω i ) 1≤i≤N , the push-forward measure m # Leb if finitely supported. Denoting by M i ∈ R d the constant value of the map m on the subdomain ω i we have

m # Leb = 1≤i≤N Leb(ω i )δ M i = 1 N 1≤i≤N δ M i .
Thus, computing the projection operator P S amounts to the numerical resolution of an optimal transport problem between the Lebesgue measure on Ω and a finitely supported measure. Thanks to recent work [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF][START_REF] Mérigot | A multiscale approach to optimal transport[END_REF][START_REF] De Goes | Blue noise through optimal transport[END_REF][START_REF] Lévy | A numerical algorithm for L 2 semi-discrete optimal transport in 3d[END_REF], this problem can be solved efficiently in dimensions d = 2, 3. We give more details in Section 5.

Remark 1.5. The idea of using optimal transport to impose incompressibility contraints has recently been exploited as a heuristic for computational fluid dynamics simulations in computer graphics [START_REF] De Goes | Power particles: an incompressible fluid solver based on power diagrams[END_REF]. From the simulations presented in [START_REF] De Goes | Power particles: an incompressible fluid solver based on power diagrams[END_REF], it seems that the scheme behaves better numerically, and it also has the extra advantage of not depending on a penalization parameter ε. However, it comes with no mathematical convergence analysis, and even its (formal) consistence is not obvious. It would therefore be interesting to extend the convergence analysis presented in Theorem 1.4 to the scheme presented in [START_REF] De Goes | Power particles: an incompressible fluid solver based on power diagrams[END_REF]. This however probably requires new ideas, as our technique of proof relies heavily on the fact that the space-discretization is hamiltonian, an assumption which does not seem to hold for the discretization of [START_REF] De Goes | Power particles: an incompressible fluid solver based on power diagrams[END_REF].

Remark 1.6. Our discretization (1.4) resembles (and derives from) a space-discretization of Euler's equations (1.1) introduced by Brenier in [START_REF] Brenier | Derivation of the Euler equations from a caricature of Coulomb interaction[END_REF]. The domain is also decomposed into subdomains (ω i ) 1≤i≤N , and one considers the set S N ⊆ S, which consists of measurepreserving maps s : Ω → Ω that are induced by a permutation of the subdomains. Equivalently, one requires that there exists a permutation s : {1, . . . , N } → {1, . . . , N } such that s(ω i ) = ω s(j) . The space-discretization considered in [START_REF] Brenier | Derivation of the Euler equations from a caricature of Coulomb interaction[END_REF] leads to an ODE similar to (1.4), but where the squared distance to S is replaced by the squared distance to S N . This choice of discretization imposes strong contraints on the relative size of the parameters τ , h N and , namely that h N = O(ε 8 ) and τ = O(ε 4 ). Such constraints still exist with the discretization that we consider here, but they are milder. In Theorem 1.4 the condition τ = o( 2) is due to the time discretization of (1.6) and can be improved using a scheme more accurate on the conservation of the Hamiltonian (1.5). However even with an exact time discretization of the Hamiltonian, the condition τ = o( ) remains mandatory, as explained at the end of Section 4.

Preliminary discussion on geodesics

To illustrate the approximate geodesic scheme we focus on the very simple example of R seen as R × {0} ⊂ R 2 . The geodesic is given by the function γ:

[0, T ] → R 2 with      γ(t) = (t, 0), t ∈ [0, T ], γ(0) = (0, 0), γ(0) = (1, 0). (2.1)
As in (1.4) we consider the solutions of the Hamiltonian system associated to:

H(m, ṁ) = 1 2 || ṁ|| 2 + 1 2 2 d 2 R×{0} (m). (2.2) That is      m(t) = 1 2 (P R (m) -m) = 1 2 2 ∇d 2 R×{0} (m), t ∈ [0, T ], m(0) = (0, h 0 ), ṁ(0) = (1, h 1 ). (2.3)
where P R (m) is the orthogonal projection from R 2 onto R × {0}. Notice that we assumed an initial error of h 0 on the initial position and h 1 on the initial velocity. In this case the solution is explicit and reads

m(t) = t, h 0 cos t + h 1 sin t . (2.4)
A convenient way to quantify how far m is from being a geodesic is to use a modulated energy related to the Hamiltonian H and the solution γ. We define E γ by

E γ (t) = 1 2 || ṁ(t) -γ(t)|| 2 + 1 2 2 d 2 R×{0} (m(t)). ( 2

.5)

A direct computation leads to

E γ (t) = h 2 0 2 + h 2 1 . (2.6)
This estimate shows that the velocity vector field ṁ converges towards the geodesic velocity vector fields γ as soon as h 0 goes to 0 faster then . Our construction of approximate geodesics for the Euler equation follow this idea. Estimates (2.6) suggests that our convergence results for the incompressible Euler equation in Theorem 1.2 is sharp.

Convergence of the approximate geodesics model

In this section we prove Theorem 1.2.

3.1. Strategy of the proof. We use a modulated energy approach. Let v be a solution of (1.1) and m a solution of (1.4) and for any t ∈ [0, T ], denote σ(t) = P S (m(t)). In other words, σ(t) is an arbitrary choice of a projection of m(t) on S. Equation (1.4) is the ODE associated to the Hamiltonian H :

M N × M N → R H(m, ṁ) = 1 2 ṁ 2 M + d 2 S (m) 2 2 .
We therefore consider a energy involving this Hamiltonian, modulated with the exact solution v:

E v (t) = 1 2 ṁ(t) -v(t, m(t)) 2 M + d 2 S (m) 2 2 . (3.1)
The core of the proof is to obtain a control on E v using a Gronwall estimate. As a first step we collect some lemmas. Lemmas 3.1 and 3.2 concern the projections Π M N and Π S and their orthogonality properties. Lemma 3.3 is necessary to ensure that the modulated energy introduced in (3.1) is well defined (the difficulty is that there is no reason that m(t, Ω) ⊆ Ω, and it is therefore necessary to extend v outside of Ω). Then we compute the derivative of (3.1) and modify its expression so as to identify terms of quadratic order, which are easier to control. This leads us to (3.7), which expresses the derivative of (3.1) as a sum of many terms. Each term is then estimated to obtain a Gronwall control. we keep track of the constants all along the proof.

3.2. Preliminary lemma. Before proving Theorem 1.2, we collect a few useful lemmas.

As before, Ω is a bounded and connected domain of R d with Lipschitz boundary.

Lemma 3.1 (Projection onto the measure preserving maps S).

Let m ∈ M = L 2 (Ω, R d ).
There exists a convex function ϕ : Ω → R, which is unique up to an additive constant, such that σ ∈ M belongs to Π S (m) if and only if m = ∇ϕ • σ up to a negligible set. Moreover, m -σ is orthogonal to the space

H div (Ω) • σ, that is ∀v ∈ H div (Ω), Ω m(x) -σ(x)|v(σ(x)) dx = 0. (3.2)
Proof. The first part of the statement is Brenier's polar factorization theorem [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF]. We first remark that

d 2 S (m) = inf s∈S m(x) -s(x) 2 dx ≥ inf π∈Π(m # Leb,Leb) x-y 2 dπ(x, y) = W 2 2 (m # Leb, Leb).
To prove the reverse inequality let ∇ϕ be the optimal transport map between m # Leb = 1≤i≤N δ M i and Leb. Let

L i = ∇ϕ -1 (M i ), by construction Leb(L i ) = 1 N .
For any i ∈ {1...N } let σ i be a measure preserving map between ω i and L i , we define a measure preserving map σ ∈ S by σ |ω i = σ i (anything can be done on the boundaries of the cells).

By construction m = ∇ϕ • σ and W 2 2 (m # Leb, Leb) = m -σ 2 .
The uniqueness of ϕ follows from the connectedness of the domain. Using a regularization argument we deduce the orthogonality relation

Ω m(x)|v(σ(x)) dx = Ω ∇ϕ • σ(x)|v(σ(x)) dx = Ω ∇ϕ(x)|v(x) = - Ω ϕdiv v(x) = 0.
Π M N (g) = N i=1 G i 1 ω i , with G i := 1 Leb(ω i ) ω i g(x)dx
and where 1 ω i is the indicator function of the subdomain ω i .

Proof. It suffices to remark that for any

m ∈ M N , m = 1≤i≤N M i 1 ω i , g|m M = Ω m(x)|g(x) dx = 1≤i≤N M i | ω i g(x)dx = m| i G i 1 ω i M Lemma 3.3. Let Ω ⊂ R d , let (V, .
) be a finite-dimensional normed vector space. There exists a linear map L :

C 1,1 (Ω, V ) → C 1,1 (R d , V ) such that for any f ∈ C 1,1 (Ω, V ), (i) Lf | Ω = f , (ii) Lf C 1,1 (R d ,V ) ≤ C Lf C 1,1 (Ω,V ) .
Proof. This lemma is a particular case of Theorem 2 in [START_REF] Fefferman | Extension of c {m,Ω} -smooth functions by linear operators[END_REF]. We also refer to [START_REF] Brudnyi | Whitney's extension problem for multivariate c {1,ω} -functions[END_REF][START_REF] Fefferman | Whitney's extension problem for c m[END_REF] for previous results.

We are now ready to prove Theorem 1.2. In the following the dot refers to the time derivative and .|. to the Hilbert scalar product on M. By abuse of notation we denote by the same name a C 1,1 function defined on Ω and its (also C 1,1 ) extension defined on the whole space R d using Lemma 3.3. The space R d is equipped with the canonical Euclidian norm, and the space of d × d matrices are equiped with the induced dual norm. All the Lipschitz constants that we consider are with respect to these two norms. Finally for a curve γ :

t ∈ [0, T ] → γ(t, •) we denote Lip [0,T] (γ) = sup t∈[0,T ] Lip(γ(t, •)). Material derivatives. Given (v, p) ∈ C 1 ([0, T ], C 1,1 (R d , R d )) × C 1 ([0, T ], C 1,1 (R d , R d ))
and X ∈ M, we define the two following functions, often called material derivatives:

D t v(t, X) = ∂ t v(t, X) + (v(t, X) • ∇) v(t, X), D t p(t, X) = ∂ t p(t, X) + v(t, X), ∇p(t, X)) . (3.3)
Remark that D t v and D t p are Lipschitz operators with

Lip [0,T] (D t v) ≤ Lip [0,T] (∂ t v) + Lip [0,T] (v) ∇v L ∞ + Lip [0,T] (∇v) v L ∞ ≤ Lip [0,T] (∂ t v) + Lip [0,T] (v) Lip [0,T] (v) + Lip [0,T] (∇v) v L ∞ (3.4) Lip [0,T] (D t p) ≤ Lip [0,T] (∂ t p) + Lip [0,T] (v) ∇p L ∞ + Lip [0,T] (∇p) v L ∞ ≤ Lip [0,T] (∂ t p) + Lip [0,T] (v) Lip [0,T] (p) + Lip [0,T] (∇p) v L ∞ . ( 3
d dt E v (t) = -m(t), v(t, m(t)) I 1 -ṁ(t) -v(t, m(t)), ∂ t v(t, m(t)) + ( ṁ(t) • ∇) v(t, m(t)) . I 2
(3.6) Using the EDO (1.4), I 1 can be rewritten as

2 I 1 = m(t) -P M N (σ(t)), v(t, m(t)) = m(t) -σ(t), v(t, m(t)) + σ(t) -P M N (σ(t)), v(t, m(t)) = m(t) -σ(t), v(t, m(t)) -v(t, σ(t)) 2 I 3
, where we have used that σ(t) -P M N (σ(t)) is orthogonal to M N and that m(t) -σ(t) is orthogonal to H div (Ω) • σ, see Lemmas 3.2 and 3.1. To handle the term I 2 we use the material derivatives defined by (3.3). Remark that Euler equations (1.1) implies that D t v(t, σ(t)) = -∇p(t, σ(t)). This leads to

I 2 = -ṁ(t) -v(t, m(t)), ∂ t v(t, m(t)) + (v(t, m(t)) • ∇) v(t, m(t)) -ṁ(t) -v(t, m(t)), ( ṁ(t) -v(t, m(t)) • ∇) v(t, m(t)) I 4 = I 4 -ṁ(t) -v(t, m(t)), D t v(t, m(t)) -D t v(t, σ(t)) I 5 + ṁ(t) -v(t, m(t)), ∇p(t, σ(t)) I 6
We rewrite I 6 as

I 6 = -ṁ(t) -v(t, m(t)), ∇p(t, m(t)) -∇p(t, σ(t)) I 7 + ṁ(t) -v(t, m(t)), ∇p(t, m(t)) = I 7 + d dt Ω p(t, m(t, x)))dx -J(t)
-

Ω ∂ t p(t, m(t, x)) -v(t,
I 3 ≤ m(t) -σ(t), v(t, m(t)) -v(t, σ(t)) 2 ≤ Lip(v(t)) m(t) -σ(t) 2 M 2 ≤ Lip [0,T] (v)E v (t). (3.8)
Furthermore

I 4 ≤ sup x∈R d ||∇v(t, x)|| ṁ(t) -v(t, m(t)) 2 M ≤ Lip [0,T] (v)E v (t), (3.9)
Where C depends only on the dimension d. To estimate I 5 and later I 8 we use that D t v and D t p are Lipschitz operators with constants given by (3.4) and (3.5). For I 5 we obtain

I 5 ≤ | ṁ(t) -v(t, m(t)), D t v(t, m(t)) -D t v(t, σ(t)) | ≤ Lip [0,T] (D t v) ṁ(t) -v(t, m(t)) M m(t) -σ(t) M ≤ Lip [0,T] (D t v)E v (t), (3.10)
where we used d S (m(t)) = m(t) -σ(t) M ≤ E v (t) and ṁ(t)-v(t, m(t)) M ≤ E v (t) to get from the second to the third line. The quantity I 7 can be bounded likewise:

I 7 ≤ | ṁ(t) -v(t, m(t)), ∇p(t, m(t)) -∇p(t, σ(t)) | ≤ Lip [0,T] (∇p)E v (t).
(3.11)

Finally to estimate I 8 and J we can assume that Ω p(t, x)dx = 0 since the pressure is defined up to a constant. Using that σ(t) is measure-preserving, this gives

Ω D t p(t, σ(t, x))dx = Ω ∂ t p(t, σ(t, x)) + v(t, σ(t, x)), ∇p(t, σ(t, x))) dx = Ω ∂ t p(t, x))dx + Ω v(t, x), ∇p(t, x)) dx = 0,
Therefore, using Young's inequality,

I 8 ≤ Ω D t p(t, m(t, x))dx - Ω D t p(t, σ(t, x))dx ≤ Lip [0,T] (D t p) m(t) -σ(t) L 1 (Ω) ≤ 1 2 ||m(t) -σ(t)|| 2 L 2 (Ω) 2 2 + C(Ω) Lip [0,T] (D t p) 2 ≤ 1 2 E v (t) + C(Ω) Lip [0,T] (D t p) 2 , (3.12)
where in this estimates and in the following estimates C(Ω) is a constant depending only on the Lebesgue measure of Ω. Similarly,

|J(t)| ≤ Ω p(t, m(t, x))) -p(t, σ(t, x)))dx ≤ Lip [0,T] (p)||m(t) -σ(t)|| L 1 (Ω) ≤ 1 2 E v (t) + C(Ω) Lip [0,T] (p) 2 . (3.13)
We finally remark that

|J(0)| ≤ Lip [0,T] (p)h N . (3.14)
Remark 3.5. The last two estimates show that we can add d dt J into the Gronwall argument. It is a general fact that the derivative of a controlled quantity can be added. This is a classical way of controlling the term of order one in the energy.

(D t v) + Lip [0,T] (∇p) + 1 2 E v (t) + C(Ω) Lip [0,T] (D t p) 2
Remark that (3.13) implies that for any K > 0,

KE v (t) ≤ KE v (t)+2KJ(t)-2KJ(t) ≤ 2KE v (t)+2KJ(t)+2KC(Ω) Lip [0,T] (p) 2 . (3.15)
Therefore, setting

C 1 = C(Ω) 4 Lip [0,T] (v) + 2 Lip [0,T] (D t v) + 2 Lip [0,T] (∇p) + 1 , C 2 = C(Ω) Lip [0,T] (D t p) + C 1 Lip [0,T] (p) , we obtain d dt (E v (t) + J(t)) ≤ C 1 (E v (t) + J(t)) + C 2 2 .
We deduce from the Gronwall inequality that for any t ∈ [0, T ]:

E v (t) ≤ (E v (0) + J(0)) + C 2 T 2 e C 1 T -J(t).
Using the estimation (3.13) one more time we obtain

E v (t) ≤ 2 E v (0) + Lip [0,T] (p)h N + C 2 T 2 e C 1 T + C(Ω) Lip [0,T] (p) 2 .
Finally, using that

E v (0) = 1 2 P M (v 0 ) -v 0 2 M + d 2 S (Id ) 2 2 ≤ h 2 N 2 + h 2 N we obtain ṁ(t) -v(t, m(t)) 2 M ≤ 2E v (t) ≤ 2 2 h 2 N 2 + h 2 N 2 2 + Lip [0,T] (p)h N + C 2 T 2 e C 1 T + C(Ω) Lip [0,T] (p) 2 (3.16) ≤ C 1 h 2 N 2 + C 2 2 + C 3 h N (3.17) where        C 1 = 2e C 1 T C 2 = C 2 T e C 1 T + C(Ω) Lip [0,T] (p) C 3 = (1 + Lip [0,T] (p))e C 1 T .
In order to estimate ṁ(t) -v(t, φ(t)) 2 M we need one additional Gronwall estimate:

ṁ(t) -v(t, φ(t)) 2 M ≤ 2 ṁ(t) -v(t, m(t)) 2 M + 2 v(t, m(t)) -v(t, φ(t)) 2 M ≤ 2E v (t) + 2(Lip [0,T] (v)) 2 m(t) -φ(t) 2 M ≤ 2E v (t) + 4(Lip [0,T] (v)) 2 m(0) -φ(0) 2 M + 4(Lip [0,T] (v)) 2 t 0 ( ṁ(s) -φ(s))ds 2 M ≤ 2E v (t) + 4(Lip [0,T] (v)) 2 h 2 N + 4T (Lip [0,T] (v)) 2 t 0 ( ṁ(s) -φ(s)) 2 M ds ≤ C 1 h 2 N 2 + C 2 2 + 4(Lip [0,T] (v)) 2 h 2 N + C 3 h N + 4T (Lip [0,T] (v)) 2 t 0 ṁ(s) -φ(s) 2 M ds (3.18)
where we used Jensen's inequality to obtain the second to last line. We conclude thanks to Gronwall inequality:

ṁ(t) -v(t, φ(t)) 2 M ≤ C 1 h 2 N 2 + C 2 2 + 4((Lip [0,T] (v)) 2 h N + C 3 )h N e 4T (Lip [0,T] (v)) 2 T ≤ C 1 h 2 N 2 + C 2 2 + C 3 h N . (3.19)
We used that and h N are smaller than C(Ω) for (3.17) and (3.19). Observe that the right-hand side of (3.17) and (3.19) goes to zero provided that h N and go to zero. This finishes the proof of Theorem 1.2.

Remark 3.6. Using (3.4) and (3.5), the constants C 1 , C 2 are bounded by:

     1 C(Ω) C 1 ≤ 1 + 4 Lip [0,T] (v) + 2 Lip [0,T] (∇p) +2 Lip [0,T] (∂ t v) + (Lip [0,T] (v)) 2 + Lip [0,T] (∇v) v L ∞ , 1 C(Ω) C 2 ≤ Lip [0,T] (p) + C 1 Lip [0,T] (∂ t p) + Lip [0,T] (v) Lip [0,T] (p) + Lip [0,T] (∇p) v L ∞ .
A close look to the explicit value of the constants C 1 , C 2 and C 1 , C 2 , C 3 , together with a diagonal argument shows that our scheme can be used to approximate solutions less regular than those supposed in Theorem 1.2. For example, it is possible to establish the following theorem: Let v, p be a solution of Euler's equation (1.1), where v is merely Lipschitz in space but where there exists (v k , p k ) k∈N a sequence of regular (in the sense of Theorem 1.2) solutions of (1.1) such that v k (0, •) -→ v(0, •) in M and Lip T (v k ) -→ Lip T (v). Then there exists N k and k , depending polynomially on the data such that ṁk (t) -v (t, m k (t)) 2

M

goes to zero as k goes to infinity, where m k is the solution of (1.6) with initial conditions m k (0) = P M N k (id) and ṁk (0) = P M N k (v k (0)) and with parameter = k . If one allows an exponential dependence on the data, it is possible to approach any solution whose velocity v belongs to the L 2 closure of the regular solutions to Euler's equation.

Convergence of the symplectic Euler scheme

In this section we prove a statement which is slightly more general than Theorem 1.4 (see Remark 4.3), and which allows a sort of a posteriori estimates. The proof follows the proof of Theorem 1.2, but one has to deal with some additional term coming from the time discretization. It combines two Gronwall estimates. The first one is a continuous Gronwall argument on each segment [nτ, (n+1)τ ], and the second one is a discrete Gronwall estimate comparing a timestep to the next one. Both steps rely on the same modulated energy.

M 0 = P M N (id), V 0 = P M N (v(0, •)). Finally let H n = H(M n , V n ) = 1 2 V n 2 M + d 2 S (M n ) 2 2 , ( 4.1
)

and κ = max n∈N∩[0,T /τ ] H n -H 0 .
Then, assuming τ ≤ and h N ≤ ε, we have

max n∈N∩[0,T /τ ] V n -v(t n , φ(t n , •)) M ≤ C 2 + h N + h 2 N 2 + τ + κ ,
where the constant C only depends on Ω, on the L ∞ norm (in space) of the velocity v(t, •) and on the Lipschitz norms (in space) of the velocity and its first derivatives ∇v(t, •), ∂ t v(t, •) and of the pressure and its derivatives p(t, •), ∇p(t, •), ∂ t p(t, •).

4.1.

Preliminary lemma. Given a solution of (1.11) and s ∈ [0, 1] and n ∈ N, we denote the linear interpolates between two timesteps nτ and (n + 1)τ by:

V n+s = V n -sτ M n -P M N •P S (M n ) 2 M n+s = M n + sτ V n+1 , (4.2)
We consider the Hamiltonian H n+s and modulated energy E n+s defined by

H n+s = 1 2 V n+s 2 M + d 2 S (M n+s ) 2 2 , E n+s = 1 2 V n+s -v ((n + s)τ, M n+s ) 2 M + d 2 S (M n+s ) 2 2 . (4.3)
We start with a lemma quantifying the conservation of the Hamiltonian.

∈ N ∩ [0, T /τ ], 1 - τ 2 2 H n+1 ≤ H n , (4.4)
H n ≤ Ce T τ -2 , (4.5)

H n+s ≤ H n + τ 2 2 H n+1 , (4.6)
Proof. The proof is based on the 1-semiconcavity of 1 2 d 2 S , see Proposition 5.2 for details. On the one hand the 1-semiconcavity of 1 2

d 2 S reads d 2 S (M n+s ) 2 2 ≤ d 2 S (M n ) 2 2 + sτ V n+1 , M n -P M N • P S (M n ) 2 + s 2 τ 2 2 2 V n+1 2 M ,
where we used that M n -P M N • P S (M n ) belongs to the superdifferential of the function d 2 S at M n , and the definition of the scheme (4.2). On the other hand, (4.2) again, leads to

V n+s 2 M 2 = V n 2 M 2 -sτ V n , M n -P M N • P S (M n ) 2 + s 2 τ 2 M n -P M N • P S (M n ) 2 2 M
Summing both equations and using (4.2) gives

H n+s ≤ H n + τ 2 s(s -1) 2 M n -P M N • P S (M n ) 2 M 2 + s 2 τ 2 2 V n+1 2 M 2 (4.7)
Taking s = 1 in (4.7) proves (4.4). The inequality (4.5) is a direct consequence of (4.4), while (4.6) follows from the combination of (4.4) and (4.7).

κ ≤ n∈N∩[0,T /τ ]0 H n+1 -H n ≤ τ 2 T e T τ -2 1 2 V 0 2 M + h 2 N 2 2 .
Using this upper bound Theorem 4.1 becomes Theorem 1.4 and the condition κ = o(1) becomes τ = o( 2). However numerically one can expect some compensation in H n and thus obtain a better "a posteriori bound" for κ in order to get rid of the strong assumption τ = o( 2 ). Figure 5.4 illustrates the conservation of the Hamiltonian in two test cases.

Notice that this estimate is not a posteriori in the usual sense since the constants in Theorem 4.1 also depend on the unknown limiting solution. The condition τ = o( ) seems mandatory for the proof techniques to work.

4.2.

The modulated energy. Remark that with the definitions of the Hamiltonian and modulated energy, we have

E n+s = H n+s -V n+s , v (n + s)τ, M n+s + 1 2 v (n + s)τ, M n+s 2 M , (4.8)
so that for any s ∈ [0, 1] and any n ∈ N,

E n+s = E n + H n+s -H n + s 0 d n+θ dθ, (4.9)
where

d n+s = d ds -V n+s , v (n + s)τ, M n+s + 1 2 v (n + s)τ, M n+s 2 M .
To evaluate d n+s , we introduce σ p = P S (M p ) and we will use the compact notation

v n+s M p = v((n + s)τ, M p ), ∂ t v n+s M p = ∂ t v((n + s)τ, M p ), ∇v n+s M p = ∇v((n + s)τ, M p ), v n+s σ p = v((n + s)τ, σ p ), ∂ t v n+s σ p = ∂ t v((n + s)τ, σ p ), ∇v n+s σ p = ∇v((n + s)τ, σ p
). We will also use a similar notation for the material derivative of the velocity and for the pressure and its derivatives. Remark 4.4. As before, the main idea of the following computation is to try to find terms of quadratic order in the expression. To control the remaining linear term we have to rewrite it as a derivative of a small quantity and add it in the Gronwall argument.

d n+s = - d ds V n+s , v n+s M n+s -V n+s , d ds v n+s M n+s + v n+s M n+s , d ds v n+s M n+s = τ -2 M n -P M N • P S (M n ), v n+s M n+s I 1 -V n+s -v n+s M n+s , τ ∂ t v n+s M n+s + d ds M n+s • ∇v n+s M n+s I 2
Recalling that σ n = P S (M n ), the term I 1 can be rewritten as

I 1 = τ -2 M n -P M N (σ n ), v n+s M n+s = τ -2 M n -σ n , v n+s M n+s + σ n -P M N (σ n ), v n+s M n+s = τ -2 M n -σ n , v n+s M n+s -v n+s σ n I 3
Here we had to control the fact that, due to the double projection, the norm of the acceleration M n -P M N • P S (M n ) 2 M is not equal to the squared distance d 2 S (M n ). We used the orthogonality property of the double projection for that purpose. On the one hand σ n -P M N (σ n ) is orthogonal to M N since M N is a linear subspace of M. On the other hand M n -σ n is orthogonal to the tangent space to S at σ n , see Lemma 3.1.

To handle the term I 2 we use the material derivatives defined in (3.3),

I 2 = -V n+s -v n+s M n+s , τ ∂ t v n+s M n+s + d ds M n+s • ∇v n+s M n+s I 2 = -V n+s -v n+s M n+s , τ ∂ t v n+s M n+s + τ v n+s M n+s • ∇v n+s M n+s -V n+s -v n+s M n+s , d ds M n+s -τ v n+s M n+s • ∇v n+s M n+s I 4 = I 4 -τ V n+s -v n+s M n+s , D t v n+s M n+s -D t v n+s σ n+s I 5 + τ V n+s -v n+s M n+s , ∇p n+s σ n+s I 6
.

We used that D t v n+s σ n+s = -∇p n+s σ n+s . We now rewrite I 6 using d ds M n+s = τ V n+1 :

I 6 = τ V n+s -v n+s M n+s , ∇p n+s σ n+s -∇p n+s M n+s I 7 +τ V n+s -v n+s M n+s , ∇p n+s M n+s = I 7 + d ds M n+s , ∇p n+s M n+s + τ V n+s -V n+1 , ∇p n+s M n+s -τ v n+s M n+s , ∇p n+s M n+s = I 7 + d ds Ω p n+s M n+s dx -J n+s -τ Ω (∂ t p n+s M n+s + v n+s M n+s , ∇p n+s M n+s )dx + (1 -s)τ 2 -2 M n -P M N • P S (M n ), ∇p n+s M n+s I 8 = I 7 + I 8 - d ds J n+s -τ Ω D t p n+s M n+s dx I 9 ,
We need to estimate all the terms in the following formula.

τ -1 I 7 ≤ V n+s -v n+s M n+s , ∇p n+s σ n+s -∇p n+s M n+s ≤ CE n+s (4.14)
For the estimation of I 8 we use σ n -P M N (σ n ), ∇p n+s M n+s = 0 to get

I 8 = (1 -s)τ 2 -2 M n -P M N • P S (M n ), ∇p n+s M n+s ≤ τ 2 -2 ||∇p((n + s)τ )|| L ∞ (Ω) ||M n -σ n || M ≤ τ 2 -1 Lip [0,T] (v) ||M n -σ n || M ≤ τ CE n + τ 2 -1 C. (4.15)
To estimate J and I 9 recall that we have assumed that Ω p(t, x)dx = 0, which implies in particular that Ω D t p(t, σ n (t, x))dx = 0. Therefore,

τ -1 I 9 ≤ Lip [0,T] (D t p)||M n+s -σ n+s || L 1 (Ω) ≤ 1 2 E n+s + C 2 . ( 4
J n+s + s 0 d n+θ dθ ≤ J n + 2τ CE n + Cτ 2 -1 (H 0 + κ) + τ C s 0 E n+θ dθ + τ CE n + τ C s 0 E n+θ dθ + τ C s 0 E n+θ dθ + τ CE n + τ 2 -1 C + τ 2 s 0 E n+θ dθ + τ 2 C ≤ J n + Cτ E n + Cτ 2 + C(H 0 + κ)τ 2 -1 (4.18) + τ C s 0 E n+θ + J n+θ dθ.
Remark that we used (3.15) to add J n+θ at the last line. Remark also that we only kept the first order terms using ≤ C. Plugging (4.18) into (4.9) we obtain

E n+s + J n+s ≤ α(s) + β θ 0 (E n+θ + J n+θ )ds) (4.19)
where α(s

) = E n + J n + H n+s -H n + Cτ E n + Cτ 2 + C(H 0 + κ)τ 2 -1 , β = τ C
so that by Gronwall lemma,

E n+1 + J n+1 ≤ α(1) + 1 0 α(s)β exp((1 -s)β)ds ≤ E n + J n + Cτ E n + Cτ 2 + C(H 0 + κ)τ 2 -1 e Cτ + H n+1 -H n + 1 0 H n+s -H n Cτ exp((1 -s)Cτ )ds R
Using Lemma 4.2 and in particular the upper bound (4.6) we find

R ≤ τ 2 2 H n+1 s 0 Cτ e Cτ (1-θ) ≤ C τ 2 2 (H 0 + κ) e Cτ -1
so that

E n+1 + J n+1 ≤ [(1 + Cτ ) (E n + J n ) (4.20) + Cτ 2 + C(H 0 + κ)τ 2 -1 +H n+s -H n + τ 2 -2 (H 0 + κ) e Cτ -1 e Cτ .
4.5. Discrete Gronwall step. From (4.20) and the discrete Gronwall inequality we deduce that for any n ∈ N ∩ [0, T /τ ],

E n + J n ≤ E 0 + J 0 + CT 2 + CT (H 0 + κ)τ -1 + H n -H 0 +τ 2 -2 (H 0 + κ) T τ e Cτ -1 (1 + Cτ ) n e CT ≤ C E 0 + J 0 + 2 + (H 0 + κ)τ -1 + κ +(H 0 + κ)τ 2 -2 e CT e CT ≤ C E 0 + J 0 + 2 + (H 0 + κ)τ -1 + κ e CT .
We used the mean value theorem to obtain the second to last line. Using (4.17) one last time and H 0 ≤ C leads us to

E n ≤ C E 0 + J 0 + 2 + τ -1 + κ + C 2 ≤ C 2 + h N + h 2 N 2 + κ + τ .
where the second line incorporates the initial error. It leads

max n∈N∩[0,T /τ ] V n -v(t n , M n ) 2 M ≤ C 2 + h N + h 2 N 2 + κ + τ .
A third Gronwall estimate, similar to the one done to obtain (3.19), concludes the proof:

max n∈N∩[0,T /τ ] V n -v(t n , φ(t n , •)) 2 M ≤ C 2 + h N + h 2 N 2 + κ + τ .. Remark 4.5.
A close look at the constant leads to a similar result as the one given in Remark 3.6: namely the convergence of the numerical scheme towards less regular solutions of the Euler's equations.

Remark 4.6. The method of the proof is robust and could easily be adapted to other numerical scheme. Any improvement to the estimate given in Lemma 4.2 (conservation of the Hamiltonian) will lead to improved convergence estimates for the numerical scheme.

Numerical implementation and experiments

5.1. Numerical implementation. We discuss here the implementation of the numerical scheme (1.11) and in particular the computation of the double projection P M N • P S (m) for a piecewise constant function m ∈ M N . Using Brenier's polar factorisation theorem, the projection of m on S amounts to the resolution of an optimal transport problem between Leb and the finitely supported measure m # Leb. Such optimal transport problems can be solved numerically using the notion of Laguerre diagram from computational geometry.

Lag i (M, ψ) = x ∈ R d | ∀j ∈ {1, . . . , N }, x -M i 2 + ψ i ≤ x -M j 2 + ψ j } .
In We denote L i := Lag i (M, ψ). Then, a function s ∈ S is a projection of m on S if and only if it maps the subdomain ω i to the Laguerre cell L i up to a negligible set, that is:

Π S (m) = {s ∈ S | ∀i ∈ {1, . . . , N }, Leb(s(ω i )∆L i ) = 0}, (5.2)
where A∆B denotes the symmetric difference between sets A and B. Moreover, the squared distance d 2 S is differentiable at m and, setting

B i = 1 Leb(L i ) L i xdx, one has d 2 S (m) = 1≤i≤N L i x -M i 2 dx, ∇d 2 S (m) = 2(m -P M N • P S (m)) with P M N • P S (m) = 1≤i≤N B i 1 ω i .
(5.3)

Proof. The existence of a vector (ψ i ) 1≤i≤N satisfying Equation (5.1) follows from optimal transport theory (see Section 5 in [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF] for a short proof), and its uniqueness follows from the connectedness of the domain Ω. In addition, the map T : Ω → {M 1 , . . . , M N } defined by T (L i ) = M i (up to a negligible set) is the gradient of a convex function and therefore a quadratic optimal transport between Leb and the measure1 N Leb(Ω) i δ M i . By Brenier's polar factorization theorem, summarized in Lemma 3.1,

s ∈ Π S (m) ⇐⇒ m = T • s a.e. ⇐⇒ ∀i ∈ {1, . . . , N }, Leb(ω i ∆(T • s) -1 ({M i })) = 0 ⇐⇒ ∀i ∈ {1, . . . , N }, Leb(s(ω i )∆L i ) = 0,
where the last equality holds because s is measure preserving. To prove the statement on the differentiability of d 2 S , we first note that the function

d 2 S is 1-semi-concave, since D(m) := m 2 -d 2 S (m) = m 2 -min s∈S m -s 2 = max s∈S 2 m|s -s 2 is convex. The subdifferential of D at m is given by ∂D(m) = {P M N (s) | s ∈ Π S (m)}, so that D (and hence d 2 S
) is differentiable at m if and only if P M N (Π S (m)) is a singleton. Now, note from Lemma 3.2 that for s ∈ Π S (m)

P M N (s) = 1≤i≤N bary(s(ω i ))1 ω i = 1≤i≤N bary(L i )1 ω i .
This shows that P M N (Π S (m)) is a singleton, and therefore establishes the differentiability of d 2 S at m, together with the desired formula for the gradient. The main difficulty to implement the numerical scheme (1.11) is the resolution of the discrete optimal transport problem (5.1), a non-linear system of equations which must be solved at every iteration. We resort to the damped Newton's algorithm presented in [START_REF] Kitagawa | Convergence of a newton algorithm for semi-discrete optimal transport[END_REF] (see also [START_REF] Mirebeau | Discretization of the 3d monge-ampere operator, between wide stencils and power diagrams[END_REF]) and more precisely on its implementation in the PyMongeAmpere library 1 . 5.1.1. Construction of the fixed tessellation of the domain. The fixed tessellation (ω i ) 1≤i≤N of the domain Ω is a collection of Laguerre cells that are computed through a simple fixedpoint algorithm similar to the one presented in [START_REF] De Goes | Blue noise through optimal transport[END_REF]. We start from a random sampling

(C 0 i ) 1≤i≤N of Ω. At a given step k ≥ 0, we compute (ψ i ) 1≤i≤N ∈ R N such that ∀i ∈ {1, . . . , N }, Leb(Lag i (C, ψ)) = 1 N Leb(Ω),
and we then update the new position of the centers (C k+1 i ) by setting C k+1 i := bary(Lag i (C k , ψ)). After a few iterations, a fixed-point is reached and we set ω i := Lag i (C k , ψ).

Iterations.

To implement the symplectic Euler scheme for (1.6), we start with M 0 i := bary(ω i ) and V 0 i := v 0 (M 0 i ). Then, at every iteration k ≥ 0, we use Algorithm 1 in [START_REF] Kitagawa | Convergence of a newton algorithm for semi-discrete optimal transport[END_REF] to compute a solution ( v 0 (x 1 , x 2 ) = (-cos(πx 1 ) sin(πx 2 ), sin(πx 1 ) cos(πx 2 ))

ψ k i ) 1≤i≤N ∈ R N to Equation (5.1) with M = M k , i.e. such that ∀i ∈ {1, . . . , N }, Leb(Lag i (M k , ψ k )) = 1 N Leb(Ω).
In Figure 1, we display the computed numerical solution using a low number of particles (N = 900) in order to show the shape of the Laguerre cells associated to the solution.

v 0 (x 1 , x 2 ) = 0.5 if x 2 ≥ 0 1 if x 2 < 0
This speed profile corresponds to a stationnary but unstable solution to Euler's equation. If the subdomains (ω i ) 1≤i≤N are computed following §5.1.1, the perfect symmetry under horizontal translations is lost, and in Figure 2 we observe the formation of vortices whose radius increases with time. This experiment involves N = 200 000 particles, with parameters τ = 0.002 and ε = 0.005, and 2 000 timesteps. As displayed in Figure 2, the hamiltonian of the system is very well preserved despite the roughness of the solution. This behaviour shows that the estimate of Lemma 4.2 might be overly pessimistic, and requires further investigation.

ρ i V k+1 i = ρ i V k i + τ 1 ε 2 (bary(Lag i (M k , ψ k ) -M k i ) + ρ i G M k+1 i = M k i + τ V k+1 i (5.5)
The computational domain is the rectangle Ω = [-1, 1] × [-3, 3], and the initial distribution of particles is given by C i = bary(ω i ), where the partition (ω i ) 1≤i≤N is constructed according to §5.1.1. The fluid is composed of two phases, the heavy phase being on top of the light phase:

ρ i = 3 if C i2 > η cos(πC i1 ) 1 if C i2 ≤ η cos(πC i1 ) ,
where η = 0.2 in the experiment and where we denoted C i1 and C i2 the first and second coordinates of the point C i . Finally, we have set N = 50 000, ε = 0.002 and τ = 0.001 and we have run 2000 timesteps. The computation takes less than six hours on a single core of a regular laptop. Note that it does not seem straighforward to adapt the techniques used in the proofs of convergence presented here to this setting, where the force depends on the density of the particle. Our purpose with this test case is merely to show that the numerical scheme behaves reasonably well in more complex situations.

Software. The software developed for generating the results presented in this article is publicly available at https://github.com/mrgt/EulerLagrangianOT

Theorem 1 . 2 .

 12 Let Ω be a bounded domain of R d with Lipschitz boundary. Let v, p be a strong solution of Euler's equations (1.1), let φ be the flow map induced by v (see (1.2)) and assume that v, p, ∂ t v, ∂ t p, ∇v and ∇p are Lipschitz on Ω, uniformly on [0, T ]. Suppose in addition that there exists a C 1 curve m : [0, T ] → M N satisfying the initial conditions

Lemma 3 . 2 (

 32 Projection onto the piecewise constant set M N ). The projection of a function g ∈ L 2 (Ω, R d ) on M N is the following piecewise constant function :

3. 4 . 8 ≤ 2

 482 Gronwall argument. Collecting estimates (3.8), (3.9), (3.10), (3.11), (3.12), we get d dt (E v (t) + J(t)) ≤ I 3 + I 4 + I 5 + I 7 + I Lip [0,T] (v) + Lip [0,T]

Theorem 4 . 1 .

 41 Let Ω be a bounded domain with Lipschitz boundary and let , τ positive numbers and let N ∈ N. Let v, p be a strong solution of (1.1), and let φ be the flow map induced by v (see (1.2)). Assume that v, p, ∂ t v, ∂ t p, ∇v and ∇p are Lipschitz on Ω, uniformly on [0, T ]. Let (M n , V n ) n≥0 be a sequence generated by (1.11) with initial conditions

Lemma 4 . 2 (

 42 Conservation of the Hamiltonian). For any s ∈ [0, 1] and n

Remark 4 . 3 .

 43 Lemma 4.2 gives an upper bound for κ in Theorem 4.1 namely

d n+s = I 3 + I 4 + I 5 + I 7 + I 8 + I 9 3 .

 3457893 Gronwall estimates on [nτ, (n + 1)τ ]. From now and for clarity we do not track the constants anymore, and C will be a constant depending only on T , Ω, Lip [0,T] (v), Lip [0,T] (p), Lip [0,T] (∇p), Lip [0,T] (D t v) and Lip [0,T] (D t p). The value of the constant C can

Definition 5 . 1 (

 51 Laguerre diagram). Let M = (M 1 , . . . , M N ) ∈ (R d ) N and let ψ 1 , . . . , ψ N ∈ R. The Laguerre diagram is a decomposition of R d into convex polyhedra defined by

Proposition 5 . 2 .

 52 the following proposition, we denote Π S (m) = {s ∈ S | m -s = d S (m)}, and for a bounded subset A ⊆ R d with positive measure we set bary(A) := 1 Leb(A) A xdx. Let m ∈ M N \ D N and define M i = m(ω i ) ∈ R d . Assume that Ω is a bounded and connected domain of R d with Lipschitz boundary. Then, there exist scalars (ψ i ) 1≤i≤N , which are unique up to an additive constant, such that ∀i ∈ {1, . . . , N }, Leb(Lag i (M, ψ)) = 1 N Leb(Ω). (5.1)

Figure 1 .) 2 2 .

 122 Figure 1. (Top row) Beltrami flow in the square, with N = 900 particles, τ = 1/50 and ε = .1. The particles are colored depending on their initial position in the square. From left to right, we display the Laguerre cells and their barycenters at timesteps k = 0, 24 and 49. The partition (ω i ) 1≤i≤N is induced by a regular grid. (Bottom row) Same experiment, but where the partition (ω i ) 1≤i≤N is optimized using the algorithm described in §5.1.1.

5. 3 .

 3 Kelvin-Helmoltz instability. For this second test case, the domain is the rectangle Ω = [0, 2]×[-.5, .5] periodized in the first coordinate by making the identification (4, x 2 ) ∼ (0, x 2 ) for x 2 ∈ [-.5, .5]. The initial speed v 0 is discontinuous at x 2 = 0: the upper part of the domain has zero speed, and the bottom part has unit speed:

Figure 2 .

 2 Figure 2. Numerical illustration of the Kelvin-Helmotz instability on a rectangle with periodic conditions (in the horizontal coordinate) involving a discontinuous initial speed. The parameters are given in §5.4.

Figure 3 .Figure 4 .

 34 Figure 3. Numerical illustration of the Rayleigh-Taylor instability occuring when a heavy fluid (in green) is placed over a lighter fluid (in red) at timesteps n = 0, 200, 400, . . . , 2000. The parameters are given in §5.4.

  Time derivative. We compute d dt E v (t) and modify the expression in order to identify terms of quadratic order. Since the Hamiltonian H( ṁ(t), m(t)) is preserved, we find

.5) 3.3. Proof of Theorem 1.2. We can now go to the proof of Theorem 1.2. Note that we need to use Lemma 3.3 to define the modulated energy E v in (3.1) since the maps m(t, •) ∈ M N can send points outside of Ω when Ω is not convex.

3.3.1.

  n+s M n+s -p n+s σ n+s dx ≤ Lip [0,T] (p)||M n+s -σ n+s || L 1 (Ω)

	Similarly
	|J n+s | = |J((n + s)τ )| ≤
	≤	1 2	E n+s + C 2 .	(4.17)
	Note also that J 0 ≤ Lip [0,T] (p)h N ≤ Ch N still holds see (3.14).
	4.4. Gronwall argument on [nτ, (n + 1)τ ]. Collecting estimates (4.11), (4.12), (4.13),
	(4.14), (4.15), (4.16) and (4.17) and integrating equation (4.10) from 0 to s we obtain
				.16)

Ω p

https://github.com/mrgt/PyMongeAmpere

Acknowledgements

We would like to thank Yann Brenier who pointed out to us the reference [8] on which this article elaborates, and for several interesting discussions at various stages of this work. We also thank Pierre Bousquet who indicated the reference [16] to us.

T.O.G was supported by ANR grant ISOTACE (ANR-12-MONU-0013) and by the Fonds de la Recherche Scientifique -FNRS under Grant MIS F.4539.16. Q.M. is supported by ANR grant MAGA (ANR-16-CE40-0014). COMMUNICATED BY EITAN TADMOR