Godunov's Method — One Tiny Shock Tube Inside Every Cell

Classify each wave by which side of the star pressure it sits on, and drop any step whose neighboring fans overlap.

Throw a hundred cells at a shock and a first-order upwind scheme still smears it across five of them. In 1959, Sergei Godunov picked the opposite move. He placed a tiny shock tube (a Riemann problem) at every cell interface, evolved its exact self-similar solution for one time step, then re-averaged. This post walks through the five-region fan structure with an exact solver, runs the Lax shock tube in twelve lines of Python, and pushes against the first-order wall that Godunov himself proved.

Move the slider across 1.0. Neighboring fans overlap and a red warning lights; below 1.0 the interface flux stays constant in time.

time step Δt / Δt_CFL = 0.70max wave speed = 1.00

Each cell interface spawns a self-similar Riemann fan. Push the time step past 1 and the fans collide — Godunov requires Δt ≤ Δx / max|λ|.

A lit warning is a step whose self-similar fans have already touched. It is the picture of $\Delta t$ past $\Delta x/\max|\lambda|$ .

1959 — One Shock Tube Per Cell Interface#

Finite-volume methods keep the state piecewise constant inside each cell. That means every interface $x_{i+1/2}$ carries a jump between the left state $\mathbf{q}_i$ and the right state $\mathbf{q}_{i+1}$ . Extend that jump to all of space — left state to $-\infty$ , right state to $+\infty$ — and you have a Riemann problem.

Godunov's move was direct. Solve this little Riemann problem at every interface, let its exact self-similar solution evolve for one step, then take the new cell averages. The compressible Euler eigenstructure (three characteristics $u$ and $u \pm c_s$ ) is handled all at once rather than split into advection plus pressure source.

The Riemann Fan: Five Self-Similar Regions#

A $\gamma$ -gas obeys

\partial_t \mathbf{q} + \partial_x \mathbf{f}(\mathbf{q}) = 0, \quad \mathbf{q} = (\rho, \rho u, \rho E)^\top

with $\rho$ density, $u$ velocity, $E$ total specific energy, and $\mathbf{f}$ the conservative flux.

Given $(\rho_L, u_L, p_L)$ and $(\rho_R, u_R, p_R)$ , the solution depends only on $\xi = x/t$ . The fan that opens at the diaphragm splits the plane into five regions.

Undisturbed left state — the left wave head hasn't reached this point.
The left wave itself — smooth ramp if it's a rarefaction, a jump if it's a shock.
Star-L region — the star pressure and velocity, plus a density tied to the left side.
Star-R region — same star pressure and velocity, a different density.
Undisturbed right state.

The 3–4 boundary is the contact discontinuity: pressure and velocity match across it, only density (and entropy) jump.

The Functional Equation in $p^{\star}$ #

There are four star-state unknowns, but solving for the star pressure alone closes the system. Toro's classical reduction gives

f_L(p^{\star}; \mathbf{q}_L) + f_R(p^{\star}; \mathbf{q}_R) + (u_R - u_L) = 0

where each $f_K$ branches on shock vs. rarefaction.

f_K(p) = \begin{cases} (p - p_K)\sqrt{\dfrac{A_K}{p + B_K}} & p > p_K \;\;\text{(shock)} \\[6pt] \dfrac{2 c_K}{\gamma - 1}\left[\left(\dfrac{p}{p_K}\right)^{\!\frac{\gamma-1}{2\gamma}} - 1\right] & p \le p_K \;\;\text{(rarefaction)} \end{cases}

with $A_K = 2/[(\gamma + 1)\rho_K]$ , $B_K = \frac{\gamma - 1}{\gamma + 1} p_K$ , $c_K = \sqrt{\gamma p_K / \rho_K}$ . The star velocity follows directly:

u^{\star} = \tfrac{1}{2}(u_L + u_R) + \tfrac{1}{2}\bigl[f_R(p^{\star}) - f_L(p^{\star})\bigr].

Newton iteration handles $p^{\star}$ in a handful of steps. Once you have it, sample the state at $\xi = x/t$ based on which side you're on and whether each wave is a shock or a rarefaction. The path from that sample to a numerical flux is written out in the flux that runs from the Riemann problem to Godunov.

The Time-Step Bound — When Fans Touch#

The essential constraint is that adjacent fans must not collide during one step.

\Delta t \le \min_i \frac{\Delta x}{\max_k |\lambda_k^{(i+1/2)}|}

The eigenvalues $\lambda_k$ are $u$ and $u \pm c_s$ . When two fans touch, the self-similar assumption breaks and conservation breaks with it.

The red warning above is this inequality, already broken.

Python: The Lax Shock Tube, Exactly#

Here is Toro's procedure in compact form. Lax tube: $(\rho_L, u_L, p_L) = (0.445, 0.698, 3.528)$ , $(\rho_R, u_R, p_R) = (0.5, 0, 0.571)$ .

import numpy as np
 
GAMMA = 1.4
 
def f_wave(p, rho_K, p_K, c_K, gamma=GAMMA):
    # one-wave pressure-velocity function (branches shock vs rarefaction)
    if p > p_K:
        A = 2.0 / ((gamma + 1.0) * rho_K)
        B = (gamma - 1.0) / (gamma + 1.0) * p_K
        return (p - p_K) * np.sqrt(A / (p + B))
    return (2.0 * c_K / (gamma - 1.0)) * ((p / p_K) ** ((gamma - 1.0) / (2.0 * gamma)) - 1.0)
 
def f_prime(p, rho_K, p_K, c_K, gamma=GAMMA):
    if p > p_K:
        A = 2.0 / ((gamma + 1.0) * rho_K)
        B = (gamma - 1.0) / (gamma + 1.0) * p_K
        return np.sqrt(A / (p + B)) * (1.0 - (p - p_K) / (2.0 * (p + B)))
    return (1.0 / (rho_K * c_K)) * (p / p_K) ** (-(gamma + 1.0) / (2.0 * gamma))
 
def solve_star(left, right, gamma=GAMMA, tol=1e-9, max_iter=80):
    rho_L, u_L, p_L = left
    rho_R, u_R, p_R = right
    c_L = np.sqrt(gamma * p_L / rho_L)
    c_R = np.sqrt(gamma * p_R / rho_R)
    # first guess: pressure average corrected by velocity jump
    p = max(1e-6, 0.5 * (p_L + p_R) - 0.125 * (u_R - u_L) * (rho_L + rho_R) * (c_L + c_R))
    for _ in range(max_iter):
        f  = f_wave(p, rho_L, p_L, c_L) + f_wave(p, rho_R, p_R, c_R) + (u_R - u_L)
        df = f_prime(p, rho_L, p_L, c_L) + f_prime(p, rho_R, p_R, c_R)
        dp = f / df
        p_new = max(1e-6, p - dp)
        if 2.0 * abs(p_new - p) / (p_new + p) < tol:
            p = p_new
            break
        p = p_new
    u = 0.5 * (u_L + u_R) + 0.5 * (f_wave(p, rho_R, p_R, c_R) - f_wave(p, rho_L, p_L, c_L))
    return p, u, c_L, c_R
 
# Lax shock tube
left  = (0.445, 0.698, 3.528)
right = (0.5,   0.0,   0.571)
p_star, u_star, c_L, c_R = solve_star(left, right)
print(f"p* = {p_star:.5f}, u* = {u_star:.5f}")
# -> p* = 3.40948, u* = 0.62556

This twelve-line Newton loop is the beating heart of one Godunov step. The exact solver was eventually replaced by Roe's 1981 linearization and later by HLL and HLLC, mostly for speed. The exact solution still serves as the reference answer when benchmarking new schemes.

Split the printed states. Left pressure 3.528, right 0.571, star pressure 3.40948: down 0.11852 from the left, up 2.83848 from the right. Star velocity 0.62556 sits 0.07244 under the left speed 0.698 and above the right speed 0. The arithmetic mean of the speeds is 0.349. The second term of the star-velocity formula pushes 0.349 up to 0.62556.

The branch test is enough. Star pressure below the left state selects the rarefaction form of $f_L$ ; above the right state, the shock form of $f_R$ . Tolerance $10^{-9}$ , cap 80, floor $10^{-6}$ . With $u_R-u_L=-0.698$ the guess starts above the mean pressure 2.0495 and reaches 3.40948. Five cells on a 100-cell mesh are 5% of the grid, and slider 1.0 is $\Delta t=\Delta x/\max|\lambda|$ .

Open the Fan Yourself#

The viz below re-runs the same Newton iteration in your browser and uses the time slider to open the self-similar fan. Switch the preset to 123 problem (two flows running away from each other into a near-vacuum) and watch the deep pressure trough form in the middle. Strong-blast collapses to almost a single shock.

presettime t = 0.140

Three waves emerge from the origin: a left-going rarefaction or shock, a contact discontinuity at u*, and a right-going wave. The vertical pink dashed line marks the initial diaphragm at x = 0.

Start with a small t and grow it slowly. Left and right waves move out at fixed speeds. In Lax, the left side opens as a smooth rarefaction and the right side jumps as a shock. The tiny density step in the middle is the contact discontinuity.

The First-Order Wall — Godunov's Theorem#

In the same 1959 paper, Godunov proved a sober corollary: no linear, monotone scheme can be more than first-order accurate. His own construction was, by his own theorem, capped.

The detour around that wall is the entire history of 1970s–80s high-resolution shock capturing. van Leer's MUSCL, Harten's TVD, Roe's linearization, Colella and Woodward's PPM, Shu and Osher's WENO — all keep Godunov's skeleton but replace the piecewise-constant subgrid model with a nonlinear reconstruction that preserves monotonicity while raising the order. The Godunov bones stayed.

When the exact fan belongs at the interface#

$p^{\star}=3.40948$ lies under the left pressure 3.528 and over the right pressure 0.571, so the left wave is a rarefaction and the right wave is a shock. $u^{\star}=0.62556$ is 0.07244 under the left speed 0.698. If every interface needs that pair, run the exact fan. If a shock five cells wide on a 100-cell mesh is acceptable, first-order upwind is enough.

The exact fan is expensive. That is why Roe, HLL, and HLLC took the seat after 1981. The reference values that remain are this $p^{\star}$ and this $u^{\star}$ . The opposite choice, smearing a shock on purpose, is the 1950 artificial viscosity that thickened shocks. Godunov captures instead of smearing, and $\Delta t$ stays tied to the fan.