A projection method for the computation of admissible measure valued solutions of the incompressible Euler equations

We formulate a fully discrete finite difference numerical method to approximate the incompressible Euler equations and prove that the sequence generated by the scheme converges to an admissible measure valued solution. The scheme combines an energy conservative flux with a velocity-projection temporal splitting in order to efficiently decouple the advection from the pressure gradient. With the use of robust Monte Carlo approximations, statistical quantities of the approximate solution can be computed. We present numerical results that agree with the theoretical findings obtained for the scheme.

1. Introduction. We consider the incompressible Euler equations, which model the motion of an inviscid, Newtonian fluid. These equations can be written as where u denotes the velocity field and p is the scalar pressure, which can be viewed as a Lagrange multiplier for the divergence-free constraint (1b). Incompressible Euler equations are a fundamental building block of fluid dynamics, and are considered a good model for flows with very low Mach number and high Reynolds number [23]. Yet, their mathematical and numerical understanding is still incomplete.
1.1. Known results. In the following, we summarise the main results about Euler equations.
Classical solutions, which require sufficient smoothness of the velocity, fail to exist in general cases, for example in the three dimensional torus. For special cases (e.g. in two dimensions, with periodic boundary conditions) classical solutions exist globally and are unique [23,2].
In general, one considers weak solutions to the Euler equations. Given a function u 0 ∈ L 2 x (D), a function u(t) ∈ L 2 x (D) is said to be a weak solution to the Euler 942 FILIPPO LEONARDI equations, if it satisfies the equations in the sense of distributions, namely: for all test functions ϕ ∈ C ∞ 0 ([0, T ) × D) d , with divϕ = 0. The pressure can be recovered from the velocity up to a constant function depending only on the time.
In [27], Yudovich showed that, in R 2 , as long as the initial vorticity is bounded in L ∞ x , there is global existence and uniqueness of weak solutions to the Euler equations. This encompasses vortex patches, where the initial data for the vorticity is the indicator function 1 A of a set A ⊂ D. The existence result was later generalised by Di Perna and Majda [13] to generic L p spaces.
When the vorticity is unbounded, Delort proved [10] the existence of global weak solutions for two dimensional Euler equations, provided that the initial vorticity is a positive Radon measure in H −1 . This particular class of initial data is interesting for practical applications, since it encompasses vortex-sheets.
In general, uniqueness has been proven to be false for Euler equations (see [24,25], or the recent [21]). Numerically, there is also evidence for non-uniqueness of solutions of the vortex-sheet problem [22].
Our work will focus on the more general framework of (admissible) measure valued solutions (MVS). In this framework (introduced in [14]) the solution is not an integrable function, but rather a parametrised measure ν (x,t) for (x, t) ∈ D × I. Measure valued solutions are shown to exist globally, in 3D, and to represent limits of interesting sequences [13] and numerical approximation schemes [20].
1.2. Numerical methods. Many numerical approximation techniques have been applied successfully to the approximation of the weak formulation of the Euler equations, even in cases where the underlying mathematical analysis in unavailable. However, numerical approximations of problems with rough initial data, may not converge (cf. e.g. [6]). In such cases, it is necessary to question, whether the numerical approximation of weak solutions is an appropriate approach [15].
Existing numerical methods include spectral methods, which, for sufficiently regular solutions, have been shown to converge to the underlying solution with spectral accuracy (exponential accuracy in the number of modes) in a spatially semi-discrete framework [5]. The convergence of spectral methods to measure valued solutions was proved in [20]. However, spectral methods are only valid for periodic domains.
Finite element methods (cf. [18,17]) are also a possible discretization technique, although they generally require some amount of regularity of the initial data.
Finite difference methods can be employed to approximate the equations. This, combined with operator splitting methods, yields efficient and robust methods, cf. [19,9,6,1]. However, no convergence result for non-smooth initial data, including in the context of measure valued solutions, is available for these type of schemes.
1.3. Aim of the paper. Our aim is to propose and analyse a variant of the finitedifference projection methods, and prove that it generates sequences that converge to an admissible measure valued solutions of the incompressible Euler equations.
The procedure used in this paper is based on a similar approach used for the computation of measure valued solutions in the context of hyperbolic conservation laws [15,16]. In the articles, the authors devised a class of schemes that converge to measure value solutions.
The ideas contained in [15,16] can be adapted to the context of incompressible Euler equations. In the context of incompressible flows, the divergence-free constraint imposes restrictions on the type of scheme and time-stepping that can be used for the approximation. For this reason, we will utilise fully discrete conservative fluxes for the nonlinear term.
In [20], the authors present a semi-discrete spectral method that is able to approximate measure valued solutions for incompressible Euler equation in the 2dimensional torus.
One of the novelties of our approach compared with that of [20] resides in the use of finite volume/finite differences as choice of deterministic solver. This, allows the treatment of more complex model problems, such as problems with non-periodic boundary conditions. Furthermore, finite difference methods can be readily used to compute solutions with signed vorticity. For example, in this paper, we consider problems with outflow boundary conditions and positive vorticity. Another advantage of finite differences/volumes is that they are suitable for simulations with very rough initial data (L 2 velocity). Finite volume methods also have the advantage of allowing an easier treatment of complex geometries. In addition, compared to [15,16,20], we perform an analysis with a fully discrete scheme, providing bounds that are independent of the simultaneous refinement of the space and time discretization.
1.4. Contents of the paper. In section 2, we will recall basic definitions of measure valued solutions for incompressible Euler equations. In section 3, we will define the numerical scheme for the approximation of measure valued solutions. Finally, in section 4, we will provide numerical experiments of convergence of the scheme to measure valued solutions.

2.
Admissible measure valued solutions. The concept of measure valued solution is based on that of (generalized) Young measure [26], namely, a probability measure parametrized in space and time.
This notion was initially introduced for conservation laws [11], subsequently extended to uniformly L ∞ -bounded approximations to (incompressible) Euler equations [14] and then generalized to arbitrary L 2 sequences [12], after a small, technical modification taking into account the possibility of concentrations.
We define the family F of functions: We define g ∞ (u) := lim r→∞ g(ru) whenever the limit exists and is unique.
[(Generalized) Young measure theorem [14]] Given a sequence {u n } n∈N of uniformly bounded functions in L 2 (D × [0, T )), there exists: such that, for a subsequence, relabelled as the original by abuse of notation, it holds: using the same notation as (2).
Based on [15], we consider a general measure as initial data for the incompressible Euler equations. Having a measure as initial datum for the problem (rather than an integrable function), further generalizes the problem allowing to model uncertainty in the initial data. This automatically provides a framework for uncertainty quantification.
Definition 2.3 (Admissible measure valued solution). Let σ 0 := ν (x,0) be an initial measure. We deonte by ξ ⊗ξ the mapping of a vector x ∈ R d into the tensor product x ⊗ x ∈ R d 2 . A measure value solution for the Euler equations is a (generalized) Young measure ν such that then we call ν an admissible measure valued solution. Here, we denote λ := λ t × dt and define E(t) := D ν (x,t) , ξ 2 dx + λ t (D) to be the mean (Kinetic) energy of the measure valued solution.
Admissible measure valued solution are not unique [21], thus, a stronger notion of admissibility is necessary to recover uniqueness. However, at least in the case where a smooth solution exists, the notion of admissible measure valued solution coincides with it [7]. The vanishing viscosity approach can be exploited to construct approximations of measure valued solutions [8] from sequences of weak solutions to the Navier-Stokes problem. Another possible approach for the construction of measure valued solutions consists in using a numerical scheme that produces a suitable sequence of approximations. The sequence can be then shown to converge, in the weak* sense, to admissible measure valued solutions.
3. Numerical scheme. In this section, we present a numerical scheme that is able to approximate (admissible) measure valued solutions. In order to compute measure valued solutions, we require two main ingredients: • an "individual solver", able to approximate the Euler equations in their classical formulation, whenever such formulation is well-posed; • a statistical integrator, capable of approximating measures by exploiting stochastic sampling procedures. In general, we cannot expect from the individual solver to converge (weakly or strongly) in the usual sense (cf. [6]). Instead, we require some stability properties allowing us to control key quantities of the approximated solutions, which, in turn, allow us to show consistency with measure valued solutions. Our individual solver is based upon the ideas of [9] and [6]. We utilise a projection method (also known as fractional step or pressure correction method) for an efficient enforcement of the divergence-free condition. The spatial discretization is enforced using a "mimetic" approach: we construct a scheme that possesses properties analogous to the continuous equation.
We consider the domain [0, 1] d × [0, T ]. For simplicity, we assume a uniform grid in both space and time, with uniform grid spacing ∆x = ∆y = 1/M > 0 in space and ∆t > 0 in time. Furthermore, we assume O(∆x) = O(∆y) = O(∆t). The velocity will be approximated by a piecewise constant function in both space and time variables (in a cell C i,j := [i∆x, (i + 1)∆x] × [j∆y, (j + 1)∆y]). We will identify these velocities with a sequence {u n i,j } i,j,n (denoted also u by abuse of notation) indexed by space variables i ∈ I and j ∈ J (I = {0, . . . , M −1}, J = {0, . . . , M −1}) and a time variable n ∈ {0, . . . , N }, T = N ∆t. We denote the space of grid functions for the choice h = ∆x = ∆y as G h .
The scheme uses a fractional step method, in which the velocity is evolved, in a prediction step (4a), without imposing the divergence-free constraint. The divergence constraint is enforced in a second projection step (4b), which also provides an evolution equation for the pressure. The scheme reads: where C(u n ,ū n+1/2 ) i,j is an approximation of the nonlinear convective term (which will be specified later). The operator D is an numerical dissipation term, which ensures the numerical stability of the scheme. The time-averaged velocityū n+1/2 := θu n + (1 − θ)u * ,n+1 (with θ ∈ [0, 1/2) as free parameter) ensures that the scheme is sufficiently dissipative. The approximation of spatial differential operators exploits finite difference (FD) and finite volume approximations (FV).

3.1.
Definitions. In this section, we introduce concrete discretisations of the spatial finite difference operators. We formulate the space discretization in 2 space dimensions, but remark that extension to 3D is straightforward.
3.1.1. Discrete differential operators. We denote by grad h a standard cell-centred central, second order discrete gradient operator defined as ∆y .
We introduce the notation ∆ h (resp. ∆ h ) for a standard five point Laplace (resp. vector Laplace) operator: the discrete L p norm in space and with the L p,q norm in space-time (with standard extensions to ∞). We also introduce the L 2 inner product (in the spatial variables) for grid functions: We introduce the notion of jump a x i+1/2,j := a i+1,j − a i,j and average a x i+1/2,j = 1 2 (a i+1,j + a i,j ) (with an analogous definition for the y-direction).
3.1.3. The (numerical) diffusion operator. The operator D is a second-order, energy dissipative diffusion operator evaluated at the time averaged quantityū n+1/2 := 1 2 (u n+1 + u n ) (this ensures energy dissipation in both space and time). We will use the 2nd order Lax-Wendroff diffusion, which has the form: with 0 < c i±·,j±· < ∞ as parameters of choice.
3.1.4. The exact projection operator. We denote by P the exact discrete projection operator of a vector field into a discretely divergence-free vector space (i.e. s.t. div h P(u) = 0) and Q := I − P (the gradient, or irrotational, part of the vector field, where I is the identity map). The discrete projection is implemented by solving: is a cell-centred five point Laplacian acting on cell-centres in checker-board formation: 3.1.5. The non-linear (convective) term. The non-linear term is evaluated at two different discrete quantities, namely, the "lagged" velocity u n and the time-averaged u n+1/2 . For those we use a second order, consistent flux term: , , We define the cinvenience notation

3.2.
Preliminaries. In this subsection, we discuss some technical propositions which will be used in later sections. First of all, we remark that the fluxes are consistent with the continuous convective terms. Moreover, the fluxes are energy conservative, in the following sense.
Lemma 3.1. The flux function is energy conservative, i.e. satisfies: provided u is discretely divergence free (w.r.t. the operator div h ). Moreover, We omit the proof. Notice that the flux function F x satisfies the following inequality: and similarly for the other arguments of F x and for the same arguments for F y .
3.2.1. Projection. We start by stating two useful lemmas on the exact projection operator.
Lemma 3.2. There exists a grid function ψ such that: Proof. Observe that (4b), together with P(u) = u, the linearity of the projection operator, and the fact Q(u) = grad h ψ for some ψ, implies the claim.
Lemma 3.3. The operators div h and grad h are self adjoint, namely: (For suitable boundaries: the lemma holds for a compactly supported velocity field, wall boundary conditions, periodic boundaries, or a compatible inflow/outflow boundary condition.) The following proposition provides a discrete energy bound (i.e. an L 2 bound for the velocity) and a bound on the gradient pressure for the scheme (4).
Proposition 3.1. The projection P is energy dissipative, namely: Proof. Take the L 2 discrete inner product of (10) with u and use the fact that u is divergence-free: The usual reformulation implies

Energy conservation.
Proposition 3.2. The approximate solution for the projection scheme satisfies the energy balance: u n+1 2 2 = u n 2 2 + 2∆t ū n+1/2 , Dū n+1/2 (12) Moreover, we have by the sign property of the diffusion term and positivity of the other terms.
Proposition 3.3. The following limit holds with a second order Lax-Wendroff diffusion operator D: We will denote these limits as space (resp. time) weak total variation bounds.
Proof. The discrete diffusion operator (of Lax-Wendroff type) provides the spatial total variation limits. Hence, using the energy estimate, we can write Summing over each time-step and exploiting the telescoping summation, we obtain Using the Lax-Wendroff diffusion (7), we obtain ∆x∆y.

FILIPPO LEONARDI
With c i,j > 0 we obtain This proves our claim.
3.3. Convergence to admissible measure valued solutions. In this section we prove that the limit of the sequence of piecewise constant functions generated by the scheme is consistent with admissible measure valued solutions in the limit when the space-time mesh is refined.
Lemma 3.4. The vector valued piecewise constant function R, defined as Proof. Consider the term: . Using the consistency of the fluxes, the last term is exactly f (u n i,j ). Then (using Equation 9): − u n i,j | + 2|u n i+1,j − u n i,j ||u n i,j |. Using the discrete weak total variation bounds, all terms converge to zero after summation over the space-time indices and integration over the probability space Ω. In fact, we have (we consider the first term only, the rest is completely analogous): using the fact that ϕ is sufficiently smooth (i.e. of class C ∞ 0 ). Here, p, q, r, s > 0 are the exponents obtained from the application of the discrete Hölder inequalities.
In the following, we will establish that the weak*-limit of the approximate solution is a measure valued solution.
For the approximation of measure valued solutions with initial data σ x , we proceed as following (cf. [15]): from the law u 0 (ω; x, t) of σ x , we draw a sufficient number of i.i.d. samples u 0 (ω i ; x), for ω i ∈ Ω, i ∈ I for some index set I. We evolve each sample with scheme (4), and find a discretization u h (ω i ; x, t), where h > 0 represent a discretization parameter. From the individual samples, we reconstruct the approximate measure valued solution as: From this approximation, we can compute statistical quantities of interest, such as mean, variance and mean mean kinetic energy.
The following theorem holds: Let σ x be an initial Young measure in L ∞ w (D; P(R d )). Then, there exists a probability space (Ω, F, P ) and a random variable u 0 (ω; x) with law σ x . For fixed ω ∈ Ω, let u h (ω; x, t) be the approximate solution with initial data u 0 (ω; x) obtained with scheme (4). Let ν h (x,t) be its law. Then there exists a subsequence (which we denote by abuse of notation by ν h (x,t) ) and a limit Young measure ν (x,t) , such that, for every divergence-free test vector field ϕ ∈ C ∞ 0 (D × [0, T )) d and q ∈ C ∞ 0 (D × [0, T )), the following holds: moreover the admissibility criterion is satisfied. The limit is therefore an admissible measure valued solution.
Proof. Denote by u h (ω; x, t) and grad h p h (ω; x, t) the approximated step vector fields generated by the numerical scheme (depending on the outcome ω and on the step size h). The fundamental theorem of Young measures provides a measure ν (x,t) , which is the narrow limit of a subsequence ν h (x,t) , the law of the approximation u h . The evolution equations for the projection scheme are condensed into (removing dependency on ω for ease of notation): We denote by C i,j the cell [x i−1/2 , x i+1/2 ] × [y j−1/2 , y j+1/2 ] and by I n the time interval [t n , t n+1 ]. Let ϕ := (ϕ 1 , ϕ 2 ) ∈ C ∞ 0 (D × [0, T )) 2 be a solenoidal test vector field and construct a sequence sequence ϕ h with limit ϕ, such that div h + ϕ h = 0. Multiplying the evolution equation by the corresponding approximate test function and by the cell-indicator function 1 Ci,j ×I n (x, t), summing over each time step n = 0, ..., N − 1 and each spatial index, and integrating in space, time, and on the probability space Ω, we obtain 0 = Dū n+1/2 · ϕ h 1 Ci,j ×I n dxdt dP (ω).
Let us consider each term of the sum individually.
• For the first term, we define: Consider: As ϕ is compactly supported in [0, T ) we have ϕ N i,j = 0, ∀i ∈ I, j ∈ J for a large enough N (i.e. for small enough ∆t). Hence (using Fubini's theorem): Thus, using the weak*-convergence and the consistency with the initial data, we obtain • The pressure term is tested against discretely divergence-free vector field, hence the term drops upon integration. • Let us now consider the convection term. We notice that the divergence for the flux form is considered as an edge-centred sequence. Thus, integration by parts with a cell-centred quantity leads to a forward gradient. The convective term becomes: using consistency, the weak*-convergence of the Young measure (by Young's theorem), Fubini's theorem, and the fact that ϕ ∈ C ∞ 0 . • For the dissipation term, consider: Dū n+1/2 · ϕ1 Ci,j ×I n dxdt dP (ω).
The first term can be bounded in the following way (assuming boundedness of the coefficients c i,j of the diffusion): Ω i∈I j∈J by the time total variation limit and smoothness of ϕ. Having space bounds on the averaged state it follows immediately that the term goes to zero. • The divergence-free constraint is treated analogously. As h → 0, the approximate measure generated by the scheme converges to a measure valued solution.
Given the previous theorem, we have a tool for the consistent approximation of measure valued solutions from a measure as initial data. For an atomic initial data (i.e. σ 0 = δ u0 ), we can define an analogous algorithm. The algorithm is the following (cf. also [15,20]) and uses a Monte-Carlo sampling procedure: 1. set a small discretization parameter δ > 0, and create a new random initial data u δ,0 (ω; x), adding a random perturbation of size δ to the initial data u 0 (x); 2. draw M > 0 i.i.d. samples from the random field arising from the perturbed u δ,0 (ω; x); 3. evolve each sample, individually, using (4), after choosing suitable space-time discretization parameters; 4. compute mean, variance and other quantities of interest form the samples. 4. Numerical experiments. The following computations are performed on the Cray XC40 Piz Dora of the Swiss National Supercomputing Center and on the ETH Zürich cluster Euler. The code was developed in C++ and exploits the parallel framework PETSc [3,4] to achieve efficient load balancing on distributed domains. 4.1. Vortex patch. As first set of experiments for computation of measure valued solutions, we consider a vortex patch. We call vortex patch an initial data whose vorticity takes the form of an indicator function of a set. The vorticity is bounded and weak solutions with vortex patches as initial data are unique and well defined in 2D and with periodic boundary conditions [27].
Numerical evidence indicates that, as δ → 0, the approximate measure valued solution converges towards a limit measure (cf. Figure 2). Moreover, the perturbed measure valued solution converges towards an atomic solution (1): i.e. towards the Dirac measure of the indicator function of the ball of radius r 0 (the initial data with δ = 0).

4.2.
Perturbed vortex sheet. We now consider a so-called vortex-sheet. The vortex-sheet is a standard benchmark commonly used in engineering. Contrary to the vortex-patch, the vortex-sheet is ill-defined and present large instabilities for small times. The initial data for the velocity is where H(y) = 1 for y > 0 and H(y) = −1 for y ≤ 0. This initial data is a steady state for the Euler equations.
We consider a small random perturbation of this initial state. To this end, we define a perturbation by setting the initial data for the velocity to be the random field: u δ 0 (ω; x, y) := P(u 0 (y + f δ (ω; 2πx)), 0), where f δ is the same perturbation function used in the case of the vortex patch.
Here P denotes the Leray projection. As boundary conditions, we impose periodic condition in the x-direction and wall boundaries in the y-direction (where u = 0). We choose, as final time, T = 1 and as domain [0, 1]×[−1, 1] to ensure the turbulence zone does not reach the wall boundaries. We perform experiments with a random space of dimension M = 20. We compute the sample mean and the sample variance of the ensemble (cf. Figure 4). The error of the mean and of the variance of the approximate solution converges to zero for any time (cf. Figure 5), even if the individual samples do not. Moreover, looking at the distributions of the components of the velocity at fixed points in time and space, the variance does not converge to zero and, in fact, the limit measure appear to be non-atomic (cf. Figure 6).

5.
Conclusion. In this paper, we presented a fully discrete Monte Carlo finite difference projection method for the approximation of measure valued solutions of the incompressible Euler equations. We showed the consistency and weak* convergence of our method. Those properties are derived from the interesting stability features of our scheme.  By applying similar techniques to [15,16] (for hyperbolic conservation laws), we were able to construct a similar finite difference/finite volume method, that is applicable in the context of incompressible flows. Furthermore we were able to extend the convergence to a fully discrete scheme, and consider the simultaneous refinement of both the space and time discretisations.
Numerical experiments indicate that vortex-patches (well defined in the weak sense) and vortex-sheets have very different underlying behaviours for the measure   valued solutions: the first appear to converge to atomic measure valued solutions (i.e. to standard weak solutions), whilst the second appear to converge to nonatomic weak solutions. Our findings strongly suggest that some initial data, even if atomic, may generate measure valued solutions that are non-atomic. Our results are consistent with the ones illustrated in [20]. By using finite volume/finite difference methods for the deterministic solver, we are able to consider boundary conditions that differ from periodic boundary conditions. All difference operators can be generalised to unstructured grids, and, combined with a finite element solver for the Poisson equation, the scheme can be generalised to arbitrary geometries.