A Lagrangian Scheme à la Brenier for the Incompressible Euler Equations

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: (1.1) As noticed by Arnold [2], 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: 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 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 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 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 [6], 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 Ω [23].
Previous work: discretization of geodesics in S. The first numerical experiments to recover generalized minimizing geodesics have been performed by Brenier in 1D [9]. He also proposed a scheme to compute the solutions of the Cauchy problem (1.3) in [5]. 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 [8], 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 [3,20,13,18]. This alternative discretization has already been used successfully to compute minimizing geodesics between measure-preserving maps in [21], 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 [8]. 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 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: 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 : 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: 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 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 [8] and requires C 1,1 regularity assumptions on the solution to Euler's equations. See also [10,12] for related works.
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 [1,4,19], see also [10,Appendix].
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: 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τ .
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 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 [7] 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 [24] between the restriction of the Lebesgue measure to Ω, denoted Leb, and its pushforward m # Leb under the map m: 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 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 [3,20,13,18], 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 [14]. From the simulations presented in [14], 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 [14]. 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 [14].
Remark 1.6. Our discretization (1.4) resembles (and derives from) a space-discretization of Euler's equations (1.1) introduced by Brenier in [8]. 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 : . The space-discretization considered in [8] 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 γ: As in (1.4) we consider the solutions of the Hamiltonian system associated to: 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 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 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 : We therefore consider a energy involving this Hamiltonian, modulated with the exact solution v: 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. Proof. The first part of the statement is Brenier's polar factorization theorem [7]. We first remark that 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 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 , ) be a finite-dimensional normed vector space. There exists a linear map L : such that for any f ∈ C 1,1 (Ω, V ), Proof. This lemma is a particular case of Theorem 2 in [16]. We also refer to [11,15] 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 γ : ) and X ∈ M, we define the two following functions, often called material derivatives: Remark that D t v and D t p are Lipschitz operators with 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.

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 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 We rewrite I 6 as Remark 3.4. The quantity I 5 + I 7 would vanish if (v, p) was a solution to the Euler equations on the whole space R d . This is not the case in our setting, as the couple (v, p) is constructed by the extension Lemma 3.3.
Collecting the above decompositions (3.6) rewrites (3.8) Furthermore 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 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: (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, where in this estimates and in the following estimates C(Ω) is a constant depending only on the Lebesgue measure of Ω. Similarly, We finally remark that (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. Remark that (3.13) implies that for any K > 0, Therefore, setting We deduce from the Gronwall inequality that for any t ∈ [0, T ]: Using the estimation (3.13) one more time we obtain Finally, using that In order to estimate ṁ(t) − v(t, φ(t)) 2 M we need one additional Gronwall estimate: where we used Jensen's inequality to obtain the second to last line. We conclude thanks to Gronwall inequality: 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: 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. 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.  (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

Finally let
and κ = max Then, assuming τ ≤ and h N ≤ ε, we have max 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, ·).

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: We consider the Hamiltonian H n+s and modulated energy E n+s defined by We start with a lemma quantifying the conservation of the Hamiltonian. H n ≤ Ce T τ −2 , (4.5)

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 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 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 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, 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 where the second line incorporates the initial error. It leads max A third Gronwall estimate, similar to the one done to obtain (3.19), concludes the proof: 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.
In 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 Proof. The existence of a vector (ψ i ) 1≤i≤N satisfying Equation (5.1) follows from optimal transport theory (see Section 5 in [3] for a short proof), and its uniqueness follows from the connectedness of the domain Ω. In addition, the map T : 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 [17] (see also [22]) 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 [13]. 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  v 0 (x 1 , x 2 ) = (− cos(πx 1 ) sin(πx 2 ), sin(πx 1 ) cos(πx 2 )) 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. 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.
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