Conservation Laws and Riemann Problems

Pdes

A CONSERVATION LAW is a PDE of the form

\partial_{t} u + \partial_{x} f (u) = 0, u (x, 0) = u_{0} (x),

expressing that a conserved quantity $u$ evolves through a flux $f (u)$ in the $x$ -direction. The integral form — $\partial_{t} \int_{a}^{b} u d x = f (u (a, t)) - f (u (b, t))$ — says that $u$ is neither created nor destroyed, only transported. For NONLINEAR $f$ , smooth initial data can develop SHOCKS in finite time (see method of characteristics), and the classical PDE breaks down. The theory of WEAK SOLUTIONS, the RIEMANN PROBLEM, and modern shock-capturing FINITE-VOLUME methods all start from this point.

Two canonical examples drive the theory:

Inviscid Burgers: $f (u) = u^{2} /2$ . The simplest nonlinear scalar conservation law; develops shocks from smooth data. Used as the testbed for almost every shock-capturing numerical method.
Euler equations of gas dynamics: a system of three conservation laws (mass, momentum, energy) for an inviscid compressible fluid. The foundation of computational fluid dynamics; the standard hyperbolic test that any production CFD code must handle.

Weak solutions

A classical (smooth) solution to $\partial_{t} u + \partial_{x} f (u) = 0$ ceases to exist once characteristics intersect. To handle the post-shock regime, we relax the requirement to WEAK SOLUTIONS: functions $u$ (allowed to have jumps) that satisfy the equation in the integral sense. Test against any smooth compactly-supported function $ϕ (x, t)$ :

\int_{0}^{\infty} \int_{- \infty}^{\infty} (u \partial_{t} ϕ + f (u) \partial_{x} ϕ) d x d t + \int_{- \infty}^{\infty} u_{0} (x) ϕ (x, 0) d x = 0.

A weak solution may have JUMP DISCONTINUITIES at curves $x = x_{s} (t)$ . Across such a discontinuity, the Rankine-Hugoniot condition relates the jump in $u$ to the jump in flux and the shock speed:

s = \frac{d x _{s}}{d t} = \frac{f ( u _{L} ) - f ( u _{R} )}{u _{L} - u _{R}} = \frac{[ f ]}{[ u ]} .

For Burgers' equation specifically, this becomes $s = (u_{L} + u_{R}) /2$ — the shock moves at the AVERAGE of the values on either side.

The catch: Rankine-Hugoniot alone admits MULTIPLE weak solutions for the same initial data. Some are physical (compressive shocks); some are mathematical ghosts (expansive "shocks" that violate the second law). The ENTROPY CONDITION (Lax, Oleinik) selects the physical solution:

For Burgers: a shock is admissible only if u_{L} > u_{R} .

Equivalently: characteristics on BOTH SIDES must run INTO the shock. The physical-shock criterion is the vanishing-viscosity limit: solutions of $\partial_{t} u + \partial_{x} f (u) = ϵ \partial_{xx} u$ as $ϵ \to 0^{+}$ always satisfy the entropy condition. Numerical methods that produce the correct entropy solution are called ENTROPY-CONSISTENT.

The Riemann problem

A RIEMANN PROBLEM is the initial-value problem with piecewise-constant initial data:

u (x, 0) = {u_{L} u_{R} x < 0 x > 0 .

This is the FUNDAMENTAL BUILDING BLOCK of modern numerical methods for conservation laws. For Burgers' equation, the Riemann solution is:

Shock case ( $u_{L} > u_{R}$ ): $u (x, t) = u_{L}$ if $x < s t$ , $u_{R}$ if $x > s t$ , with shock speed $s = (u_{L} + u_{R}) /2$ . Two characteristics converging.
Rarefaction case ( $u_{L} < u_{R}$ ): $u (x, t) = u_{L}$ if $x < u_{L} t$ , $x / t$ if $u_{L} t < x < u_{R} t$ , $u_{R}$ if $x > u_{R} t$ . A smooth FAN of characteristics filling the wedge.

The two cases are RELATED but qualitatively different: shocks are discontinuous, rarefactions are smooth. The entropy condition picks the right one in each case. Notice the SELF-SIMILAR structure: the solution depends only on $x / t$ , a consequence of the scale invariance of Burgers' equation under $(x, t) \to (λ x, λ t)$ .

Godunov's method

The seminal idea (Sergei Godunov, 1959). Discretize $u$ as PIECEWISE CONSTANTS on a finite-volume grid: in each cell $[x_{j - 1/2}, x_{j + 1/2}]$ , the value $u_{j}^{n}$ is the cell average. At each time step, between adjacent cells with values $u_{j}^{n}$ and $u_{j + 1}^{n}$ , the local situation is exactly a RIEMANN PROBLEM with left state $u_{L} = u_{j}^{n}$ and right state $u_{R} = u_{j + 1}^{n}$ . SOLVE that Riemann problem analytically, read off the flux at the cell interface, and update via:

u_{j}^{n + 1} = u_{j}^{n} - \frac{Δ t}{h} (F_{j + 1/2} - F_{j - 1/2}) .

The flux $F_{j + 1/2}$ is the value of $f (u^{*})$ where $u^{*}$ is the Riemann solution at the cell interface. For Burgers:

If $u_{L} > u_{R}$ (shock): the shock moves at speed $s = (u_{L} + u_{R}) /2$ . If $s > 0$ , the interface sees $u_{L}$ and the flux is $u_{L}^{2} /2$ ; otherwise the flux is $u_{R}^{2} /2$ .
If $u_{L} < u_{R}$ (rarefaction): if both sides have the same sign, the flux is the "upwind" one. If they straddle zero, the flux is the SONIC value $f (0) = 0$ .

The Godunov flux is the COARSEST upwind reconstruction — first-order accurate in space but FORMALLY ENTROPY-CONSISTENT because it solves the local Riemann problem exactly. Higher-order schemes (Lax-Wendroff, MUSCL, WENO) extend this by using piecewise-LINEAR or piecewise-POLYNOMIAL reconstruction within each cell while keeping the Riemann-flux idea at cell interfaces.

Code: Riemann problems for Burgers

# Godunov scheme for the inviscid Burgers equation u_t + (u²/2)_x = 0.
# Tests:
#   (1) Shock-tube Riemann problem (u_L = 1, u_R = 0) — exact shock at
#       speed s = (u_L + u_R)/2 = 0.5; shock position at T = 1 is x = 0.5.
#   (2) Rarefaction-wave Riemann problem (u_L = 0, u_R = 1) — exact
#       solution u(x,t) = max(0, min(x/t, 1)), a fan opening at the origin.

import numpy as np

def godunov_flux(uL, uR):
    """Exact-Riemann flux for f(u) = u²/2."""
    if uL > uR:                          # compressive: shock
        s = 0.5 * (uL + uR)
        return 0.5 * uL**2 if s > 0 else 0.5 * uR**2
    else:                                # expansive: rarefaction
        if uL >= 0: return 0.5 * uL**2   # supersonic right
        if uR <= 0: return 0.5 * uR**2   # supersonic left
        return 0.0                       # sonic point at u = 0

def godunov_burgers(u_init, T, cfl=0.4):
    N = len(u_init)
    h = 4.0 / (N - 1)                    # we'll use x in [-2, 2]
    u = u_init.copy()
    dt = cfl * h / max(np.max(np.abs(u)), 1e-6)
    nsteps = int(T / dt)
    dt = T / nsteps
    for _ in range(nsteps):
        F = np.array([godunov_flux(u[i], u[i+1]) for i in range(N - 1)])
        u_new = u.copy()
        u_new[1:-1] = u[1:-1] - (dt / h) * (F[1:] - F[:-1])
        u = u_new
    return u

N = 401
x = np.linspace(-2, 2, N)
T = 1.0

# ─── Test 1: shock-tube Riemann problem ────────────────────────────────
u0_shock = np.where(x < 0, 1.0, 0.0)
u_shock  = godunov_burgers(u0_shock, T)

# Locate the numerical shock as the midpoint between u > 0.95 and u < 0.05
idx_high = np.where(u_shock > 0.95)[0]
idx_low  = np.where(u_shock < 0.05)[0]
shock_num = 0.5 * (x[idx_high[-1]] + x[idx_low[0]])

print(f"(1) Shock-tube Riemann (u_L=1, u_R=0), T = {T}:")
print(f"    Analytical shock position:  x_s = 0.5")
print(f"    Numerical shock position:   x_s = {shock_num:.4f}")
print(f"    Error: {abs(shock_num - 0.5):.4f}")

# ─── Test 2: rarefaction-wave Riemann problem ──────────────────────────
u0_rare = np.where(x < 0, 0.0, 1.0)
u_rare  = godunov_burgers(u0_rare, T)
exact = np.clip(x / T, 0, 1)
exact[x < 0] = 0
err = np.max(np.abs(u_rare - exact))

print(f"\n(2) Rarefaction Riemann (u_L=0, u_R=1), T = {T}:")
print(f"    Exact solution: u(x,t) = max(0, min(x/t, 1))")
print(f"    ||u_num - u_exact||_inf = {err:.4f}")
print(f"    (Godunov resolves the fan; small error is from grid spacing.)")

Output:

(1) Shock-tube Riemann (u_L=1, u_R=0), T = 1.0:
    Analytical shock position:  x_s = 0.5
    Numerical shock position:   x_s = 0.4950
    Error: 0.0050

(2) Rarefaction Riemann (u_L=0, u_R=1), T = 1.0:
    Exact solution: u(x,t) = max(0, min(x/t, 1))
    ||u_num - u_exact||_inf = 0.0513
    (Godunov resolves the fan; small error is from grid spacing.)

Two outcomes. The shock-tube Riemann problem (compressive, $u_{L} = 1 > u_{R} = 0$ ) produces a clean numerical shock at $x = 0.495$ vs the analytical $x = 0.5$ — 0.005 error from 401 grid points over $[- 2, 2]$ , i.e. essentially exactly resolved at this discretization. The shock speed is captured exactly by Rankine-Hugoniot; the small position error is from the finite cell width. The rarefaction problem (expansive, $u_{L} = 0 < u_{R} = 1$ ) produces a smooth fan profile, with max error 0.05 (5%) — first-order Godunov is somewhat diffusive at the fan edges, and higher-order reconstruction (MUSCL with limiters) would resolve this more sharply.

Beyond Burgers: systems of conservation laws

The scalar Burgers theory generalizes to SYSTEMS:

\partial_{t} u + \partial_{x} f (u) = 0, u \in R^{p} .

The Jacobian $A (u) = \partial f / \partial u$ has $p$ real eigenvalues (assuming the system is HYPERBOLIC), each giving a characteristic FIELD. The Riemann problem decomposes into a sequence of $p$ waves — each either a shock, a rarefaction, or a contact discontinuity. The full Riemann solution is the composition of these.

For the EULER EQUATIONS (compressible gas dynamics with conserved variables $(ρ, ρ u, E)$ ), the three characteristic fields are: a sound wave moving LEFT ( $u - c$ with $c$ the speed of sound), a CONTACT discontinuity moving at fluid speed $u$ , and a sound wave moving RIGHT ( $u + c$ ). The Riemann problem is the classical "shock tube": initial data with different gas states on the two sides produces a left-moving rarefaction, a middle contact, and a right-moving shock — the structure observed in shock-tube experiments and the standard CFD test case (the SOD shock tube).

Exact Riemann solvers for the Euler equations are non-trivial (involving a nonlinear iteration to find the pressure in the "star region" between the left and right waves), so most practical CFD codes use APPROXIMATE Riemann solvers — Roe, HLL, HLLC, AUSM — that capture the essential wave structure while being cheap to evaluate at each cell interface. The choice of approximate Riemann solver is a major tuning knob in production CFD codes.

The MUSCL / WENO upgrade

First-order Godunov is robust but diffusive. The MUSCL framework (van Leer 1979) extends it to second order by piecewise-LINEAR reconstruction within each cell, with SLOPE LIMITERS (minmod, van Leer, superbee, MC) that prevent spurious oscillations near discontinuities. WENO (weighted essentially non-oscillatory, Shu 1989) reaches arbitrarily high order via weighted polynomial reconstruction. Both are the standard tools for modern shock-capturing CFD: keep the Riemann-flux idea, just feed it more accurately-reconstructed states.

Connections

Method of characteristics — the analytical foundation; conservation laws are the natural setting for the shock-formation story.
The three classical equations — the wave equation is a LINEAR conservation law ( $f (u) = c u$ ); nonlinear conservation laws are its hyperbolic cousins where the wave speed depends on the solution.
Finite differences — upwind FD is a degenerate case of Godunov-type finite volume.
Finite elements — discontinuous-Galerkin methods unify finite-volume shock capturing with finite-element high-order accuracy.
Method of lines — once spatial discretization is done, time integration of the resulting ODE system follows the usual MOL framework.