A semi-implicit finite volume scheme for dissipative measure-valued solutions to the barotropic Euler system

A semi-implicit in time, entropy stable finite volume scheme for the compressible barotropic Euler system is designed and analyzed and its weak convergence to a dissipative measure-valued (DMV) solution [E. Feireisl et al., Dissipative measure-valued solutions to the compressible Navier-Stokes system, Calc. Var. Partial Differential Equations, 2016] of the Euler system is shown. The entropy stability is achieved by introducing a shifted velocity in the convective fluxes of the mass and momentum balances, provided some CFL-like condition is satisfied to ensure stability. A consistency analysis is performed in the spirit of the Lax's equivalence theorem under some physically reasonable boundedness assumptions. The concept of K-convergence [E. Feireisl et al., K-convergence as a new tool in numerical analysis, IMA J. Numer. Anal., 2020] is used in order to obtain some strong convergence results, which are then illustrated via rigorous numerical case studies. The convergence of the scheme to a DMV solution, a weak solution and a strong solution of the Euler system using the weak-strong uniqueness principle and relative entropy are presented.


Introduction
Hyperbolic systems of conservation laws arise in the modelling of various physical phenomena where finite speeds of propagation and conservation of quantities are involved.The most notable among them is, invariably, the non-linear Euler system of gas dynamics which represents the basic conservation laws of mass, momentum and energy.A common trait of non-linear conservation laws is that discontinuities in the form of shock waves may develop in finite time, even if the given initial data was sufficiently regular.As a result, the concept of weak or distributional solutions to conservation laws was introduced as an appropriate replacement.But, it is well known that these weak solutions may not be unique.As conservation laws usually model the real world, a selection criterion in the form of a relevant entropy condition is usually enforced in order to rule out the non-physical solutions.The enforced entropy condition is most often a mathematical equivalent of the second law of thermodynamics.This approach has yielded successful results in the scalar case; see [4,5,42].However, in light of recent results by Chiodaroli et al. [8] and De Lellis et al. [10], the uniqueness of entropy weak solutions fails to hold in the multidimensional case for the barotropic Euler system, and the result is extended to the complete Euler system by Feireisl et al. [21].The ill-posedness of the Euler system seems to stem from the lack of compactness of entropy weak solutions, in the sense that even a bounded sequence of solutions may develop oscillations or concentrations or both.
Inspired by the previously stated results, there was a renewed interest in the concept of measure-valued (MV) solutions, introduced by DiPerna [13] in the context of general conservation laws.This was further extended to incompressible fluid flows by DiPerna and Majda [14] and by Neustupa [45] for the compressible barotropic Euler and Navier-Stokes systems.The main advantage of MV solutions is that they can be identified as cluster points or limits of oscillatory approximated sequences; see the monograph [44].The measure in a MV solution refers to the Young measures being used to describe the possible oscillations in the approximate sequence.In the present work, we focus on the so-called dissipative measure-valued (DMV) solutions, first introduced by Feireisl et al. [18] for the compressible Navier-Stokes system and later on for the Euler system in [6,19].The class of DMV solutions can be seen as a natural closure in the weak topology of the solution set to the underlying system, with the dissipative solutions being the barycenters of the associated Young measure.The major interest in this approach is because the phenomenon of "weak-strong uniqueness" seems to hold in the class of DMV solutions, i.e. a DMV solution and a strong solution emanating from the same initial data seemingly coincide provided the latter exists.We refer the interested reader to [6,18,19,24,33,52] for an exhaustive survey of the same.
Due to the existence of discontinuities in weak solutions, numerically approximating solutions of general conservation laws is challenging.Conventional finite difference approximations are known to be inappropriate near discontinuities where the PDE does not hold.The finite volume methods which captures discontinuities with high resolution present an attractive, robust and cost effective alternative to classical finite difference schemes.Over the past few decades a great deal of progress has been made about the finite volume methods and today, a variety of these methods are available.We refer the reader to [17,31,43,51] for an extensive survey of numerical methods for general hyperbolic problems.The convergence of finite volume methods to weak solutions of the compressible Euler system has been studied extensively in [34][35][36][37][38].The consistency analysis, in the spirit of Lax-Wendroff, is carried out across all the above mentioned references under strong assumptions, which may fail to hold in general.However, in the framework of DMV solutions, provided some consistency and stability estimates are established, one can show that the numerical solutions generated by the finite volume method converge, albeit in the weak sense, to a DMV solution of the Euler system under some physically reasonable boundedness assumptions.The current state of the art literature pertaining to the convergence of finite volume methods to DMV solutions and in general, MV solutions, includes [20,[22][23][24][25][26][27][28] and the references therein.
Due to the convergence of the numerical solutions being weak and the limit being a DMV solution, there are difficulties in effectively capturing the limit solution.A celebrated result by Komlós [41] asserts that any sequence of uniformly bounded  1 -functions {  } ∈N admits a subsequence {   } ∈N such that the arithmetical averages 1  ∑︀  =1    converge pointwise a.e. to a function  in  1 , with any further subsequence of {   } ∈N also enjoying the same property.Balder [2] adapted this result to a measure-valued framework and introduced the novel concept of -convergence for sequences of Young measures.Using this, one can show that the arithmetical averages of the numerical solutions converge strongly to the expectations (or barycenters) of the Young measure associated to the DMV limit; see [23][24][25] for more details.

Aims and scope of the paper
The aim of the present paper is to design and analyze a semi-implicit in time, entropy stable finite volume method for the compressible barotropic Euler system and show its weak convergence to a DMV solution.The main motivation in choosing a semi-implicit setup is that there is a considerable reduction in the computational complexity when implementing the scheme, as it avoids the need to invert large and dense matrices which is a common occurrence in fully-implicit schemes.Further, the restrictive stability conditions that are typical of fully-explicit schemes are also relaxed.The entropy stability is achieved by the introduction of a shifted velocity in the convective fluxes of the mass and momentum balances, with the shifted velocity being proportional to the pressure gradient, provided some CFL-like condition is satisfied.A consistency analysis is performed in the spirit of the Lax's equivalence theorem under some boundedness assumptions, which reflect the numerical solutions staying in a non-degenerate region.The limiting process is then carried out in which we show that the sequence of numerical solutions contains a subsequence that generates a Young measure, and the subsequence further converges to the expectations (or barycenters) of said Young measure in the weak-* sense.We then employ the techniques of -convergence to obtain the convergence of the Cesàro averages of the numerical solutions to the expectations of the Young measure in any   .We also deduce the convergence of the -Wasserstein distance,   , between the Cesàro averages of the Dirac masses centered at the numerical solutions and the Young measure in any   , along with the convergence of the absolute first variance or the  1 -deviations of the numerical solutions, which reflects the probabilistic nature of the limit solution.Said convergences are then verified via rigorous numerical case studies, where we exhibit the convergence of the scheme to a DMV solution, a strong solution and a weak solution of the Euler system using the tools of the weak-strong uniqueness principle and relative entropy.The convergence to a weak solution is illustrated using a very specific setup wherein we consider a Riemann problem in one space dimension with a parameterized pressure law.Then, we simulate the pressure going to 0 by taking decreasing values of the above parameter.The vanishing pressure causes the formation of a -shockwave, which is mathematically represented by a Dirac measure, leading to the solution being interpreted as a distributional solution of the system.

Barotropic Euler system
We consider the following initial value problem for the barotropic Euler system in (0,  ) × Ω which reads where the variables  = (, ) and  = (, ) are the density and velocity of the fluid respectively.Here,  = () denotes the pressure and for the sake of simplicity we consider the isentropic pressure law with  > 0 and  > 1 denoting the ratio of specific heats also known as the adiabatic constant.We also set as the so-called pressure potential or the internal energy per unit volume.Note that   along with  satisfies  ′  () −   () = (). (1.4) In the following, we consider the domain Ω ⊂ R  ,  = 1, 2 or 3, to be bounded together with the space-periodic boundary conditions.We recall from [38] the following energy estimates satisfied by regular solutions of (1.1a) and (1.1b) Proposition 1.1 (A priori energy estimates).Regular solutions of (1.1a) and (1.1b) satisfy: (1) A renormalization identity: (2) A kinetic energy identity: (3) The total energy identity/entropy identity: where  = 1 2 || 2 +   ().

Dissipative measure-valued solutions
Dissipative measure-valued (DMV) solutions can be seen as a generalization to the concept of weak solutions of the barotropic Euler system.We follow the formulation proposed in [24] in order to recall the notion of a DMV solution.The interested reader can refer to alternative formulations presented in [6,19,22,23] and the references therein.
The phase space for the barotropic Euler system is defined as and we let (ℱ) denote the space of all proability measures on ℱ.
∈ Ω.Further, the following integral identity

Weak-strong uniqueness principle
Though DMV solutions are much weaker objects when compared to weak solutions, the strong solutions of the Euler system are stable in the class of DMV solutions, i.e. any DMV solution and strong solution of the Euler system that emanate from the same initial data necessarily coincide, provided the latter exists.This fundamental property is labelled as the "weak-strong uniqueness principle" in literature; see [6,18,19,24,33,52] for more details.The key tool in proving this principle is the so-called relative entropy functional, which is defined as where ,  are test functions mimicking a strong solution of the barotropic Euler system.
Remark 1.4.Though (1.18) is not a metric as it is not symmetric, the following still hold: We now recall the statement of the theorem following [24].
Theorem 1.5 (Weak-strong uniqueness principle).Let a parameterized family of probability measures  = { , } (,)∈(0, )×Ω be a DMV solution of the barotropic Euler system in the sense of Definition 1.2, with initial data [ 0 ,  0 ] and associated concentration defect measures C  and R  .
Let ,  be a strong solution of the barotropic Euler system belonging to the class

Velocity stabilization
As mentioned earlier, the motivation for this work is to introduce a finite volume scheme to approximate the barotropic Euler system and to study its stability in the sense of energy stability, i.e. to obtain a discrete equivalent of (1.7).In order to enforce the stability at a numerical level, we adopt the formalism introduced in [9,15,16,32,48] wherein a stabilization term in the form of a shifted velocity is introduced in the mass and momentum balances to yield the following modified system: Analogous to Proposition 1.1, we can show that the solutions to the modified system satisfy the following a priori energy estimates.
Proposition 1.6 (A priori energy estimates for the modified system).Regular solutions of (1.20a) and (1.20b) satisfy (1) a renormalization identity: (1.21) (2) a kinetic energy identity: (3) the entropy balance: Therefore, at least at the continuous level, we can see that formally setting  = ∇  ,  > 0, gives us We thus observe that introducing a shift in the velocity indeed has a stabilizing effect, allowing us to obtain the desired entropy inequality.

The finite volume method
We wish to numerically approximate the velocity stabilized Euler system (1.20a), (1.20b) and to this end, we start with a space domain Ω such that the closure of Ω is a union of closed rectangles ( = 2) or closed parallelepipeds ( = 3).

Mesh and unknowns
We introduce a tessellation  of Ω ⊂ R  , known as the primal mesh, consisting of possibly non-uniform closed rectangles ( = 2) or closed parallelepipeds ( = 3) such that Ω = ∪ ∈ , where  is called a control volume.The -dimensional Lebesgue measure of  will be denoted by ||.The set of all edges ( = 2) or faces ( = 3) (hereafter commonly referred to as edges) of all control volumes  ∈  is denoted by ℰ and by ||, we denote the ( − 1)-dimensional Lebesgue measure of  ∈ ℰ.By ℰ int , ℰ ext and ℰ(), we denote the subsets of all internal edges, external edges, i.e. edges lying on Ω, and the edges of a control volume  ∈  , respectively.Any two cells ,  ∈  share a common edge denoted by  = |.For  ∈ ℰ(),  , denotes the unit normal to  pointing outwards from .By   we mean  ≤  for a constant  > 0, independent of the mesh parameters.
The mesh size ℎ  is defined by where ℎ  = diam() and we assume that there exists a constant  > 0 (independent of other mesh parameters) such that ℎ  ≤ inf{ℎ  :  ∈  }.
By ℒ  (Ω) ⊂  ∞ (Ω), we denote the space of all scalar-valued functions constant on each cell  ∈  with the associated projection map Π  :  1 (Ω) → ℒ  (Ω) defined as where   denotes the indicator function of .Analogously, we define ℒ  (Ω; R  ), the space of all vector-valued piecewise constant functions with the projection operator being defined componentwise.
For  ∈ ℒ  (Ω),   denotes the constant value of  on the cell  ∈  .Further, if we have  = | ∈ ℰ(), we define the average value of  across  as The discretization is collocated in the sense that the discrete unknowns are all associated to the same location.The discrete unknowns are thus (  ) ∈ for the density and (  ) ∈ for the velocity and we denote the pressure

. Discrete mass flux and discrete differential operators
We introduce the discrete mass flux and the discrete differential operators on the function space ℒ  (Ω) defined above.
where we define The positive and negative parts of  , will be specified later once we give a description of the numerical scheme.
Definition 1.8 (Discrete gradient and divergence).In the case of a collocated grid, we follow [37] in defining the discrete gradient and divergence operators so that the classical grad-div duality still holds at the discrete level.
The discrete divergence operator div  : while the discrete gradient operator (1.30) We now report the following lemma from [37].

The scheme
Let us consider a discretization 0 =  0 <  1 < • • • <   =  of the time interval [0,  ] and let  =  +1 −   for  = 0, 1, . . .,  − 1 be the constant timestep.We initialize the scheme by taking the initial approximations to be the average values of our initial data  0 and  0 respectively, i.e.
We now consider the following fully discrete scheme for 0 ≤  ≤  − 1: where   =   −  +1 denotes the stabilized velocity and we set  +1  = (∇   +1 )  with  > 0, a parameter to be determined later; cf.(1.24).The positive and negative parts of the term  , appearing in the mass flux  , are defined so as to bring about an upwind-bias and maintain the signs split in  , .
The discrete momentum update (1.35b) along with the discrete mass update (1.35a) yields the following update of the velocity component which is used in establishing the entropy stability. where

Main results and organization
The proposed scheme (1.35a) and (1.35b) involves solving a non-linear system at each time step, mainly due to the presence of non-linear stabilization terms in the mass balance.However, once we compute the updated density  +1 from (1.35a), the updated velocity  +1 can be computed in an explicit manner from the momentum balance (1.35b).The existence of a discrete solution to the scheme is established following a standard topological degree argument, analogous to the treatment in [1,30,46,47], which uses classical tools from topological degree theory in finite dimensions; see [11].For the sake of completeness, we just state the existence result.
Once the existence of a discrete solution is known, we can define the notion of a numerical solution of the scheme as follows.
Definition 1.12 (Numerical solution of the scheme).Let {(  ,   ) ∈ ℒ  (Ω) × ℒ  (Ω; R  ) :  = 0, . . .,  − 1 ;   > 0} be a family of discrete solutions to the scheme (1.35a) and (1.35b).Then a numerical solution of the scheme corresponding to the space-time discretization ( , ) is defined to be the pair of piecewise constant functions Though Theorem 1.11 guarantees the positivity of density, we cannot claim or prove the existence of suitable a priori bounds that guarantee the numerical density being bounded below, away from zero or being bounded above.The density being bounded uniformly from above and below is crucial in the analysis carried out later.Hence, we impose the following as our principal working hypothesis; see also [22,24].Hypothesis 1.13.There exist constants  ≥  > 0 (independent of the mesh parameters) such that 0 <  ≤   , ≤ . (1.39) We can now state the main result concerning the entropy stability of the proposed scheme and we defer the proof to Section 2. Theorem 1.14 (Local in-time entropy inequality).Any numerical solution to the scheme (1.35a) and (1.35b) satisfies the following entropy inequality where )︁ .
In addition, we also get the following global entropy estimate.
Theorem 1.15 (Global entropy estimate).Assume that the hypothesis of Theorem 1.14 holds.Then, there exists a constant  > 0 such that the discrete solutions (  ,   ) 0≤≤ satisfy the following global entropy inequality: It is convenient to perform the consistency analysis of the scheme on the dual mesh.However, we remark that the dual mesh is not needed to implement the scheme.In order to perform the said analysis, in Section 3, we introduce the dual mesh  corresponding to the primal mesh  .We also introduce a reconstruction operator ℛ  which transforms a function which is piecewise constant on the control volumes of the primal mesh into a function which is piecewise constant on the dual cells; see [37] for more details.We now present a result related to the weak convergence of these dual mesh reconstructions and postpone the proof to Section 3. Lemma 1.16 (Weak convergence of dual mesh reconstructions).Let ( () ) ∈N be a sequence of space discretizations such that the mesh size ℎ  () → 0 as  → ∞.Let 1 ≤  < ∞ and let  ∈   (Ω).Let { () } ∈N be a uniformly bounded sequence in   (Ω) such that  () ∈ ℒ  () (Ω) for each  ∈ N and  () ⇀  in   (Ω) as  → ∞.Let ℛ  () denote the reconstruction operator on the mesh  () .Then, ℛ  ()  () ⇀  in   (Ω) as  → ∞.
Before we state the consistency result, we would like to remark that as a consequence of the timestep restriction imposed by Theorem 1.14,  → 0 as ℎ  → 0. We now present the weak consistency result of our scheme and give a proof of the same in Section 3.
Theorem 1.17 (Consistency of the numerical scheme).Let ( , ) be a given space-time discretization of our domain [0,  ] × Ω.Let (  , ,   , ) denote a numerical solution to the scheme (1.35a) and (1.35b) in the sense of Definition 1.12 with the initial data given by (1.34).In addition to Hypothesis 1.13, we make the following assumption: -There exists a constant  > 0 (independent of the mesh parameters) such that |  , | ≤ . (1.42) Then, the scheme (1.35a) and (1.35b) is consistent with the weak form of the Euler equations (1.1a) and (1.1b), i.e.
for any  ∈  ∞  ([0,  ) × Ω; R  ).Here,  , and  , denote the reconstructions of   , and   , on the dual mesh and the consistency errors ℛ mass  , and ℛ mom  , are such that We finally present the key result related to the weak convergence of the scheme and prove it in Section 4.
Theorem 1.18 (Weak convergence).Let ( () ,  () ) ∈N be a sequence of space-time discretizations such that ℎ () = ℎ  () → 0 as  → ∞.Let ( () ,  () ) ∈N ,  () =   () , () and  () =   () , () , be the sequence of numerical solutions generated by the scheme.Suppose that the hypotheses of Theorems 1.14 and 1.17 are satisfied.Let  () =  ()  () .Then, there exists a subsequence (not relabelled) such that where  = { , } (,)∈(0, )×Ω is a DMV solution of the barotropic Euler system (1.1a) and (1.1b) (in the sense of Def.1.2) with C  = 0 and R  = 0. (1.46) Once the weak convergence of the numerical solutions is achieved, we apply the techniques of -convergence in order to better visualize the limit solution.The key idea is that averaging the sequence of numerical solutions "compactifies" the sequence, which in turn allows us to obtain a strong limit (possibly for a subsequence).The following result can be obtained, and the proof follows analogous lines as in [24,25].
Theorem 1.19 (-convergence).Let the assumptions of Theorem 1.18 hold.Then, we have the following convergences (passing to another subsequence if needed).
(i) Strong convergence of the Cesàro averages (iii)  1 -convergence of the deviations or the first variance Finally, in Section 5, we present the results of the numerical case studies performed, and exhibit the convergences claimed in Theorem 1.19.

Entropy stability of the scheme
This entire section is devoted to the proof of Theorem 1.14, and we begin by proving the discrete analogue of Proposition 1. 6.In what follows, by ,  we denote the interval [min{, }, max{, }].
Theorem 2.1.Any numerical solution to the scheme (in the sense of Def.1.12) satisfies the following discrete identities for each 0 ≤  ≤  − 1: -a discrete renormalization identity: where the remainder term  +1 , is given by where the remainder term  +1 , is given by -a discrete entropy identity: (2.5) Proof.See Appendix A.1 for the proof of (2.1) and Appendix A.2 for the proof of (2.3).Identity (2.5) then follows as the sum of (2.1) and (2.3).
We can now give the proof of Theorem 1.14.
Proof of Theorem 1.14.Multiplying the discrete entropy identity (2.5) by || and then summing over all  ∈  , we obtain wherein we have used the conservativity of the flux and the discrete duality (1.32).Note that ∑︀ ∈ || +1 , is unconditionally non-negative.The convexity of the map  ↦ → || 2 allows us to apply the Jensen's inequality in the velocity update (1.37) to yield (2.7) We can estimate the first term on the right hand side using the Cauchy-Schwarz inequality to get Combining (2.4) along with the estimates (2.7) and (2.8) and using them in (2.6), we finally obtain (2.9) The first term on the right hand side is non-positive under the timestep restriction The proof of Theorem 1.15 follows due to (2.9).

Consistency formulation of the scheme
The goal of this section is two-fold: one to establish the weak convergence of dual mesh reconstructions (Lem.1.16) and the other to rigorously show that the proposed scheme is indeed consistent with the weak formulation of the barotropic Euler system (Thm.1.17).In particular, we show that the numerical solutions satisfy the weak formulation of the continuous Euler system modulo some perturbation terms that vanish in the limit ℎ  → 0. As stated earlier, it is convenient to perform the consistency analysis on the dual mesh and for that purpose, we recall the notion of a dual mesh along with a discrete gradient defined therein; see [29,37] for further details.

Useful estimates
A few basic inequalities used in the numerical analysis that follows are listed below and we refer the interested reader to [17,20,24,30] for more details.Assuming that  ∈  ∞  (Ω),  ∈  ∞  (Ω; R  ) and 1 ≤  ≤ ∞, the following hold: (3.4)

Dual mesh reconstructions
Following [37], we introduce the notion of a reconstruction operator in order to mainly simplify the proof of the consistency result and recall a stability result of the same.Definition 3.2 (Reconstruction operator).Given a primal mesh  of Ω and a corresponding dual mesh , a reconstruction operator is a map ℛ  : ℒ  (Ω) → ℒ  (Ω) defined as where q = We now establish the following estimate, which is analogous to the result ( [30], Lem.3.5).
We can now give the proof of Lemma 1.16.
Then, noting that ℛ  Π   () = ℛ   () as  () ∈ ℒ  () (Ω) for each  ∈ N, we get Every term on the right hand side can be estimated using a combination of Lemmas 3.3, 3.4 and the estimates provided by (3.4) as follows: Combining all the above estimates, we get < , ∀  ≥ .
We now have all the necessary tools required to prove Theorem 1.17.
Proof of Theorem 1.17.First, we prove the consistency of the mass balance (1.35a).Let  ∈  ∞  ([0,  ) × Ω) and let   , denote its interpolate on the primal grid given by (3.1).Multiply the mass balance (1.35a) with ||   and sum over all  ∈  and 0 ≤  ≤  − 1 to obtain where )︁    . (3.9) Analogous to Theorem 3.5 of [36] and Theorem 4.5 of [38], we can write  1 as Combining the terms, we obtain where The terms  1 ,  3 are readily estimated using the strong convergence of the interpolates   , , ð  , to ,    respectively along with the presumed assumptions.To estimate  2 , we use the estimate provided by (3.4) to obtain Finally,  2,2 can be simplified in the same manner as  2,1 and further, using Lemma 3.3 and the pressure gradient estimate (2.11) along with Hypothesis 1.13, we get Thus, we obtain that Next, we show the consistency of the momentum balance (1.35b).Let  ∈  ∞  ([0,  ) × Ω; R  ) and let   , denote its interpolate on the primal grid given by (3.1).Take the dot product of the momentum balance (1.35b) with ||   and sum over all  ∈  and 0 ≤  ≤  − 1 to obtain where Again, proceeding as in Theorem 3.5 of [36] and Theorem 4.5 of [38], we can write  1 as Now, we use the discrete grad-div duality (1.32) to rewrite  3 as Analogous to the mass balance, we split  2 as  2 =  2,1 +  2,2 with the simplified form of  2,1 given by and  2,2 reads Combining the terms, we obtain where 1 ,  4 are readily estimated due to the strong convergence of the interpolants to the respective classical variants along with the presumed assumptions. 2 and  3 are handled using the estimates provided by (3.4) as follows: Finally,  2,2 can be estimated using Lemma 3.3 and the pressure gradient estimate (2.11) along with Hypothesis 1.13 to yield Thus, we obtain |ℛ mom  , | = (ℎ  , √ ) and this completes the proof.
Passing to the limit  → ∞ along with the use of Lemma 0 +   ( 0 ), proving that  is indeed a DMV solution of the Euler system which concludes the proof.

Numerical results
The goal of this section is to present the results of the numerical case studies.The implementation of the scheme (1.35a) and (1.35b) is done as follows.The mass balance (1.35a) is solved first in order to obtain the updated density  +1 , which is done using an iterative Newton method due to the presence of nonlinear terms.Then, the momentum equation (1.35b) is evaluated explicitly to get the updated velocity  +1 .
In accordance with Theorem 1.14, note that the timestep  needs to satisfy and  ≥ ( +1  ) −1 .Both these conditions are implicit in nature and since the former involves the flux as well, they are difficult to implement in practice.However, we remark that if the following implicit restriction holds Therefore, we get that 3 2   ≥ 1  +1  and as a consequence, choosing satisfies the necessary condition  ≥ ( +1  ) −1 .We now state a sufficient timestep restriction which is easy to implement in practice for the scheme (1.35a) and (1.35b).
Proof.The proof follows exactly as given in Proposition 3.2 from [9] and hence we omit the details.
Remark 5.2.Although the sufficient timestep condition (5.3) is still implicit in nature, we implement it explicitly during the computations as also done in [1,15].
In all the numerical experiments that follow, we consider a uniform Cartesian mesh.Let   denote the solution computed on a mesh of  points for a 1 problem and a grid of  ×  points for a 2 problem.Let   and Ũ respectively denote the Cesàro average of the numerical solutions and their first variance, i.e.
Let   be the reference solution computed on the finest possible mesh.In the experiments that follow, the reference solution is always computed using an explicit Rusanov scheme.Analogous to [24,25], we compute the four errors where ‖•‖ denotes the  1 -norm and   , = 1
Remark 5.3. 1 above denotes the 1-Wasserstein distance and it is defined as follows: Given probability measures  and  on a subset  ⊂ R  , the 1-Wasserstein distance between them is defined as in [24].
where ( × ) denotes the set of probability measures on  ×  and  1 ,  2  represent the projections (marginals) of the measure  ∈ ( × ), i.e.  1 () = ( × ) and  2 () = ( × ) for any subset  of .In the experiments that follow, this distance is computed using the wasserstein distance function available in the scipy.statspackage of the open source library SciPy, cf.[40].
In Figure 1, we present the pseudocolor plots of the density  and the momentum  = ( 1 ,  2 ), with  = 1024, which clearly shows the resolution of the circularly expanding shock wave.In Figure 2, we present the plots of the errors  1 , . . .,  6 for the density and the momentum.We observe that classical convergence of individual numerical solutions is not obtained, in particular due to the increase in  1 error for the momentum components.However, the approach using -convergence yields the convergence of the Cesàro averages along with the convergence of the  1 -deviations and the Wasserstein distance, leading us to conclude that the numerical solutions indeed converge to a DMV solution of the Euler system.Also, the claims made regarding the convergence of the Cesàro averages and the Wasserstein distance in any   -norm seem to hold, described by the errors  5 and  6 respectively.The values of the errors  1 , . . .,  6 of  along with their convergence rates are presented in Table 1.As an additional observation, one can see that the errors of the momentum components end up coinciding due to the symmetric nature of the problem.The maximum value of  observed for this problem was 1.5.
Remark 5.4.For the above problem, the  1 error in the case of density is decreasing, indicating the strong convergence of the sequence of numerical solutions { () } ∈N in  1 .In accordance with Proposition 6.13 of [49], the measures  , generated by the sequence of numerical solutions of the above problem can be decomposed as  , =  (,) ⊗  , for a.e.(, ) ∈ (0,  ) × Ω, where  denotes the  1 -limit of the sequence { () } ∈N and { , } (,)∈(0, )×Ω is the Young measure generated by the sequence { () } ∈N .

Kelvin-Helmholtz instability
We consider the classical problem of Kelvin-Helmholtz instability.The set up is same as that of [50].The computational domain is set as [−0.5, 0.5] × [−0.5, 0.5] with periodic boundaries everywhere.The initial density reads and the pressure law is given by () =  5/3 .A cartesian shear flow is set up in the  1 -direction with velocity    Such a configuration is known to be unstable at all wavelengths.We seed the instability by applying a small velocity perturbation in the  2 -direction given by where we set  = 0.025 and  = 1/6.For the chosen density ratio of 2 : 1, the time after which the instabilities will start to develop is around  = 0.35; see [50] for details.Hence, we set the final time as  = 0.4 and also set  = 1024.
In Figure 3, we present the density profiles for successive mesh refinements and we can see that the instabilities are just starting to develop when  = 512.In Figure 4, we present the error plots for the density  and the momentum  = ( 1 ,  2 ).The numerical results clearly show that all the errors are decreasing upon refining the mesh, leading us to conclude that the Young measure generated by the numerical solution is the family of Dirac masses centred at the limit solution.We also compute the  1 -norms of the relative entropy of the numerical solutions with respect to the reference solution, for each , in order to determine whether the limit   solution is a strong solution or a weak solution of the Euler system.The results are presented in Table 2 and we observe that the norm of the relative entropies are also decreasing, indicating that the limit is a strong solution of the Euler system in accordance with the weak-strong uniqueness principle; cf.Remark 1.4 and Theorem 1.5.The maximum value of  ∼ 1.76 for this problem.
The pressure law is taken as () =  2  1.4 ,  > 0. It has been shown in [7] that as  → 0, the isentropic Euler system with the above pressure law reduces to the zero pressure equations of sticky particles.The zero pressure model admits vacuum states and -shockwaves due to coincident and resonant characteristic fields.Consequently, any two-shock Riemann solution tends to a -shock of the zero pressure Euler equations and the density between the two shocks tends to a weighted Dirac-measure which in turn forms the -shock; see [7] for more details.We compute the numerical solutions for  = 1, 10 −2 , 10 −3 , 10 −5 , with final time  = 0.2 and  = 2048.For the sake of brevity, we present the plots of the density when  = 1024 , for the cases of  = 1 and  = 10 −5 in Figure 5. Clearly, we see the formation of classical shock wave corresponding to  = 1, whereas the figure for  = 10 −5 indicates the formation of vacuum states and a -shock.We did not observe any instability or the breaking down since the present scheme is positivity preserving.Furthermore, the proposed scheme is able to effectively capture non-classical -shock type discontinuities.
In Figures 6, 7, 8 and 9, we present the error plots of the density  and the momentum , for  = 1, 10 −2 , 10 −3 and 10 −5 respectively.For each , all the errors decrease upon refining the mesh, indicating that the limit of the numerical solutions for each case is a weak solution of the Euler system.Therefore, we also compute the  1norms of the relative entropy in each case to check whether convergence towards a strong solution of the system is achieved or not.The results are presented in Figure 10 and we observe that for  = 1, the limit is indeed a strong solution of the Euler system, cf.Remark 1.4 and Theorem 1.5.For the cases of  = 10 −2 , 10 −3 and 10 −5 , it seems that we have no convergence towards a strong solution due to the disordered nature of the plots, which is in line with the theoretical expectations.The maximum value of  observed was  ∼ 7.5, 93.05, 98.18, 98.29 for  = 1, 10 −2 , 10 −3 , 10 −5 respectively.Remark 5.5.We note that in the above problem, corresponding to  = 10 −5 , we observe the violation of theoretical estimates assumed in Theorem 1.17, i.e. related to the boundedness of the numerical density.This indicates that the present scheme is robust and that it can perform in the degenerate cases where vacuum is formed and the density is unbounded above.

. 10 )
while the second term is non-positive under the condition  ≥ 1

Figure 1 .
Figure 1.Density and momentum profiles of the cylindrical explosion problem with  = 1024 at time  = 0.25.

Figure 2 .
Figure 2. Error profiles for the cylindrical explosion problem.

Figure 3 .
Figure 3. Density profiles of the Kelvin-Helmholtz problem for different mesh sizes.

Figure 4 .
Figure 4. Error profiles of the Kelvin-Helmholtz problem.

Table 1 .
Errors and convergence rates of  for the cylindrical explosion problem.

Table 2 .
Relative entropy values for the Kelvin-Helmholtz problem for different mesh sizes.
with  ′  ( +1  ) to yield  1 +  2 = 0, where To simplify  1 , we employ a Taylor expansion to get Taylor expansion to rewrite the second term in (A.2) and simplifying the resulting expression yields the following form for  2 :   )︀)︀  +  +1  (div    )  + Taking the dot product of the above equation with  +1     and using the algebraic identity ( − ) •  = (|| 2 − || 2 − | − | 2 )/2 in the first and second terms yields Here, we have used    =    + (∇   +1 )  to rewrite the pressure gradient term.Multiplying the discrete mass balance (1.35a) by |   | 2 /2, summing it up with (A.4) and simplifying the resulting expression gives us the kinetic energy identity (2.3).