DISCRETE ENTROPY INEQUALITIES VIA AN OPTIMIZATION PROCESS

. The solutions of hyperbolic systems may contain discontinuities. These weak solutions verify not only the original PDEs, but also an entropy inequality that acts as a selection criterion determining whether a discontinuity is physical or not. Obtaining a discrete version of these entropy inequalities when approximating the solutions numerically is crucial to avoid convergence to unphysical solutions or even unstability. However such a task is difficult in general, if not impossible for schemes of order 2 or more. In this paper, we introduce an optimization framework that enables us to quantify a posteriori the decrease or increase of entropy of a given scheme, locally in space and time. We use it to obtain maps of numerical diffusion and to prove that some schemes do not have a discrete entropy inequality. A special attention is devoted to the widely used second order MUSCL scheme for which almost no theoretical results are known.


Introduction
Many physical phenomena can be described by means of a hyperbolic system, also called a system of conservation laws.Some famous hyperbolic systems include the Lighthill-Whitham-Richards or Aw-Rascle models for traffic flow, the shallow water equations of Barr de Saint-Venant, the Euler equation for fluid dynamics and the inviscid magnetohydrodynamics equation.
This class of partial differential equations (PDEs) does not contain any regularization terms such as diffusion or dispersion.Their solutions typically develop discontinuities in finite time.These discontinuities are observed in traffic jams, during floods caused by dam breaks, at hydraulic jumps, or in aeronautics.The PDE should be understood in the weak sense to allow such discontinuous solutions.Doing so makes it possible to construct infinitely many discontinuous solutions for the same initial data.An additional criterion should consequently be included to select only the physical weak solution.It generally takes the form of an entropy (or energy) inequality and is related to the second law of thermodynamics.
Discretizing the PDE to obtain a numerical approximation of the solution can be done in several ways.In this paper we focus on finite volume schemes which are well adapted to the low regularity of the solution and built around the idea of conservation laws.In the design of such schemes, it seems important that the entropy also decreases at the numerical level.This condition ensures that the scheme will not converge towards a nonphysical solution of the PDE.The loss of entropy in each cell during one time step is called the numerical diffusion.Discrete entropy inequalities have mainly been obtained for first order schemes in space and time.Realistic codes use high-order discretizations and splitting techniques and usually incorporate ideas and knowledge from the well-understood first-order framework.For this reason, they probably verify a discrete entropy inequality in most cases.However, no explicit formulas are known in practice.
This paper proposes to quantify the numerical diffusion with an a posteriori minimization technique where the scheme is used as a black box.Our primary goal is to obtain maps of numerical diffusion which quantify in space and time the loss of entropy coming from the choice of discretization.As a practical motivation and a long term objective, let us mention numerical oceanic circulation models where the numerical diffusion is related to artificial changes of salinity, density and temperature between two adjacent distinct water masses, and is usually called spurious mixing.This phenomenon is identified as a major issue in numerical cores for climate application [13].We refer the reader to [5,18] for the quantification of spurious mixing in simplified configurations and [17] in a realistic setting.
This paper provides a mathematical insight on the quantification of numerical diffusion in realistic codes but is still far away from oceanic applications.Maps of numerical diffusion are obtained by minimizing a functional which takes into account the consistency of the numerical entropy fluxes and the fact that the entropy should decrease at each time step.This minimization couples all the cells of the mesh, but we also propose a local and cheap quantification that gives qualitatively good results.A different perspective allows us to construct the worst initial data in terms of entropy by minimizing a different but related functional.We apply this procedure to prove that no discrete entropy inequality exists for most of the versions of the widely used MUSCL approach with a 2 steps Runge-Kutta time discretization.This was suspected in [2].A notable exception is the limitation in the entropy variables with a HLL first order scheme [1].
We now present the mathematical framework on hyperbolic PDEs and their discretization with finite volume schemes, with an emphasis on discrete entropy inequality which is at the core of this work.Our main results and the outline of the paper are described hereafter.

Fundamentals on discrete entropy inequalities
Consider a hyperbolic system in 1 dimension (1D) in space    (, ) +    ( (, )) = 0,  ∈ R,  ∈ R + , where the vectorial unknown  belongs to some convex domain Ω ⊂ R  .The flux  : Ω → R  is  1 -regular and its Jacobian matrix  is diagonalizable with real eigenvalues.We are only interested in weak solutions of (1) that additionally satisfy the entropy inequality (or energy inequality) where the entropy  : Ω → R is strictly convex.The entropy flux  is linked to the entropy  through the relation on their Jacobian matrices  = .Such hyperbolic systems arise in particular in the modeling of nonviscous flows.In this paper, we consider the scalar ( = 1) Burgers equations related to the Lighthill-Whitham-Richards model for traffic flows and the Euler equations of inviscid gas dynamics for which  = 3.We now turn to the numerical discretization of (1) with a finite volume technique.For the sake of simplicity a space interval [, ] with periodic boundary conditions is considered throughout the paper.A finite number of cells  is fixed with size ∆ =  −   and we note  −1/2 =  + ( − 1)∆ for  ∈ 1,  .The points  =  1/2 and  =   +1/2 are identified and all the space subscript are considered modulo  .For the sake of simplicity, every family () ∈ 1, of quantities indexed by the set of cells will be denoted by (  )  .
We also consider a discretization in time 0 =  0 <  1 < . . .  < . . .A Courant-Friedrichs-Lewy (CFL) condition is imposed at each time step.It reads, for some CFL  ∈ (0, 1) depending on the scheme, where is the spectral radius of the Jacobian matrix  (︀

𝑗
)︀ .The time step varies with the time iteration but for the sake of simplicity, we will denote it independently of  as  +1 −   = ∆.
A finite volume scheme writes where the vectors    and   +1/2 correspond to the numerical approximations of the mean of the exact solution  and flux  () Equations ( 1) and ( 2) are understood in a weak sense to allow discontinuous solutions.Inequality (2) does not hold for every discontinuity but selects only entropy satisfying shocks.The fact that the scheme ( 4) is stable and computes the entropy solution, with only entropy satisfying shocks, is strongly related to the existence of a numerical counterpart of ( 2) at the discrete level, called a discrete entropy inequality In the spirit of Definition 1, we introduce the notion of entropy satifying scheme.
Definition 2. Consider a consistent finite volume scheme of   ∈ N cells to the left and   ∈ N to the right.It is a consistent entropy satisfying scheme if there exists a numerical entropy flux function  such that the following two conditions are satisfied.

. , 𝑢) .
There exist several first order schemes (4) with a stencil   =   = 1 for which an explicit formula for  yielding to nonpositive diffusion    is known: -for scalar equations  = 1, monotone schemes, see [14]; -for hyperbolic systems  ≥ 2 the Godunov and HLL schemes are entropy satisfying; see [14].In the specific case of gas dynamics, the HLLC scheme and some relaxation or kinetic schemes are also entropy satisfying see [3,20] and references therein.
For schemes of order larger than 1 the specific form of (5) seems out of reach for hyperbolic systems and the question is still largely open.Some works present results in that direction, mainly for second order schemes.
-In [4,9] and [1], inequality ( 5) is slightly modified.The schemes are either difficult to implement or there is no guarantee that they capture entropy solution.-The local discrete entropy inequality ( 5) is replaced by the weaker global condition )︀ in [10] for the multilayer shallow water equations, in [20] and [16] for gas dynamics and in [6] for a conservation law with nonconvex flux.
-A different approach consists in using a second order scheme and to go back to first order if (5) does not hold.
The MOOD technique (see [7,11]) was initially developed to ensure the positiveness of some quantities like the density and the pressure, as well as some discrete maximum principle.This method was later extended in [2] to ensure the discrete entropy inequalities (5) hold.However it is limited to the gas dynamics (23).
On the other hand many schemes of order 2 or more are employed in the applications for their realistic results.They are typically designed to be of high order when the solution is smooth and include corrections to ensure stability, such as limiters or explicit numerical diffusion.Positivness and lack of spurious oscillations are also often taken into account.However, the existence of discrete entropy inequality often remains an open question, as is the case for the Piecewise Parabolic Method (PPM) [8] or the (Weighted) Essentially Non Oscillatory method (ENO/WENO) [15,21,26].
In this work, we are concerned with the a posteriori determination of the numerical entropy fluxes   +1/2 and numerical diffusions    of Definition 2. We do not attempt to find an analytical choice of .Instead, we follow a different approach and construct an adequate functional  that is minimized to obtain and the corresponding a posteriori numerical diffusion The first objective is to be able to quantify the numerical diffusion in settings where finding an explicit formula for  is out of reach.The functional is built so that the two constraints of Definition 2, namely the nonpositivity of the diffusion and the consistency of the fluxes, are satisfied "as much as possible" and can be applied to any explicit finite volume scheme satisfying Definition 1.This a posteriori approach offers a novel method to construct numerical entropy fluxes.There are strong links between Definition 2 of entropy satisfying scheme and the fact that the minimum of the functional  is 0. Our main results are summarized in the diagram below.
Existence of an explicit  such that ∀ (   )  , ∀ ∈ 1,  ,    ≤ 0 (Defs. 2 and ( 6)) The minimization always reaches its lowest bound 0 For this specific (   )  , for all  ∈ 1,  , The scheme is not entropy satisfying: it is impossible to find an explicit  such that Definition 2 holds The paper is organized as follows.The construction of the functional  and the minimization procedure are introduced in Section 2. In Section 3, we prove several theoretical results about the minimizers of the functional  , including the two first implications of the previous diagram.The first one ensures that the minimization procedure is efficient on the class of entropy satisfying schemes where an expression for the numerical entropy fluxes  +1/2 in (5) can be found in the literature.The second one is more practical, since it is applied on a specific initial data for a generic scheme and allows to visualize the distribution of numerical diffusion by making the best choice for  +1/2 .We also prove a Lax-Wendroff theorem adapted to our minimization procedure in the sense that if a numerical scheme converges and we can find numerical entropy fluxes such that )︂ = 0, then the limit is an entropy solution.
The minimization procedure is illustrated on various testcases in Section 4 both for the Burgers equation and for the Euler equations of gas dynamics.The third implication of the previous diagram is nothing but the contraposition of the first one.We exploit it in Section 5 to build initial data that cannot satisfy any discrete entropy inequality.As an illustration, we prove that most variants of the MUSCL scheme with a second-order Runge-Kutta time discretization for the gas dynamics are not entropy satisfying.

Minimization procedure
The minimization is carried out at some fixed time iteration  which is referred to as "the initial data".In other words is fixed by the user and  +1  is then obtained with Scheme (4).The superscript  plays no particular role and one can fix, e.g.,  = 0.
In this section, the numerical entropy flux is computed by a minimization procedure (7).The function  : R  → R is defined by the sum of two contributions  =   +   .The first part   gathers the contributions of positive numerical diffusion Indeed the quantity )︁ corresponds to the numerical diffusion (6) in the -th cell if the numerical entropy fluxes are given by  +1/2 =  +1/2 .A negative value is consistent with the Definition 2 of an entropy satisfying scheme; thus only positive contributions are kept in the function   .Overall, we have   (︀ ; (   )  )︀ = 0 if and only if Inequality (5) holds with the choice  = .It remains to take into account the consistency property of the numerical entropy fluxes of Definition 2, which is the role of the second part of the functional.We make use of a priori consistency bounds on the numerical entropy fluxes   +1/2 ≤   +1/2 ≤   +1/2 depending on   cells on the left and   on the right.These bounds are defined in Lemma 1 and satisfy the following consistency property This motivates the choice In the particular case where vanishes if and only if   +1/2 =  (), meaning that  is consistent.The function  is nonnegative,  1 -regular and convex with respect to , but the set where it vanishes has no reason either to exist or to be reduced to a single point.

𝑘
. It yields bounds on the individual numerical entropy flux   +1/2 that can be constructed from the initial data (   )  and depends on the chosen numerical scheme.
Proof.The flux at the interface  + 1/2 is independent of the values in cells  −   and smaller, and in cells  +   + 1 and larger.Thus, we modify the initial data by extending it by   −  +1 on its left and by   +  on its right ũ,+1/2 which corresponds to (10).We update this modified initial data with the chosen finite volume scheme ℱ and obtain û,+1/2  .
The interest of the modified initial data is that the fluxes at interfaces  −   −   + 3/2 (and before) and  +   +   − 1/2 (and after) are given by consistency: Now, as the scheme is entropy satisfying On the other hand the flux at interface  + 1/2 remains unchanged G,+1/2 We eliminate the other numerical entropy fluxes by summation: )︀ )︁ .
We can now bound   +1/2 from above and below and conclude.Remark 1.In the case of two points scheme   =   = 1, Lemma 1 correspond to the notion of "interface entropy inequality" in Definition 2.7 of [3].In this specific case it implies the desired discrete entropy inequality (5) at the cost of a time step ∆ twice smaller Proposition 2.9 of [3].

Main results
In this section, we state some properties of the minimizers of  and revisit some standard concepts of finite volume scheme theory from our a posteriori point of view.

Entropy dissipation and zero minimization
Proposition 1.Consider a finite volume scheme (4) that admits a discrete entropy inequality (5) for some numerical entropy fluxes Proof.This property follows from the definition of the functional.If (5) holds, then )︂ = 0. Then the corresponding scheme verifies The scheme is also consistent in the sense that if Let us insist that this result only shows that there exist numerical entropy fluxes yielding to nonpositive numerical diffusion for the particular choice of initial data (︀    )︀  , which does not mean that the scheme is always entropy satisfying.
which is exactly (5).The term   is also zero, thus is constant equal to    , as well as

Discrepancy between global minimizers
When the set where the functional  vanishes is nonempty, it most likely contains several solutions.We first quantify how close they are from each other.
)︂ = 0, and that the numerical flux ℱ used in (4) is consistent and Lipschitz regular.Then for all , When the scheme is known to satisfy a given discrete entropy inequality and when the solution is smooth, this result shows that the difference between the numerical entropy flux found in the literature and the one returned by the minimization procedure is of order ∆ 2 .Note that the consistency property alone gives an order ∆.
Proof.Suppose that the minimization procedure has two different global minimizers , where By convexity of the entropy , )︁ .
If the numerical flux is Lipschitz regular, the last sum is . A first order expansion at point    of the other contributions is The result then follows from  = .

Lax-Wendroff theorem
One of the main theoretical results about numerical schemes for systems of conservation laws is the Lax-Wendroff theorem, which states that if a numerical scheme converges in a certain sense, then the limit is a weak solution.In addition, if the scheme satisfies relevant discrete entropy inequalities, then the limit is an entropy solution.
The latter statement usually requires the numerical entropy flux to be a continuous and consistent function  of the neighboring approximations.This is not the case in this work, since the numerical entropy fluxes   +1/2 are obtained through a minimization procedure and we do not use a entropy flux function .However, it is possible to adapt the Lax-Wendroff theorem to our framework.This extension somehow means it is possible to replace the classical Definition 2 of entropy satisfying scheme by the requirement "for all initial data (︀    )︀  , there exists  , such that . The latter condition differs on the treatment of the consistency property but leads to an identical Lax-Wendroff theorem.Note that proving such a result has no reason to be easier than finding an explicit function  and is possibly more difficult as  has   +   + 1 variables while  is defined on R  .
Since the time discretization will be important in this section, the objects related to the optimization problem (7) at time   , i.e., with initial data (︀

𝑗
)︀  , will be noted with a superscript .The following Theorem is given on a bounded space domain [, ] with periodic boundary conditions since it is the framework of this paper.It can easily be generalized on the unbounded space domain R.
Furthermore, if for a given entropy pair (, ), for all discretization ∆  = (∆  , ∆  ) and for all  ∈ N, there exists a family Proof.For the sake of simplicity, the subscript  in ∆  will be omitted all along the proof.The proof of the convergence of  Δ toward a weak solution of ( 1) is exactly the same as in the original Lax-Wendroff theorem.

A posteriori quantification of the numerical diffusion: numerical results
We now apply the minimization of the functional  to obtain maps of numerical diffusion in different settings1 .The numerical entropy fluxes  , ±1/2 are the results of the minimization (7) and we define the corresponding a posteriori numerical diffusion according to (6) by In practice the minimization is performed using the blackbox fminunc of MATLAB [23], with the initial guess )︁ .The values for  ,  are often small, in particular on fine meshes or when the solution is smooth.However, these values remain far away from the machine precision estimated at 10 −16 .

Continuous solution of the Burgers equation
In this first test case, we consider the scalar ( = 1) Burgers equation Figure 2. Quantification of the numerical diffusion on Testcase (22) for the entropy satisfying Rusanov scheme and the entropy violating Roe scheme.
The choice of () =  2 is arbitrary; for the scalar Burgers equation any convex entropy  could be considered.The initial data is and we use periodic boundary conditions.The exact solution is piecewise linear, and for  < 2 3 is given by The space interval [−2, 2] is discretized with 50 or 100 cells, and the functional  is minimized at the last iteration at  = 0.4.The parameter  in the Courant-Friedrichs-Levy condition ( 3) is set to  = 0.5.We compared the entropy satisfying Rusanov scheme and the entropy violating Roe scheme with a forward Euler march in time, see the appendix for the definitions of the corresponding fluxes ℱ.The stencil is   =   = 1.The minimization procedure detects the two previous behaviors:  ,  remains nonpositive everywhere for Figure 3. Influence of the time discretization on the numerical diffusion for the MUSCL scheme Testcase (22).
the Rusanov scheme, while it has large strictly positive values localized around the stationary entropy creating shock at the sonic point for the Roe schemes.
The two plots of Figure 3 show the results for the MUSCL second order flux in space based on a Rusanov flux and a minmod limiter.We illustrate here the well-known fact that increasing the order in space without increasing the order in time is detrimental in terms of entropy, whereas increasing both orders is efficient.The reader will find more specific details related to the MUSCL scheme in the following section.Together with a forward Euler march in time, ( 5) is grossly false in the sense that the total entropy increases with time: . Combined with a second order Runge-Kutta time stepping, the a posteriori procedure finds a discrete entropy inequality (5).
To conclude this section, let us illustrate Proposition 3. The Rusanov scheme with an Euler method is entropy satisfying.With the choice of  given in the appendix (A.2) the discrete entropy inequality (5) holds.We denote by    the associated nonpositive numerical diffusion (6).We compare the total discrepancy ∆ ∑︀ ⃒ ⃒  ,  −    ⃒ ⃒ between this quantity and the one obtained through minimization (20) and obtained as expected that is decreases one order faster than the total diffusion ∆ ∑︀ )︀⃒ ⃒ as the mesh becomes finer and finer.Note that this quantity does not depend on the numerical entropy fluxes   +1/2 .

Gas dynamics
We now turn to the Euler equations.The unknown is  = (, , ), where  is the density of the fluid,  is its velocity and  its total energy.This system reads The pressure force  is related to  and  through an ideal gas equation of state )︂ where  ∈ (1,3] .
Both the density and the pressure should remain nonnegative, thus }︂ .
There exists an infinite number of entropy inequalities for these equations, see [1].We consider the classical inequality on the specific entropy The minimization method for quantifying the numerical diffusion does not depend on the number of unknowns , since the quantification of the numerical diffusion only concerns the scalar equation on the entropy evolution (24).We focus on a widely chosen scheme of order two in space and time, namely the Van Leer version of the MUSCL scheme [25].
In this scheme the piecewise constant in space approximation (︀

𝑗
)︀  is replaced by a reconstructed piecewise affine data.For scalar equations  = 1, the reconstruction procedure heavily relies on the fact that the exact solutions of (1) verify a maximum principle property and are total variation diminishing (TVD).A family of functions called limiters is considered to determine the slopes in each cell [14] and allows to keep these features at the discrete level.Both properties are lost for hyperbolic systems  ≥ 2. Limiters are also used but several choices are possible.We investigate the effects of some of them in this section.
It remains to decide how the evolution in time of the reconstructed initial data is performed.A first possibility is to compute the exact flux of (1) by solving generalized Riemann problems.This is only feasible for some particular hyperbolic systems and very costly from a computational point of view.This path is followed in [4] for scalar equations and in [9] for systems.A more convenient approach is to combine the reconstruction in space with a second order method in time (typically a second order Runge-Kutta method (RK2)) to obtain a precise enough approximation of the fluxes.Each individual substep is based on a first order solver.
In any case the discrete entropy inequalities found in the literature for the MUSCL approach differ from (5).The quantity  (︀

𝑗
)︀ is replaced with a linear approximation in equation (1.9) of [4] for scalar equations or with a nonlinear entropy diminishing projection operator in Theorem 2.9 of [9] for systems.In [1],  (︀

𝑗
)︀ is replaced by a convex combination of three terms that not only depends on    but also on   −1 and   +1 , see equations (2.7) and (2.10) of [1].None of these modified entropy inequality is sufficient to prove a Lax-Wendroff theorem.Several numerical simulations [2] even indicate that the MUSCL+RK2 scheme may converge to incorrect solutions.
We first reproduce the Sod tube testcase of Toro [24]  We stick to periodic boundary conditions, so there is another discontinuity at  = −1.It creates a 1-shock, a 2-contact discontinuity and a 3-rarefaction wave.
On Figure 5 we compare the first order HLLC scheme Section 10.4.2 of [24] and the Roe scheme without entropy fix Section 11.2 of [24] with a forward Euler time stepping.Then we consider the second order MUSCL scheme with a RK2 time stepping.The slopes are limited on the primitive variable (, , ) and the underlying first order scheme is the HLLC scheme.We compare the results obtained with a minmod limiter and a superbee limiter.
The a posteriori quantification of the numerical diffusion gives once again positive values near the stationary nonphysical shock created by the Roe scheme.It also detects the overcompressive behavior of the superbee limiter, with a spike of positive numerical diffusion located on the oscillations in density around the central contact discontinuity.The superbee limiter is often too strong and may prevent the scheme from converging, see [2].For second order schemes, the numerical diffusion is located around one or two spikes around each discontinuity.This depends on the initialization of the optimization algorithm, see Figure 8 below.
Then we consider a testcase where the solution does not contain a shock, but only a contact discontinuity.The velocity and the pressure are initially constant, equal to 0.1 and 1.The initial density is The final time is 2 seconds, the CFL number is  = 0.75.On Figure (6), we compare the HLL scheme and the HLLC scheme and find unsurprisingly that the latter is much less diffusive.As usual, we observe this fact  .Densities (first line) and numerical diffusion (second line) for Testcase (26).The discontinuity is a slowly moving contact.
through the higher level of detail captured by the HLLC scheme at the final time.The quantification of numerical diffusion confirms it quantitatively, with lower and more localized  , values.The same holds for their second order extensions with a MUSCL procedure, using either the HLL or the HLLC flux as the underlying first order scheme.The slope limitation is on the conservative variables (, , ) and we use a minmod limiter.

A naive a priori quantification of the numerical diffusion
In this section we turn to the practical consideration of the computational cost of the minimization presented above.The functional   couples all the cells of the mesh because each numerical entropy flux  +1/2 appears in two consecutive terms of the sum.The cost grows quadratically with the meshsize, and becomes larger than the minute for meshes with more than 4000 cells, see Figure 7, left.
This observation disqualifies the use of this procedure in realistic codes were the quantification of    would be useful.In ocean global circulation models (OGCMs) for example, the meshes contain billion of cells and the numerical diffusion under scrutiny here is strongly linked with the question of spurious mixing identified as a key issue in such codes [13].This explains why several works deal with its quantification.Some of them ( [5,18] and others) are directly related to the mathematical theory of discrete entropy inequality in the sense that specific choices of numerical entropy fluxes are imposed in (6) to defined the numerical diffusion    .Here we propose a complementary approach: while these works focus on   +1/2 at the risk of making a suboptimal choice and quantifying    incorrectly, we present a more robust choice of    at the cost of forgetting the numerical entropy fluxes  +1/2 .As a consequence the theoretical results of the previous section are almost completely lost; however the quantitative results are very satisfactory and the computational cost is drastically Figure 7. Computational cost of the minimization procedure (left) and of the a priori guess (28) at the first iteration of Testcase (26).diminished, see Figure 7, right.We believe this fastest quantification could be used for some applications in the future, including the quantification of spurious mixing in OGCM or as a metric for mesh refinement.
In the derivation of the consistency part of the functional and )︁ .
An important point is that    and    are computationally affordable, see right of Figure 7. Indeed,   +1/2 and   +1/2 are computed with the scheme (4) on a small initial data centered around the interface  + 1/2 with   +   cells on its left and on its right, see Lemma 1.We observe that the lower bound ( 27) is particularly interesting for two reasons.First, if the left hand side is positive it indicates that (5) cannot hold on this particular cell.We detail and capitalize on that idea in Section 5. Second, even though it is much lower than  ,  , the variations of    are similar.As a final stage, we perform a naive renormalization by a constant coefficient  and define   ,  =    in such a way that the total amount of numerical diffusion is correct: It yields to the a priori quantification of the numerical diffusion On Figure 8 we compare   ,  and  ,

𝑗
. The a posteriori quantification  ,  depends on the initialization of the minimization procedure and produces narrow spikes of numerical diffusion.The a priori quantification   ,  gives smoother results, with a diffusion spread out the whole length of the discrete discontinuity.The computation of this quantity only requires the bounds   +1/2 and   +1/2 and does not require any optimization.We insist on the fact   ,  is pertinent if one is only interested in the numerical diffusion and not in the numerical entropy fluxes, which are not provided by this approach.Remark 3. One could think that the upper bound D  is also interesting.In particular if this quantity was nonpositive in all cells it would be a strong indication that the scheme is entropy satisfying.However, we observed that D  is in most cases strictly positive, even for first order entropy satisfying scheme.It makes this quantity far less interesting both qualitatively and theoretically.

An entropy stress test for numerical schemes
In this section we go back to the classical question of the existence of discrete entropy inequality (5).We explore further the implication of the a priori consistency bounds   +1/2 ≤   +1/2 ≤   +1/2 .It leads to a novel criterion and an algorithm which allows us to detect that entropy inequality ( 5) is not verified.

Worst initial data in terms of entropy
where  0 1/2 ,  0 1/2 and  0 −1/2 are given by ( 9) with  = 0, then the finite volume scheme does not have a discrete entropy inequality.
Proof.By contraposition, let us assume that the scheme has a discrete entropy inequality.Therefore, for all initial data (︀
Therefore the opposite of (30) is satisfied.This is true in particular for the initial data of the form (29).
As a consequence, to determine if a scheme is entropy satisfying or not, we introduce the new functional )︂ .
For a hyperbolic system of  equations, it has (  +   + 1) unknowns in Ω.Following Proposition 4, if the minimization of ℰ returns a strictly negative result then the scheme is not entropy satisfying.In other words this new optimization constructs the worst initial data in terms of entropy.
Remark 4. The first part of (30) is similar to the definition (27) of  0 0 .We already knew that if this quantity is positive the scheme is not entropy satisfying.This result goes further in that direction by restricting the size of the data to   +   + 1 values and by shifting the optimization procedure on the initial data.   of a counterexample.The results for the two versions of the HLLC schemes are very similar and we include only the version based on the wave-speed estimates.Similarly the Rusanov and HLL schemes behave similarly and we only include the HLL scheme.We used the fmincon of MATLAB [23] to minimize the functional ℰ under the constraints (31).The histograms of the values corresponding to (30) for those counterexamples are given on Figure 11.We only kept counterexamples where this value is greater than −10 −10 which is several order of magnitude larger than the machine precision estimated at 2.2 × 10 −16 .
The influence of the choice of the set variables which are limited is striking on Figure 9, with more counterexamples for the conservative variables (, , ) than for the primitive variables (, , ).There are even fewer counterexamples for the entropy variables (, , ).This hierarchy is less clear on the zoom of Figure 10.The distribution on counterexamples around || = 0 depends on the numerical choices.
The counterexample isolated on Figure 12 has a very small total variation and a null velocity.We checked that it remains a counterexample for smaller and smaller timestep (or equivalently for smaller and smaller CFL number).This may indicate that the limit ∆ → 0 in (5) does not hold for the MUSCL+RK2 scheme based on a HLL Riemann solver with a limitation on the conservative variables.
On the other hand we see many counterexamples with large Mach number ||/ √︁   on Figure 9.The chosen first order scheme is not to blame for the lack of discrete entropy inequality.Indeed in that case the exact flux at interface  + 1/2 is either

𝑗+1
)︀ depending on the sign of the velocity.Most numerical schemes reproduce that.
Eventually Functional (30) remains nonnegative for all the random initializations when the limitation is in the entropy variable (, , ) and the first order scheme is HLL or Rusanov.To explore further if this scheme is entropy satisfying, we relaxed the constraint on the total variation and searched for counterexamples in the much larger domain This search was unsuccessful with 57 000 random initializations.Interestingly it holds for the overlimited version of the MUSCL scheme of [1] but also for the simpler original scheme.It supports the following conjecture.
(1) a MUSCL scheme in space based on a HLL first order scheme with a minmod limiter on the entropy variables (, , ) (2) the RK2 march in time described in the Appendix verifies a discrete entropy inequality for a CFL number  = 0.1.The same type of numerical experiment has been performed for the scalar Burgers equation (21).Three entropy satisfying two points schemes (Rusanov, Osher and Godunov) and three non entropy satisfying schemes (Roe, Lax-Wendroff and Mac Cormack) have been tested.We performed 30 000 optimization of the functional ℰ and found no negative values for the entropy satisfying schemes, and counterexamples for the non entropy satisfying schemes, which supports the following conjecture.
Conjecture 2. Consider the RK2+MUSCL scheme for the Burgers equation ( 21) with a minmod limiter on  and the RK2 march in time described in the appendix.For a Courant number up to 1, this scheme is entropy satisfying when the chosen underlying first order scheme is entropy satisfying, and non entropy satisfying when it is not.)︂ .
The HLL scheme (first order) This scheme is also a simple two waves approximate Riemann solver.We present the scheme for gas dynamic.Following Section 10.3 of [24] [24].With −  =   = (  ,   ) the HLL scheme is the Rusanov scheme.
The HLLC scheme for gas dynamics (first order) This scheme is a commonly used approximate Riemann solver for (23).This solver is exact on some important particular solutions of (23), the stationary contact discontinuities.We implemented two versions of the schemes and obtained similar results.The first one corresponds to Paragraph 10.4.2 of [24].The extremal wave speed are estimated with (10.49) of [24] and the pressure and velocities are constant in the "star region".The second one is based on a pressure estimate and corresponds to Paragraph 10.6, variant 1 of [24].

The Roe scheme (first order)
In the case of scalar conservation laws, the numerical scheme is We refer the reader to [22] and Chapter 11 of [24] for the presentation of the Roe scheme for the Euler equations (23).This scheme does not have a discrete entropy inequality since it preserves entropy violating shock ( (  ) =  (  ),  (  ) >  (  )).In can also produce negative pressure [12], in which case the simulation fails entirely.
The MUSCL scheme (first order) The MUSCL procedure is a commonly used procedure to obtain a second ordre scheme in space, and can be easily combined with a second order time scheme such as a 2 step Runge-Kutta method.It is a 4-points flux with   =   = 2.
For scalar conservation laws, the procedure is the following.where ℱ is a first order two points scheme.
In the case of hyperbolic system with  > 1, this strategy is mimicked componentwise.For the Euler equation (23) we followed the strategy of [1].The reconstruction with a minmod limiter can be applied on the three conservative variables (, , ).It is also common to reconstruct the primitive variables (, , ) or in the )︀ 2

•
A variation of the discrete entropy inequality ( 5) is obtained in [1] but is not sufficient to obtain a Lax-Wendroff theorem.
A second order Runge-Kutta method in time )︀

•
For scalar conservation laws  = 1, the maximum of the wave speeds decreases with time, thus the CFL restriction for the computation of ū is more restrictive than the one for ū.This is not true when  ≥ 2, and we adopt the time evolution of equation (2.12) of [1], with the slight modification that ∆ 2 cannot exceed ∆ 1 .

Figure 1 ..
Figure 1.Modification of the initial data away from the stencil.Consistency can be used to compute the fluxes at interfaces  + 3/2 −   −   and  − 1/2 +   +   and the value of   +1/2 then the solution  satisfies the entropy inequality (2) in the sense of distributions on [, ] × (0, +∞).

Figure 4 .
Figure 4. Total numerical diffusion and total discrepancy between two numerical diffusions on the Testcase (22).

Figure 6
Figure 6.Densities (first line) and numerical diffusion (second line) for Testcase(26).The discontinuity is a slowly moving contact.

Proposition 4 .
Consider a numerical flux ℱ with a stencil of   points on the left and   on the right.If there exists an initial data The scalar quantity  (  ,   ) should be large enough.For scalar equation, we set  (  ,   ) = max (| ′ (  ) |, | ′ (  ) |) and for the equations of gas dynamic (23) we take  (  ,   ) = max we introduce the sound speeds   = √︀   /  and   = √︀   /  and the Roe averages ā = √     + √     √   + √   and v = √     + √     √   + √   .The wave speed are estimates as   = min(  −   ,   −   , v − ā) and   = min(  +   ,   +   , v + ā) and the flux is given by equation (10.21) of [24]he beginning of each time step, the constant value   in cell   = [  − Δ 2 ,   − Δ 2] is replaced by the affine function ↦ →    +    ( −   ) ,where   ∈ R is a slope, determined in such a way that no new extrema are created and the scheme remains total variation diminishing.A common choice is Other limiters are possible, see Section 13.7.3of[24].The MUSCL flux is given by To achieve second order in time one can use a two step Runge-Kutta method in time.Consider a numerical flux ℱ with a stencil of   points to the left and   points to the right.Two substeps are computed If the numerical flux ℱ is entropy satisfying with a numerical entropy flux , the same holds for F with the numerical entropy flux +1 , . . ., ū  +1 , . . ., ū +