What You Give Up the Exact Riemann Solution For — the Roe Solver and the Entropy Fix#

The 1-D Euler Riemann problem has an exact solution. One nonlinear equation for the star pressure, one Newton loop, done. Production compressible codes still skip that solve at every face. The answer is in hand. It gets thrown out. Cost. Past an ideal gas, the exact solution may not exist at all.

Move the left and right density and velocity sliders. The sqrt(rho) weights and the three wave speeds $\hat{u}-\hat{c}$, $\hat{u}$, $\hat{u}+\hat{c}$ recompute on the face.

When $\hat{u}-\hat{c}<0<\hat{u}+\hat{c}$, the wave fan straddles the face. Subsonic. Push both velocities hard to one side and the three waves share a sign. The state goes supersonic.

Why the Exact Solution Vanished From Practice#

An exact Riemann solver (Godunov) finds the root of a nonlinear equation at every face. With millions of cells, that root-finding loop runs millions of times per step. Worse, the exact solution is bound to the ideal-gas equation of state. Move to a real gas or a two-phase mixture and the exact solution no longer exists.

Approximate Riemann solvers sidestep this. They linearize the nonlinear Riemann problem locally. One set of algebraic formulas gives the flux with no iteration. The trade is a small loss of exactness.

Collapsing the Jacobian Into a Single Average#

Roe's idea is simple. Approximate the flux difference across a face by a single constant matrix $\hat{A}$ acting on the two states $q_L$, $q_R$.

Here $q$ is the conserved-variable vector and $f$ is the physical flux. This condition is decisive. If the two states are connected by a single wave (a shock or a contact), then $\hat{A}$ propagates that wave exactly, even at finite amplitude. So the solver is approximate, yet exact in front of a pure shock or contact.

The catch is how to choose $\hat{A}$. A careless average breaks the condition above.

Where the sqrt(rho) Average Comes From#

Roe's answer is an average weighted by the square root of density.

Here $u$ is velocity, $H$ is total specific enthalpy, and the carets mark face-averaged values. The Roe sound speed follows as $\hat{c}=\sqrt{(\gamma-1)(\hat{H}-\hat{u}^2/2)}$.

Why the square root of density? Write the state $q$ and flux $f$ in the components of the parameter vector $W=\sqrt{\rho}\,(H,u,1)^T$ and both become perfectly quadratic. Being quadratic, the derivatives of $q$ and $f$ with respect to $W$ are linear, and the path integral that builds $\hat{A}$ evaluates exactly. What falls out is precisely this sqrt(rho)-weighted average.

The weight denominator is $\sqrt{\rho_L}+\sqrt{\rho_R}$. If one density is four times the other, its square root is twice as large, so that side's velocity and enthalpy pull twice as hard on the average. An arithmetic mean does not satisfy $\hat{A}(q_R-q_L)=f(q_R)-f(q_L)$ on a finite jump. Without sqrt(rho) the wave speeds are already wrong in front of a pure shock. The sound speed is $\hat{c}=\sqrt{(\gamma-1)(\hat{H}-\hat{u}^2/2)}$. At $\gamma=1.4$, $\gamma-1=0.4$. What remains after subtracting the kinetic energy $\hat{u}^2/2$ from the total enthalpy is the square of the sound speed. The place where the three wave speeds change sign is the same three-wave split as the Euler eigensystem on an arbitrary normal.

Split Into Waves, Then Add Them Back#

Assembling the flux takes three steps. First decompose the jump $\Delta q = q_R - q_L$ into a sum of the three eigenvectors $\hat{K}_k$; the size of each component is the wave strength $\alpha_k$. Then advect each wave upwind according to the sign of its eigenvalue $\hat{\lambda}_k$.

The first term is a centered average; the second is upwind dissipation weighted by the eigenvalue magnitudes. That structure ports straight to code. The test is the Shu–Osher problem: a Mach 3 shock races into a sinusoidal density field and leaves high-frequency structure behind it. It is exactly where Roe's low numerical dissipation earns its keep.

import numpy as np
 
gamma = 1.4
 
def phys_flux(U):                       # conserved U=(rho, rho*u, E) -> physical flux
    rho = U[0]; u = U[1] / rho; E = U[2]
    p = (gamma - 1) * (E - 0.5 * rho * u * u)
    return np.array([rho * u, rho * u * u + p, u * (E + p)])
 
def roe_flux(UL, UR, delta):            # Roe approximate Riemann flux (+ Harten entropy fix)
    rhoL, rhoR = UL[0], UR[0]
    uL, uR = UL[1] / rhoL, UR[1] / rhoR
    pL = (gamma - 1) * (UL[2] - 0.5 * rhoL * uL * uL)
    pR = (gamma - 1) * (UR[2] - 0.5 * rhoR * uR * uR)
    HL = (UL[2] + pL) / rhoL; HR = (UR[2] + pR) / rhoR
    sL, sR = np.sqrt(rhoL), np.sqrt(rhoR)      # sqrt(rho) weights
    u = (sL * uL + sR * uR) / (sL + sR)        # Roe-averaged velocity
    H = (sL * HL + sR * HR) / (sL + sR)        # Roe-averaged enthalpy
    c = np.sqrt((gamma - 1) * (H - 0.5 * u * u))   # Roe-averaged sound speed
    rr = sL * sR
    drho, dp, du = rhoR - rhoL, pR - pL, uR - uL
    alpha = np.array([(dp - rr * c * du) / (2 * c * c),   # wave strengths alpha_k
                      drho - dp / (c * c),
                      (dp + rr * c * du) / (2 * c * c)])
    lam = np.array([u - c, u, u + c])          # eigenvalues (wave speeds)
    K = np.array([[1, u - c, H - u * c],       # right eigenvectors
                  [1, u,     0.5 * u * u],
                  [1, u + c, H + u * c]])
    al = np.abs(lam)
    small = al < delta                         # entropy fix: floor on |lambda|
    al[small] = (lam[small] ** 2 + delta ** 2) / (2 * delta)
    diss = (al * alpha) @ K
    return 0.5 * (phys_flux(UL) + phys_flux(UR)) - 0.5 * diss
 
def run_shu_osher(N=400, tmax=1.8, cfl=0.4, delta=0.1):
    x = np.linspace(0, 10, N); dx = x[1] - x[0]
    rho = np.where(x < 1, 3.857143, 1 + 0.2 * np.sin(5 * x))   # shock + sinusoidal density
    u   = np.where(x < 1, 2.629369, 0.0)
    p   = np.where(x < 1, 10.33333, 1.0)
    U = np.array([rho, rho * u, p / (gamma - 1) + 0.5 * rho * u * u])
    t = 0.0
    while t < tmax:
        r = U[0]; v = U[1] / r; pp = (gamma - 1) * (U[2] - 0.5 * r * v * v)
        dt = cfl * dx / np.max(np.abs(v) + np.sqrt(gamma * pp / r))
        dt = min(dt, tmax - t)
        F = np.zeros((3, N + 1))
        for i in range(1, N):
            F[:, i] = roe_flux(U[:, i - 1], U[:, i], delta)
        F[:, 0] = phys_flux(U[:, 0]); F[:, N] = phys_flux(U[:, N - 1])
        U[:, 1:N - 1] -= dt / dx * (F[:, 2:N] - F[:, 1:N - 1])
        t += dt
    return x, U[0]
 
x, rho = run_shu_osher()
print(f"t=1.8  min rho={rho.min():.3f}  max rho={rho.max():.3f}")   # -> min~0.81  max~4.08

Roughly forty lines make a complete compressible solver. No Newton iteration, no exact Riemann solve. delta is the entropy-fix parameter covered next.

The Shu–Osher initial field is $\rho=3.857143$, $u=2.629369$, $p=10.33333$ for $x<1$, and $\rho=1+0.2\sin(5x)$, $u=0$, $p=1$ on the right. $N=400$ on $[0,10]$, CFL 0.4, $t=1.8$, $\delta=0.1$. The printout is a density minimum near 0.81 and a maximum near 4.08. The left density 3.857 rises a little toward 4.08, and a trough of the sine falls toward 0.8. That range, with no Newton loop, is what this run holds.

Roe Breaks the Entropy Condition#

The Roe solver treats every wave as a jump. Even a rarefaction is approximated as a train of small shocks. Usually that is harmless. But when a rarefaction contains a sonic point, things change. At that point an eigenvalue $\hat{\lambda}_k$ changes sign, passing through $0$.

When $|\hat{\lambda}_k|\to 0$, the upwind dissipation for that wave disappears. The scheme locks in a stationary expansion shock instead of opening a smooth fan. It satisfies the Rankine–Hugoniot condition but violates the entropy condition — an unphysical solution.

The cure is Harten's entropy fix. It puts a floor under the eigenvalue magnitude near $0$.

At $|\hat{\lambda}_k|=0$ the fix formula equals $\delta/2$, because $(0+\delta^2)/(2\delta)=\delta/2$. If $\delta$ is 0 the floor is 0 and the upwind dissipation vanishes. The code's delta=0.1 lays 0.05 on that floor. Five to ten percent of $(|\hat{u}|+\hat{c})$ is the width before a contact smears too far.

Raise the $\delta$ slider from 0 on the Burgers expansion $u_t+(u^2/2)_x=0$ with $u_L=-0.5$ and $u_R=1.0$. Watch whether the blue kink at $x=0$ moves. The sonic point $u=0$ sits in the middle of the interval.

At $\delta=0$ a blue kink — the stationary expansion shock — sits at $x=0$. Raise $\delta$ and the numerical curve slides down onto the amber exact rarefaction fan. Five to ten percent of the local $(|\hat{u}|+\hat{c})$ is usually safe. Too large, and contacts smear out.

Cheaper: HLL and HLLC#

If Roe feels heavy, keep fewer waves. HLL (Harten–Lax–van Leer) keeps only the two acoustic waves. It drops the middle contact wave.

Here $S_L$, $S_R$ are estimates of the outermost left and right wave speeds. HLL is robust and preserves positive density well. But it cannot keep contacts sharp — there is no middle wave.

The C in HLLC stands for the central Contact wave. It restores the wave HLL discarded: three waves, four constant-state regions. The three solvers sit on a spectrum of wave count.

The production default is usually HLLC. It keeps contacts alive while making density and pressure positivity easy to enforce. Roe has excellent resolution, but the entropy fix is mandatory, and you must watch for the carbuncle phenomenon on strong grid-aligned shocks.

Leave δ at zero and the expansion shock stays on the face#

The sqrt(rho) weight comes from the parameterization $W=\sqrt{\rho}\,(H,u,1)^T$ that makes $q$ and $f$ quadratic. Without that path integral the jump condition on $\hat{A}$ breaks at finite amplitude. Leave $|\hat{\lambda}_k|\to 0$ untouched at a sonic point and Rankine–Hugoniot still holds while entropy does not. Roe without Harten's floor locks the expansion shock as a stationary solution.

How many waves to keep is the next decision. HLL keeps two acoustic waves and smears the contact. HLLC restores the contact and four constant states. Roe goes to five waves in 3-D, but the entropy fix is mandatory, and a strong shock aligned with the grid still needs a carbuncle check. At low Mach the artificial-diffusion scales split on their own. The Hope–Collins note that splits low-Mach Roe into three scalings is that split.