Two-point and interface fluxes

The statement “the discrete entropy cannot grow” is a theorem only under conditions, and the default configuration does not meet all of them.

With the volume fluxes kennedy_gruber and central there is no entropy statement for any interface flux.

What the proof for llf and es_roe covers, and what no proof in Gyre covers:

• It is semi-discrete (exact time integration), for the inviscid terms, on periodic or interior faces, and it assumes positive density and pressure at every node. Positivity is not enforced by the scheme; a state that leaves the admissible set ends the argument. (The steady solver’s line search rejects updates with non-positive density or pressure; a solve that finds no admissible update ends with the reason pseudo-time step collapsed.)
• Wall and far-field boundary conditions are imposed through the interface flux with a constructed exterior state; they carry no discrete entropy proof.
• The viscous terms (BR2 with gradients of velocity and temperature rather than of the entropy variables) and the shock-capturing viscosity carry no discrete entropy proof.
• The steady solver is an implicit pseudo-time iteration; nothing above bounds its iterates.

So the default configuration (ranocha volume flux, roe interface flux) gives entropy-conservative volume terms — aliasing inside an element cannot produce entropy of the wrong sign — and an interface dissipation that removed entropy in every state of the verification study, but not a proven entropy inequality. Choose surface_flux: llf when the proven property matters more than the lower dissipation of Roe.

• consistency: $\mathbf F^*(\mathbf U,\mathbf U;\mathbf n) = \mathbf F(\mathbf U)\cdot\mathbf n$;
• conservation: $\mathbf F^*(\mathbf U_L,\mathbf U_R;\mathbf n) = -\mathbf F^*(\mathbf U_R,\mathbf U_L;-\mathbf n)$, so what leaves one element enters its neighbour even though each face is evaluated from both sides;
• full upwinding in supersonic flow, under a condition on the wave-speed estimates, not on the Mach numbers of the two traces. For flow from left to right, HLLC returns the left flux when its left wave speed $\min(u_{n,L} - c_L,\ \hat u_n - \hat c) \ge 0$, and Roe when the Roe-averaged normal Mach number $\hat u_n/\hat c \ge 1.1$ (hats denote the Roe average). Between $1$ and $1.1$ the slow acoustic eigenvalue $\hat u_n - \hat c$ lies inside the band $|\lambda| < 0.1\,\hat c$ of Harten’s fix, where $|\lambda|$ is replaced by $(\lambda^2/\delta + \delta)/2$, so the Roe flux differs slightly from the upstream flux there. For nearly equal traces the Roe average is the trace state and the condition is the trace Mach number. For very different traces it is not: the Roe-averaged sound speed contains the jump of the tangential velocity, so two traces that are each supersonic along the normal can have a subsonic Roe average, and then neither flux is the upstream flux. The unit test of this property runs at Mach 2.5 with nearby traces.

An evaluated example ($\gamma = 1.4$, normal along $x$; executed by tests/test_docs_content.py): left trace $\rho = 1$, $u = 1$, $v = 3$ at normal Mach number 2, right trace $\rho = 1.3$, $u = 1.1$, $v = -3$ at normal Mach number 1.5. The Roe-averaged normal Mach number is 0.711, and the largest component of the difference to the left flux is $6.7\cdot10^{-2}$ for roe and $5.4\cdot10^{-2}$ for hllc. With $v = \pm 1$ instead (Roe-averaged normal Mach number 1.36) both return the left flux to round-off.

Entropy dissipation is not in this list of properties. For llf and es_roe it follows from the construction; for roe and hllc it is an observation: in the study below the discrete entropy rate was negative for one perturbed periodic state in 2-D and one in 3-D, each on a warped mesh. That is evidence for those states, not a guarantee for every state.

The second table tests the interface fluxes alone. For a pair of states the entropy produced at a face is $(\mathbf w_R - \mathbf w_L)\T\mathbf F^* - (\psi_R - \psi_L)\cdot\mathbf n$; an entropy-stable flux makes it non-positive for every pair, which is the property the proof rests on. The pairs are random and far apart (densities and pressures over a decade, Mach numbers up to about three), far outside what a resolved computation presents to a face.

Let $\mathbf A = \partial(\mathbf F\cdot\mathbf n)/\partial\mathbf U$ with eigenvalues $\boldsymbol\Lambda = \operatorname{diag}(u_n - a,\ u_n,\dots,u_n,\ u_n + a)$ and right eigenvectors $\mathbf R$. Roe’s solver subtracts $\tfrac12\mathbf R|\boldsymbol\Lambda|\mathbf R^{-1}(\mathbf U_R - \mathbf U_L)$ from the mean of the two fluxes; nothing guarantees the sign of the entropy this produces. Barth’s observation is that the eigenvectors can be scaled so that $\mathbf R\,\mathbf T\,\mathbf R\T = \partial\mathbf U/\partial\mathbf w$ with a positive diagonal $\mathbf T$. Writing the jump in entropy variables then gives a dissipation

for any positive mean state at which $\mathbf R$, $\boldsymbol\Lambda$ and $\mathbf T$ are evaluated. Subtracting $\tfrac12\mathbf D$ from an entropy-conservative flux therefore gives a flux whose entropy production at a face is non-positive: the same argument as for llf, with a matrix in place of the scalar $\lambda_{\max}$. For nearby states $\mathbf w_R - \mathbf w_L \approx (\partial\mathbf w/\partial\mathbf U)(\mathbf U_R - \mathbf U_L)$ and $\mathbf D$ tends to Roe’s $|\mathbf A|(\mathbf U_R - \mathbf U_L)$.

Gyre uses the mean state of Winters, Derigs, Gassner and Walch: logarithmic means $\rho^{\ln}$ and $\beta^{\ln}$ of the density and of $\beta = \rho/(2p)$, the arithmetic mean of the velocity, $\hat p = \{\!\{\rho\}\!\}/(2\{\!\{\beta\}\!\})$, $\hat a^2 = \gamma\hat p/\rho^{\ln}$, $\overline{u^2} = 2\abs{\{\!\{\mathbf u\}\!\}}^2 - \{\!\{\abs{\mathbf u}^2\}\!\}$, $\hat h = \gamma/\big(2\beta^{\ln}(\gamma - 1)\big) + \tfrac12\overline{u^2}$, and

The tangential waves are summed without constructing tangent vectors, using $\sum_k \mathbf t_k\mathbf t_k\T = \mathbf I - \mathbf n\mathbf n\T$, so one kernel serves two and three dimensions. The acoustic eigenvalues carry the same Harten fix as roe, which only increases $|\lambda|$ and so does not affect the inequality. With this mean state a contact discontinuity at rest is left exactly alone (no mass flux through it), where llf diffuses it with the speed of sound; unit tests check the inequality for random pairs, the limit $|\mathbf A|$, the symmetry of the operator and the stationary contact.

Inside an element the polynomial resolves the solution, and adding dissipation there would cost accuracy. Jumps at interfaces measure under-resolution: they are $\mathcal O(h^{p+1})$ for resolved features and large only where the mesh cannot represent the solution. Interface dissipation proportional to the jump is therefore automatically high-order where the solution is smooth and stabilising where it is not, while the entropy-conservative volume term guarantees that the volume terms — where aliasing errors arise — neither produce nor destroy entropy (given positive density and pressure).