PERFECT DERIVATIVES, CONSERVATIVE DIFFERENCES AND ENTROPY STABLE COMPUTATION OF HYPERBOLIC CONSERVATION LAWS

. Entropy stability plays an important role in the dynamics of non-linear systems of hyperbolic conservation laws and related convection-diﬀusion equations. Here we are concerned with the corresponding question of numerical entropy stability — we review a general framework for designing entropy stable approximations of such systems. The framework, developed in [28, 29] and in an ongoing series of works [30, 6, 7], is based on comparing numerical viscosities to certain entropy-conservative schemes. It yields precise characterizations of entropy stability which is enforced in rarefactions while keeping sharp resolution of shocks. We demonstrate this approach with a host of second– and higher–order accurate schemes, ranging from scalar examples to the systems of shallow-water, Euler and Navier-Stokes equations. We present a family of energy conservative schemes for the shallow-water equations with a well-balanced description of their steady-states. Numerical experiments provide a remarkable evidence for the diﬀerent roles of viscosity and heat conduction in forming sharp monotone proﬁles in Euler equations, and we conclude with the computation of entropic measure-valued solutions based on the class of so-called TeCNO schemes — arbitrarily high-order accurate, non-oscillatory and entropy stable schemes for systems of conservation laws.

1. Prologue: On perfect derivatives and conservative differences.Let {u ν = u(x ν )} be a gridfunction based at gridpoints x ν := ν∆x with mesh size ∆x = 1/N.The standard centered differencing, , is a conservative approximation of u x (x ν ) in the sense that ν D ∆x u ν ∆x amounts to the boundary terms, just like u x dx is, namely Moreover, u ν D ∆x u ν is also a conservative approximation of the perfect derivative uu x ≡1 2 (u 2 ) x -indeed, a telescoping sum yields quadratic boundary term.
The latter fact is less obvious.For example, u 3 ν D ∆x u ν is not a conservative approximation of the corresponding perfect derivative u 3 u x ≡ 1 4 (u 4 ) x , since ν u 3 ν D ∆x u ν lacks telescopic cancellation.Instead, we propose to consider the centered differencing This is a second-order approximation of u x (x ν ).It is clearly conservative -D ∆x u ν ∆x amounts to the boundary terms of u 4 ν+1 − u 4 ν / u 3 ν+1 − u 3 ν .Moreover, summation by parts shows that u 3  ν D ∆x u ν is also conservative, The centered differencing (1) is rather unusual, and it raises three questions.
• Is there a general recipe for construction of such conservative differences for any multiplier, F (u)u x ≡ F (u) x ?
• How can we utilize such recipes?This paper provides an answer to these questions, by surveying the series of results in [28,29,30,6,7] and concluding with the development of entropic measure-valued computations in the recent work [8].
2. Systems of conservation laws with entropic extension.We consider onedimensional hyperbolic systems of conservation laws of the form which govern the balance of the conservative variables 1 u(x,t) = (u 1 (x,t), . . .,u n (x, t)) and their fluxes f (u) = ( f 1 (u), . . ., f n (u)) .We consider the cases of a pure Cauchy problem, Ω = R, or the periodic case over the torus, Ω = T; in either case, there are no contribution from the boundaries.The system (2) is hyperbolic in the sense that the n × n Jacobian matrix A(u) := ∂ u f (u) has real eigenvalues.The study of such systems was motivated, to a large extent, by the canonical example of Euler equations.
Euler equations (3) are augmented with yet another conservation law which is expressed in terms of the specific entropy S := ln(pρ −γ ), The last equality follows by formal manipulations of (3), stating the conservation of the entropy η(u) = −ρS in terms of the entropy flux F (u) = −ρqS.This motivates the notion of entropy pairs for general system of conservation laws.
Entropy function.A convex function η : R n → R is an entropy function associated with (2) if there exists an entropy flux F : R n → R such that the following compatibility relation holds (here and below the prime denotes the gradient w.r.t. to specified variable, X (u The existence of such compatible entropy pair allows us to proceed with the following formal manipulation Thus, the pair η(u), F (u) forms a conservative extension of (2), in complete analogy to the conservation of physical entropy in Euler equations (4).The convexity of η = η(•) signifies a non-trivial conserved quantity, beyond the obvious conserved linear combinations c • u.Thus for example, the judicious minus sign in ( 4) is chosen to make the corresponding Euler's entropy, η(u) = −ρS, a convex entropy function of the conservative variables ρ, m and E.
Entropy variables [Godunov (1961) [10], Mock (1980) [23]].Define the entropy variables v ≡ v(u) := η (u).Thanks to the convexity of η(u), the mapping u → v is one-to-one and hence we can make the (local ) change of variables u = u(v), which puts the system (2) in an equivalent symmetric form Here, u(•) and f (•) become the temporal and spatial fluxes in the independent entropy variables, 2 v, and the system (7) becomes symmetric in the sense that the Jacobians of these fluxes are, namely symmetric.Indeed, a straightforward manipulation of the compatibility relation (5) implies Consequently, the Jacobians, H (v) and B(v), are the symmetric Hessians of φ(v), and respectively, ψ(v).The so-called entropy potential flux, ψ(v), will play a significant role in our discussion below.
Multi-dimensional systems.We consider hyperbolic systems of nonlinear conservation laws in several space dimensions, Here, the conservative variables u are balanced by multidimensional fluxes f (u) = f (1) , . . ., f (d) where f ( j ) (u) : R n → R n .The system is hyperbolic if the eigenvalues of the n × n symbol, j ξ j A j (u), A j (u) := ∂ u f ( j ) (u), are real for all real ξ = (ξ 1 , . . ., ξ d ).The corresponding question of an entropic extension amounts to searching non-trivial conservative extensions.Arguing formally along the lines of (6), we conclude that η(•) is an entropy function if pre-multiplication by its gradient, η (u) , preserves the structure of 'perfect gradients': In other words, (η, F) is an entropy pair associated with (9) if a convex entropy η : R n → R and the corresponding entropy flux F = F (1) , . . ., F (d) : R n → R d satisfy the compatibility relations The formal manipulation involving (11) yields the conservation of the non-trivial entropic extension, η(u Strong vs. entropic weak solutions.Strong solutions of (2) are interpreted as pointwise values u(x,t).The generic phenomena associated with these nonlinear equations is the breakdown of their strong solutions at a finite time, after which one must admit weak solutions, [3].These weak solution can be observed in terms of their (sliding) averages -for example, in the one-dimensional case, u(x,t) u(y,t)dy, are governed by the balance law, Weak solutions reflect the balance between their spatial averages on the left and the temporal averages (of fluxes) on the right at all finite scales (∆x, ∆t).In this context of weak solutions, one cannot proceed with the formal manipulations (6), (10) which led to the entropy equality η(u) t + ∇ x • F(u) = 0. Instead, among the many weak solutions, physically relevant solutions are identified as those realized as a vanishing viscosity limit which lead to an entropy inequality [10], [13, §7], [9], A weak solution of ( 9) is entropic if it satisfies the entropy inequality (13) for all admissible entropy pairs (η, F) associated with (9).This notion of entropy solution is the cornerstone for the theory of hyperbolic systems of nonlinear conservation laws.We refer the reader to the pioneering contributions of Lax [14,15], and the comprehensive book [3].
3. Entropy conservative fluxes and entropy stable schemes.We are interested in computation of approximate entropy solutions in terms of their cell averages u ν (t) ≈ u(x ν ,t).To simplify matters, we consider the semi-discrete limit of the onedimensional weak (12), ∆t ↓ 0, We consider the corresponding class of semi-discrete conservative schemes of the form serving as consistent approximations to systems of conservation laws (2).Here, u ν (t) denotes the discrete solution along the grid line (x ν ,t).At the heart of matter are the numerical fluxes, f ν+ 1 2 = f (u ν−p+1 , . . ., u ν+p ), which approximate the differential flux, , and in particular, consistent with the differential flux, The framework of conservative difference schemes (14) was initiated in the seminal paper of Lax & Wendroff [21].
Remark that the numerical flux involves a stencil of 2p neighboring grid values centered at half-indexed gridpoints, f (•, •, . . ., •) , and as such, could be clearly distinguished from the (same notation of) the differential flux tagged at integer indexed gridpoints, f (•) f ν .There are several desirable and often competing properties which are sought in the design of such numerical fluxes; we will list some of them later on.We begin with the question of entropy stability of such schemes.Let (η, F) be an entropy pair associated with the system 2. We ask whether the scheme (14a) is entropy-stable with respect to such a pair, in the sense of satisfying a discrete entropy inequality analogous to the entropy inequality η(u) t + F (u) x 0, Here, F ν+ 1 2 = F (u ν−p+1 , . . ., u ν+p ) is a consistent numerical entropy flux, F (u, u, . . ., u) = F (u).In the particular case that equality holds in (15), we say that the scheme (14) is entropy-conservative.
The answer to this question of entropy stability provided in [28] consists of two main ingredients: (i) the use of the entropy variables and (ii) the comparison with appropriate entropy-conservative schemes.We conclude this section with a brief overview.
Making the changes of variables, u ν = u(v ν ), the scheme (14a) recasts into an equivalent form expressed in terms of the discrete entropy variables with a numerical flux Fix an entropy pair (η, F) associated with (2).We say that the scheme ( 16) is entropy-conservative if the discrete analogue of (6) holds, so that the total amount of entropy is conserved in time, ν η(u ν (t))∆x = ν η(u ν (0))∆x.In other words, we seek entropy-conservative fluxes, denoted f * We conclude that is an entropy-conservative numerical flux if pre-multiplication by the entropy gradient, η (u), preserves the structure of 'perfect differences' in the sense that This raises the question of gridfunctions {X ν }, which admit the form of a 'perfect difference' in the sense of a generic representation, . This is precisely the question of conservative differences raised in the prologue, in complete analogy with preserving the structure of 'perfect gradients' in the differential framework (11).
Expressed in terms of the entropy variables, v ν = η (u ν ), entropy conservation (19) Indeed, the 'if and only if' stated in (20) follows from the fact that the difference between the two terms, on the left and on the right, is a perfect difference.We conclude that the numerical flux f * is a perfect difference.Specifically, the following identity holds, [28] (here and below we use ∆X ν+ 1 2 to abbreviate the jump across the cell, ∆X ν+ where ).This brings us to the following.Theorem 3.1.[Tadmor (1987) [28]].Let (η, F) be an entropy pair associated with the one-dimensional system of conservation laws (2), with the corresponding entropy variables v = η (u).
(i) [Entropy conservative scheme].The difference scheme ( 16) is entropy-conservative so that (17) holds, if its numerical flux, (ii) [Entropy stable schemes].Consider a numeral flux f ν+ 1 2 of the form Here, f * is an entropy-conservative flux and D ν+ 1 2 is any positive definite symmetric matrix.Then the resulting scheme ( 16) is entropy stable so that (15) holds. Proof.
satisfies ( 22) holds then the two terms on the right of ( 21) = 0, and we end with an entropy-conservative scheme so that (17) holds.Moreover, a numerical flux of the form ( 23)) satisfies Thus, the two terms on the right of ( 21) are negative, 0, and we end with an entropy stable scheme so that (15) holds 4. Entropy conservative schemes -examples of scalar conservation laws.We discuss the entropy stability of scalar schemes of the form We begin by noting that in the scalar case, all convex η's are entropy functions, and the corresponding entropy-conservative fluxes are given by [Linear equations] We begin with the linear case u t + au x = 0 and the quartic entropy η(u) = u 4 /4.That is, we seek a "perfect difference", D ∆x u ≈ u x such that both D ∆x u ν and u 3 ν D ∆x u ν are conserved.This is precisely the example discussed in the opening prologue of this paper.In this case, with f (u) = u, the entropy variables := u 3 and entropy potential ψ(u) = 3  4 u 4 yield the second-order perfect differencing (1), which we express in terms of the numerical flux u * Observe that this is indeed a second-order accurate difference approximation Thus, the difference scheme , is a second-order difference approximation of u t + au x = 0 which conserves both u ν (t)∆x and u 4 ν (t)∆x.
We turn to several nonlinear examples.
Example 4.2.[Toda flow] Consider the equation u t + (e u ) x = 0 augmented with exponential entropy pair, (e u ) t + (e 2u ) x = 0.The entropy variable associated with η(u) = e u are (u) = e u , the entropy potential is ψ( , and we end up with the entropy-conservative flux: This leads to the dispersive centered scheme,interesting for its own sake, e.g., [17,20,4] ( 1 2 u 2 ) x = 0, augmented with the quadratic entropy ( Numerical viscosity.We continue our discussion with a focus on the quadratic scale entropy η(u) = 1 2 u 2 where the entropy variables coincide with the conservative variables, = u.According to (22), the entropy-conservative schemes are uniquely determined by the numerical flux f * Recall that ψ (u) = f (u), and a further integration by parts of f * The resulting entropy-conservative scheme then takes the viscosity form [27] The entropy stability portion of Theorem 3.1 can now be restated in the following form, [28].
[28] The conservative scheme is entropy-stable if it contains more viscosity than the entropy-conservative scheme (26), in the sense that The proof is straightforward -the numerical flux associated with ( 27) can be expressed as and entropy stability follows from (23) with The corollary above enables to verify the entropy stability of first-second-order accurate schemes.A host of examples can be found in [29] and we quote here the prototype example of Example 4.4.[Lax-Wendroff viscosity (1960) [21]] We consider the genuinely nonlinear case, where f (u) is, say, convex.A quadratic entropy stability is sufficient in this case, to single out the unique physically relevant solution [2].To see how much viscosity is required in this case, we use the fact that the f is increasing, leading to the upper bound of The resulting viscosity coefficient on the right is the second-order Lax-Wendroff viscosity proposed in [21], It follows that the scheme ( 27),( 28) is entropy stable, 1 2 5. Entropy conservative schemes -systems of conservation laws.Our study of entropy stability is based on comparison with entropy-conservative schemes.
In the scalar case, entropy-conservative schemes are unique (for a given entropy pair).For systems, there are many possible choices for numerical fluxes which meet the entropy conservation requirement (22).In this section we present a systematic derivation of entropy-conservative schemes developed in [29], which enjoys an explicit, closed-form formulation using path integration in phase-space.To this end, we first define the piecewise-smooth path of integration as follows.At each cell, consisting of two neighboring values v ν and v ν+1 , we let r and ending at v n+1 the mapping u → v is one-to-one, the path is mirrored in the usual phase space of conservative variables, u = u ν and ending with Equipped with this notation we turn to our next result.
Theorem 5.1.[Tadmor (2003) [29]] Fix an entropy pair (η, F) and let ψ denote the corresponding entropy flux potential (8b).Then, the difference scheme is a conservative approximation consistent with (2), which is entropy-conservative so that (21) is reduced to the equality equality (17) with entropy-consistent numerical flux F * 5.1.An entropy-conservative approximation of Euler equations.We demonstrate the above approach in the context of Euler equations, (3).We seek a conservative approximation of Euler equations which respects the additional entropy equality, η(u) t + F (u) x = 0, for the entropy pair (η, F) = (−ρS, −ρqS).This is precisely the recipe sought in the prologue -namely, we seek an entropy-conservative flux, f * with no artificial numerical viscosity, so that we end with the precise entropy balance, η(u ν (t))∆x = η(u ν (0))∆x.
Integration in phase-space.An entropy-conservative flux for Euler equations using the recipe of theorem 5.1 was derived in [30].The entropy function η(u) = −ρS induces the entropy variables The corresponding entropy flux potential amounts to Next, we choose an approximate Riemann path in phase-space, v j+1 = v j + j , ∆v ν+ 1 2 r j , where {r j } 3 j=1 are three linearly independent directions along the eigensystem of the Jacobian A(u), { j } 3 j=1 are the corresponding orthogonal system, and {m j } 3 j=1 are the intermediate values of the momentum along the path.We end up with an explicit form of Entropy conservative fluxes -Euler eqs.: Figure 1 shows the computations of entropy-conservative Euler solutions in [30].No artificial numerical viscosity is present: except for the minimal amount of entropy decay (∼ 10 −3 ) due to dissipation of the 4th-order Runge-Kutta time integrator, the simulations in figure 2, reflect the conservation of entropy.

An affordable recipe
, is an entropy-conservative flux.An "affordable" entropyconservative flux for Euler equations was derived in by Ismail & Roe in [12] by clever manipulation of these algebraic relations.Expressed in terms of the normalized , the entropy-conservative flux, f * given by the explicit recipe,

Multidimensional systems of conservation laws.
Well balanced shallow-water equations.We consider the 2D shallow water equations, which govern the motion of shallow-water with height h and velocity field q = (q 1 , q 2 ) , driven by the convective fluxes, f ( j ) = hq j , hq 1 q j + 1 2 gh 2 δ 1 j , hq 2 q j + 1 2 gh 2 δ 2 j , and balanced by the prescribed bottom topography b(x), The entropy function is the total energy, E(u) = 1 2 (gh(h + b) + h|q| 2 ).Observe that the shallow-water fluxes are quadratic in z := (h, ) .This enables a straightforward "affordable" algebraic approach for satisfying the energy conservative compatibility relation (22) ).Here we use the usual indexing of two-dimensional grid-functions attached to grid points x ν, µ := x 1ν , x 2µ .Using the average values, z ν+ 1 2 := 1 / 2 z ν + z ν+1 , one finds the x 1 -entropy-conservative flux [6] Similar expression applies for the conservative flux f (2) * ν, µ+ 1 2 in the x 2 -direction.We end up with the energy conservative scheme These schemes are effective in computing steady solutions of shallow-water -for example, the steady state of lake at rest, H (x) := h(x, •) + b(x) = Const.;q = 0, as well as other equilibria states, q • ∇ x q + g∇ x H = 0 shown in figure 3. Indeed, the energy conservative scheme (32) recovers the precise energy balance, in terms of the energy fluxes F ( j ) (u) = 1 2 hq j |q| 2 + ghH .A simulation of a perturbed two-dimensional lake at rest using 600 × 300 gridpoints.Left column: Energy-conservative scheme (32) with first-order numerical viscosity; right column: Energy-conservative scheme (32) with second-order numerical viscosity (outlined in section 6.2).
6. Entropy stable schemes for systems of conservation laws.We recall that the entropy inequality η(u) t + ∇ x • F(u) 0 is imposed as a stability condition which excludes non-physically relevant shock discontinuities, [15].In particular, the entropy decay follows η(u(x,t 2 ))dx η(u(x,t 1 ))dx, t 2 > t 1 .The question is to quantify the inequality, namely -how much entropy decay will suffice?"physically relevant" entropy decay could be dictated by various mechanisms.We mention the most important two: (i) [Physical diffusion].The canonical example of the conservative Euler equations vs. the entropy decay dictated by Navier-Stokes equations is considered in section 6.1 (ii) [Numerical viscosity].According to theorem 3.1, one can add any amount of numerical viscosity to enforce entropy stability.The goal is to add a judicious amount of vanishing viscosity so that in the resulting scheme admits additional desirable and often competing properties of high-resolution and non-oscillatory behavior.This is the topic of section 6.2 6.1.Physical viscosity -a "faithful" approximation of Navier-Stokes equations.The one-dimensional Navier-Stokes equations (NSe) read The viscosity and heat dissipation of the right (expressed in terms of the temperature θ ∼ p/ρ > 0) dictate an entropy dissipation We abbreviate these NSe writing , where encodes the viscosity and heat amplitudes (λ, µ, κ) .We use the entropy-conservative flux, f * to discretize the convective term on the left of (34) and standard centered differencing for the diffusion term on the right: this yields the semi-discrete scheme We now turn to examine the entropy balance of (35).Pre-multiply by the entropy gradient η (u ν ) = v ν : the entropy-conservative flux f * can be decomposed into conservative sums and perfect differences which amount to We end up with the following entropy balance of (35) -a precise discrete analogue of the entropy balance in NSe derived in [30, theorem 3.6 Here, F ∆     6.2.Numerical viscosity -the class of higher-order entropy stable TeC-NO schemes.Numerical viscosity was introduced by von Neumann in the early days of scientific computation [24], as a basic paradigm to stabilize the computation of shock discontinuities4 Peter Lax is a chief architect who developed this paradigm, [21,15,16].The key question is how to tune the numerical viscosity, D ν+  , in order to satisfy a competing set of desired requirements.In this section we review the recent construction of highly-accurate, non-oscillatory entropic TeCNo schemes [7] based on vanishing viscosities tuned to entropy-conservative fluxes.We begin by setting numerical flux of the form ( 23) where D ν+ 1 2 0 is a matrix viscosity coefficient at our disposal.The corresponding difference scheme reads The expression on the right is a discretized diffusion term of vanishing order ∼ ∆x(Dv x ) x .Recall that according to theorem 3. ), we find the entropy stability statement with entropy flux Next, we need to address the question of accuracy: how to tune the {D ν+ 1 2 }'s to achieve highly accurate scheme, say, order of accuracy p? the entropy-conservative fluxes outlined in sections 4-5 are second-order accurate, f * 2 .Richardson extrapolation enables to upgrade these entropy conservative fluxes of any order [22] To maintain entropy stability while keeping the high-order accuracy one can augment these pth-order entropy-conservative fluxes with a judicial amount of highorder numerical viscosity so that Finally, we need to come into terms with a third constraint of maintaining a nonoscillatory behavior of our scheme, namely, the presence of shock discontinuities should not create spurious oscillations.This enforces viscosity coefficients cannot be too small: at the entropy-conservative limit of D ν+ 1 2 ≡ 0, one observes the spurious oscillations present in the numerical simulations reported in figures 1.This brings us to the class of TeCNO schemes -the Entropy Conservative fluxes based on ENO reconstruction schemes developed in [7], which are able to achieve the three competing properties of (i) entropy stability; (ii) arbitrarily high-order accuracy; and (iii) non-oscillatory.Here is a bird's eye view of this class of schemes.
We begin with given cell averages {u ν (t) = u ν (t)} across the computational cell ] at time-level t.Using the ENO procedure, developed in [11,26], one can reconstruct the pointvalues u + ν+ 1 2 (t) and u − ν+ 1 2 (t) on the left and right edges of the computational cells {C ν }, with high-accuracy so that the jump across each interface is of order u + The ENO reconstruction is performed on the scalar components of the rescaled entropy variables w := R v. The sign property for these rescaled variables ends up with the desired entropy stability 7. Epilogue: Computation of measure-valued solutions.

Kelvin-Helmholtz instability.
There is a rather complete stability theory for one-dimensional systems of conservation laws, whose entropic solutions are realized by the vanishing viscosity limit, u t + f (u ) x = u x x , [1].The situation is different, however, in more than one-space dimension.Consider fore example, the Kelvin-Helmholtz instability -the 2D Euler equations (33) subject to initial state which consists of three layers u 0 (x) = .Nevertheless, no matter how small is, there is a lack of convergence.To this end, set = 0.01, fix ω and keep refining the computational mesh.The density ρ computed in figure 7 does not seem to settle as the number of gridpoints increases  We note that these computations are carried out by highly-accurate, non-oscillatory entropy stable TeCNO schemes [8].Lack of convergence is not numerical instability but rather, a property inherent from the underlying conservation laws.In his 2009 Gibbs lecture, [18], Lax noted that "Just because we cannot prove that compressible flows with prescribed initial values exist doesn't mean that we cannot compute them".The question arises, therefore, what information is encoded in our computations?7.2.Computation of entropy measure valued solutions.Recall that in the passage from strong to weak solutions, we had to give up on the certainty of pointvalues, and instead, one can observe only cell averages.A new paradigm of measure valued solutions for conservation laws was introduced by DiPerna, [5] in which one seeks a family of (weighted) measures {ν x, t } such that the conservation law and its entropic extension read Thus, in entropic measure valued (EMV) solutions we give up the certainty of any one solution (either strong or weak) and instead, we only observe the averages in configuration space.The measure ν x, t quantifies the probability of assigning of value ν x, t , u at a given space-time (x,t) The careful computation of EMVs for 2D Euler equations was carried out in [8].The point we make here is that these faithful computations require high-resolution without artificial numerical viscosity.This is precisely what the TeCNO schemes provide.The use of perfect differences prevails.They enable us to explore the question of what computed quantities are encoded in unstable 2D Euler computations.Indeed, the instability is not due to the computation but is sought to be part of the underlying model.These are shown in figure 7 -the lack of convergence of Kelvin-Helmholtz computations, and the apparent convergence of the pdf realization of their entropy measure valued solutions in figure 8, which was extracted from accurate computation over large ensembles.We refer the interested reader to [8].

j ν+ 1 2 n 2 =
j=1 be an arbitrary set of n linearly independent n-vectors, and let δ j k .Next, we introduce the intermediate states,

Figure 4 .
Figure 4. NSe for Sod's problem: entropy-conservative fluxes with viscosity but no heat conduction; 4000 spatial grid points.

Figure 5 .
Figure 5. NSe for Sod's problem: entropy-conservative fluxes with heat conduction but no viscosity; 4000 spatial grid points.

Figure 6 .
Figure 6.NSe for Sod problem: entropy-conservative fluxes with viscosity and heat conduction; 4000 spatial grid points.
1, any positive definite viscosity coefficient D > 0 will induce entropy stability: repeating the same arguments we used in the derivation of NSe (35)-(36) (with ∆d ν+ 1