The Same Code Gave 0.79 and 2.00 at Once — What Verification Cases Actually Measure

A second-order scheme reported 0.79#

The verification list for the compressible code had three entries. Sod's shock tube, isentropic vortex convection, shock–vortex interaction. I started with Sod.

I refined the grid from 100 to 800 cells in four steps and measured the $L_1$ error in density. The convergence order came out as 0.79. MUSCL reconstruction, HLLC flux, two-stage SSP Runge–Kutta in time. The design order is 2.

I spent two days digging through the reconstruction and the boundary conditions. There was no bug. The 0.79 is correct behavior.

This post shows why, with numbers. It also watches the same code report 2.001 on a different test, and sorts out what each of the three verification cases can and cannot measure.

Convergence order is a property of the solution, not the scheme#

Order proofs come out of a Taylor expansion. Using one requires the solution to be differentiable that many times. At a discontinuity, $\rho'$ does not exist. When the premise breaks, so does the conclusion.

So "second-order scheme" is shorthand for "second order on smooth solutions." Drop the qualifier and the sentence becomes false. Raise N and only the smear width shrinks; press show Linf and the same profile is read on a different slope.

The width of the smeared band narrows with the grid. The height across the jump does not. The $L_1$ slope and the $L_\infty$ slope split on one screen. The numerical solution is modeled as the exact profile smeared over a width proportional to the cell size.

The area a smeared discontinuity leaves in $L_1$#

A shock-capturing scheme resolves a discontinuity by spreading it over a few cells. The number of cells $\alpha$ stays roughly fixed as the grid is refined. The physical smear width is $\alpha h$.

Here $\Delta\rho$ is the density jump across the discontinuity, $\alpha$ the number of cells the smear covers, and $h$ the cell size. The height is fixed and only the width shrinks. The area scales with $h$. That is first order.

Sod's solution carries two discontinuities: the shock and the contact. They smear differently.

A shock has a compressive mechanism that keeps sharpening it. Its smear width stays pinned at a few cells, giving $O(h)$. A contact has no such mechanism. Numerical diffusion accumulates over time and the width shrinks more slowly than $h$. The known estimate for monotone schemes is roughly $O(h^{2/3})$.

The measured 0.79 is the mixture of the two. It sits between 2/3 and 1.

Why $L_\infty$ stopped at 0.083#

Read the very same calculation in the maximum norm and the number stops moving at all.

Somewhere inside the smeared region there is always one cell where the numerical and exact values differ by a fixed fraction of $|\Delta\rho|$. That fraction $\beta$ does not care about the grid.

The measured values were 0.0678 → 0.0830 → 0.0834 → 0.0857. The printed column is $6.7837\times 10^{-2}$, $8.3037\times 10^{-2}$, $8.3428\times 10^{-2}$, $8.5667\times 10^{-2}$, with orders −0.292, −0.007, −0.038. By N=200 the error is already stuck near 0.083, and refining to 800 never pushes it past 0.0857. A smear width pinned to a few cells is the same seat as the artificial viscosity that was written in to blur a shock. There a viscosity term fixed the width. Here the norm fixes the height.

The $L_1$ column moves differently. At N=100 it is $6.5089\times 10^{-3}$, at 200 it is $3.7579\times 10^{-3}$ with order 0.792, at 400 it is $2.1461\times 10^{-3}$ with 0.808, and at 800 it is $1.2420\times 10^{-3}$ with 0.789. It wanders near 0.79 and never clears 1. On the smooth wave, the unlimited $L_1$ on the same four grids is $2.0251\times 10^{-3}$, $5.0589\times 10^{-4}$ (2.001), $1.2630\times 10^{-4}$ (2.002), $3.1562\times 10^{-5}$ (2.001). With minmod on, the column is $5.6574\times 10^{-3}$, $2.0125\times 10^{-3}$ (1.491), $5.7467\times 10^{-4}$ (1.808), $1.5736\times 10^{-4}$ (1.869). At N=100 the ratio of those two errors is the 2.8 already noted. At N=800 the unlimited $3.1562\times 10^{-5}$ and the limited $1.5736\times 10^{-4}$ still sit side by side. Order 1.869 is the climb back toward 2, not a change in the design order.

Same code, same scheme, smooth solution#

The compressible Euler equations also admit a smooth exact solution that is easy to write down. Uniform pressure, uniform velocity, and a perturbation in density alone. It is a pure contact wave, so the shape is carried along unchanged.

Here $\rho_0$ is the initial density profile, $u_0$ the uniform velocity, and $p_0$ the uniform pressure. Run one full period on a periodic domain and compare with the initial field. This is the job the two-dimensional isentropic vortex does.

I convected a Gaussian bump one period at $u_0 = 1$. With the limiter off and a plain central slope, the orders came out 2.001, 2.002, 2.001. The design order shows up untouched. Steepen that bump toward a top hat and the scheme line stays put while the measured order settles toward 1. The slider for that sits right after the printed table.

Where minmod shaves 2.001 down to 1.491#

Production codes run with the limiter on. Minmod zeroes the slope whenever the two differences disagree in sign.

Here $\sigma_i$ is the reconstruction slope in cell $i$. At the top of the bump the left difference is positive and the right one is negative. The test fires and $\sigma_i = 0$. That single cell drops to first-order reconstruction.

The measured orders were 1.491 → 1.808 → 1.869. They are climbing back toward 2. The order did not change; the first-order error from one extremum simply has not washed out yet. On the same grid, turning the limiter off makes the error 2.8 times smaller.

The practical consequence is a single rule. Verify order with the limiter off or with an extremum-preserving limiter, and run production with the limiter on. Mix the two settings in one table and neither row can be trusted.

One solver, four numbers#

I ran Sod and the smooth wave through the same solver. The references are the exact Riemann solution and a translated Gaussian.

import numpy as np
G = 1.4
 
def riemann_star_pressure(rl, ul, pl, rr, ur, pr):
    al, ar = np.sqrt(G*pl/rl), np.sqrt(G*pr/rr)
    def fk(p, r, u, a, pk):
        if p > pk:
            A, B = 2.0/((G+1)*r), (G-1)/(G+1)*pk
            f = (p-pk)*np.sqrt(A/(p+B))
            df = np.sqrt(A/(B+p))*(1.0-0.5*(p-pk)/(B+p))
        else:
            f = 2*a/(G-1)*((p/pk)**((G-1)/(2*G))-1.0)
            df = 1.0/(r*a)*(p/pk)**(-(G+1)/(2*G))
        return f, df
    p = 0.5*(pl+pr)
    for _ in range(80):
        f1, d1 = fk(p, rl, ul, al, pl)
        f2, d2 = fk(p, rr, ur, ar, pr)
        dp = (f1+f2+(ur-ul))/(d1+d2)
        p = max(p-dp, 1e-10)
        if abs(dp) < 1e-13: break
    f1, _ = fk(p, rl, ul, al, pl)
    f2, _ = fk(p, rr, ur, ar, pr)
    return p, 0.5*(ul+ur)+0.5*(f2-f1)
 
def sample_exact_riemann(x, t, rl, ul, pl, rr, ur, pr):
    ps, us = riemann_star_pressure(rl, ul, pl, rr, ur, pr)
    al, ar = np.sqrt(G*pl/rl), np.sqrt(G*pr/rr)
    s = x/t
    rho = np.empty_like(s)
    for i, si in enumerate(s):
        if si <= us:                                  # left of contact
            if ps > pl:                               # left shock
                rsl = rl*((ps/pl+(G-1)/(G+1))/((G-1)/(G+1)*ps/pl+1))
                ss = ul-al*np.sqrt((G+1)/(2*G)*ps/pl+(G-1)/(2*G))
                rho[i] = rl if si < ss else rsl
            else:                                     # left fan
                asl = al*(ps/pl)**((G-1)/(2*G))
                if si < ul-al: rho[i] = rl
                elif si > us-asl: rho[i] = rl*(ps/pl)**(1/G)
                else:
                    c = 2/(G+1)+(G-1)/((G+1)*al)*(ul-si)
                    rho[i] = rl*c**(2/(G-1))
        else:
            if ps > pr:
                rsr = rr*((ps/pr+(G-1)/(G+1))/((G-1)/(G+1)*ps/pr+1))
                ss = ur+ar*np.sqrt((G+1)/(2*G)*ps/pr+(G-1)/(2*G))
                rho[i] = rr if si > ss else rsr
            else:
                asr = ar*(ps/pr)**((G-1)/(2*G))
                if si > ur+ar: rho[i] = rr
                elif si < us+asr: rho[i] = rr*(ps/pr)**(1/G)
                else:
                    c = 2/(G+1)-(G-1)/((G+1)*ar)*(ur-si)
                    rho[i] = rr*c**(2/(G-1))
    return rho
 
def to_prim(U):
    r = U[0]; u = U[1]/r; p = (G-1)*(U[2]-0.5*r*u*u)
    return r, u, p
 
def hllc_flux(rL, uL, pL, rR, uR, pR):
    aL, aR = np.sqrt(G*pL/rL), np.sqrt(G*pR/rR)
    EL = pL/(G-1)+0.5*rL*uL**2; ER = pR/(G-1)+0.5*rR*uR**2
    sL = np.minimum(uL-aL, uR-aR); sR = np.maximum(uL+aL, uR+aR)
    sM = (pR-pL+rL*uL*(sL-uL)-rR*uR*(sR-uR))/(rL*(sL-uL)-rR*(sR-uR))
    FL = np.array([rL*uL, rL*uL**2+pL, uL*(EL+pL)])
    FR = np.array([rR*uR, rR*uR**2+pR, uR*(ER+pR)])
    UL = np.array([rL, rL*uL, EL]); UR = np.array([rR, rR*uR, ER])
    UsL = rL*(sL-uL)/(sL-sM)*np.array([np.ones_like(rL), sM, EL/rL+(sM-uL)*(sM-pL/(rL*(sL-uL)))])
    UsR = rR*(sR-uR)/(sR-sM)*np.array([np.ones_like(rR), sM, ER/rR+(sM-uR)*(sM-pR/(rR*(sR-uR)))])
    F = np.where(sL >= 0, FL, np.where(sM >= 0, FL+sL*(UsL-UL),
        np.where(sR >= 0, FR+sR*(UsR-UR), FR)))
    return F
 
def minmod_slope(q, limited):
    dL = q[1:-1]-q[:-2]; dR = q[2:]-q[1:-1]
    if not limited:
        return 0.5*(dL+dR)
    return np.where(dL*dR <= 0, 0.0, np.where(np.abs(dL) < np.abs(dR), dL, dR))
 
def residual(U, h, periodic, limited):
    V = np.vstack(to_prim(U))                      # 3 x (N+4)
    sl = np.zeros_like(V)
    for k in range(3):
        sl[k, 1:-1] = minmod_slope(V[k], limited)
    VL = V[:, 1:-1]+0.5*sl[:, 1:-1]                # right face of cell i
    VR = V[:, 2:]-0.5*sl[:, 2:]                    # left  face of cell i+1
    F = hllc_flux(VL[0], VL[1], VL[2], VR[0], VR[1], VR[2])
    R = np.zeros_like(U)
    R[:, 2:-2] = -(F[:, 1:-1]-F[:, :-2])/h
    return R
 
def fill_ghost(U, periodic):
    if periodic:
        U[:, :2] = U[:, -4:-2]; U[:, -2:] = U[:, 2:4]
    else:
        U[:, 0] = U[:, 2]; U[:, 1] = U[:, 2]
        U[:, -1] = U[:, -3]; U[:, -2] = U[:, -3]
 
def evolve(U, h, tend, periodic, limited, cfl=0.4):
    t = 0.0
    while t < tend-1e-14:
        fill_ghost(U, periodic)
        r, u, p = to_prim(U[:, 2:-2])
        dt = cfl*h/np.max(np.abs(u)+np.sqrt(G*p/r))
        dt = min(dt, tend-t)
        U1 = U+dt*residual(U, h, periodic, limited); fill_ghost(U1, periodic)
        U = 0.5*(U+U1+dt*residual(U1, h, periodic, limited))
        t += dt
    fill_ghost(U, periodic)
    return U
 
def order_table(ns, errs):
    return [None]+[np.log(errs[i-1]/errs[i])/np.log(ns[i]/ns[i-1]) for i in range(1, len(ns))]
 
def run_sod(ns, tend=0.2):
    e1, ei = [], []
    for N in ns:
        h = 1.0/N
        x = (np.arange(N)+0.5)*h
        r = np.where(x < 0.5, 1.0, 0.125); u = np.zeros(N)
        p = np.where(x < 0.5, 1.0, 0.1)
        U = np.zeros((3, N+4))
        U[:, 2:-2] = np.vstack([r, r*u, p/(G-1)+0.5*r*u*u])
        U = evolve(U, h, tend, False, True)
        re = sample_exact_riemann(x-0.5, tend, 1.0, 0.0, 1.0, 0.125, 0.0, 0.1)
        d = np.abs(U[0, 2:-2]-re)
        e1.append(np.sum(d)*h); ei.append(np.max(d))
    return e1, ei
 
def run_wave(ns, limited, tend=1.0):
    e1 = []
    for N in ns:
        h = 1.0/N
        x = (np.arange(N)+0.5)*h
        prof = lambda z: 1.0+0.5*np.exp(-60.0*(np.mod(z-0.5, 1.0)-0.5)**2)
        r = prof(x); u = np.ones(N); p = np.ones(N)
        U = np.zeros((3, N+4))
        U[:, 2:-2] = np.vstack([r, r*u, p/(G-1)+0.5*r*u*u])
        U = evolve(U, h, tend, True, limited)
        e1.append(np.sum(np.abs(U[0, 2:-2]-prof(x-tend)))*h)
    return e1
 
 
def rankine_hugoniot(m1):
    """standing-shock jump ratios for the shock-vortex setup"""
    rr = ((G+1)*m1**2)/((G-1)*m1**2+2)
    pr = (2*G*m1**2-(G-1))/(G+1)
    m2 = np.sqrt((1+0.5*(G-1)*m1**2)/(G*m1**2-0.5*(G-1)))
    return rr, pr, m2
 
NS = [100, 200, 400, 800]
s1, si = run_sod(NS)
wl, wc = run_wave(NS, True), run_wave(NS, False)
o_s1, o_si, o_wl, o_wc = (order_table(NS, e) for e in (s1, si, wl, wc))
fmt = lambda o: "    -  " if o is None else f"{o:7.3f}"
print("Sod shock tube (two discontinuities), t = 0.2")
print("    N      L1 error   order     Linf error  order")
for i, N in enumerate(NS):
    print(f" {N:4d}   {s1[i]:.4e} {fmt(o_s1[i])}   {si[i]:.4e} {fmt(o_si[i])}")
print()
print("smooth density wave (exact = translated Gaussian), t = 1.0")
print("    N      L1 minmod  order     L1 unlim.   order")
for i, N in enumerate(NS):
    print(f" {N:4d}   {wl[i]:.4e} {fmt(o_wl[i])}   {wc[i]:.4e} {fmt(o_wc[i])}")
ps, us = riemann_star_pressure(1.0, 0.0, 1.0, 0.125, 0.0, 0.1)
print(f"\nSod star state   p* = {ps:.6f}  u* = {us:.6f}")
rr, pr, m2 = rankine_hugoniot(1.2)
print(f"RH at M1 = 1.2   rho2/rho1 = {rr:.4f}  p2/p1 = {pr:.4f}  M2 = {m2:.4f}")
amp = (G-1)*25.0/(8*G*np.pi**2)
for r in (2.5, 5.0):
    print(f"vortex tail r = {r:.1f}  exp(1-r^2) = {np.exp(1-r*r):.3e}  dT(beta=5) = {amp*np.exp(1-r*r):.3e}")

Output:

Sod shock tube (two discontinuities), t = 0.2
    N      L1 error   order     Linf error  order
  100   6.5089e-03     -     6.7837e-02     -  
  200   3.7579e-03   0.792   8.3037e-02  -0.292
  400   2.1461e-03   0.808   8.3428e-02  -0.007
  800   1.2420e-03   0.789   8.5667e-02  -0.038
 
smooth density wave (exact = translated Gaussian), t = 1.0
    N      L1 minmod  order     L1 unlim.   order
  100   5.6574e-03     -     2.0251e-03     -  
  200   2.0125e-03   1.491   5.0589e-04   2.001
  400   5.7467e-04   1.808   1.2630e-04   2.002
  800   1.5736e-04   1.869   3.1562e-05   2.001
 
Sod star state   p* = 0.303130  u* = 0.927453
RH at M1 = 1.2   rho2/rho1 = 1.3416  p2/p1 = 1.5133  M2 = 0.8422
vortex tail r = 2.5  exp(1-r^2) = 5.248e-03  dT(beta=5) = 4.747e-04
vortex tail r = 5.0  exp(1-r^2) = 3.775e-11  dT(beta=5) = 3.415e-12

Raise the sharpness slider $m$ and the bump goes from a Gaussian toward a top hat, while the measured order slides from near 2 toward 1.

The unlimited line that printed 2.001 drifts toward 1 as $m$ grows. show Linf takes the same run from 2 down to 0. The scheme line does not change.

The values $p^{*} = 0.303130$ and $u^{*} = 0.927453$ are the textbook Sod star state, which puts the shock speed at 1.7522. Whether the code reproduces those numbers is what Sod's shock tube is really looking at.

The job the isentropic vortex holds#

The standard case for measuring order in two dimensions is isentropic vortex convection. Lay a uniform stream over a periodic $[0,10]^2$ domain and add a vortex perturbation on top.

Here $\beta$ is the vortex strength, $r$ the distance from the vortex center, and $(x_c, y_c)$ the center coordinates. Temperature follows from the isentropic relation.

This field is an exact solution that translates at the free-stream velocity. It is smooth and has no discontinuity, so measuring order on it is legitimate.

The trap is in the tails. The perturbation decays like $e^{1-r^2}$ but never reaches zero. Clip it at the domain boundary and the "exact" solution stops being exact.

On $[0,10]^2$ the distance to the boundary is 5. There $e^{1-r^2}$ is $3.8 \times 10^{-11}$, and with $\beta = 5$ the temperature perturbation is $3.4 \times 10^{-12}$. That is comfortably safe for a second-order scheme. Shrink the domain to $[0,5]^2$ and at $r = 2.5$ the same quantities become $5.2 \times 10^{-3}$ and $4.7 \times 10^{-4}$. A high-order scheme hits that floor while it is still refining, and the order collapses there.

Shock–vortex does not measure order#

The third case has a different character. A standing $M = 1.2$ shock sits on the left, the right state follows from the Rankine–Hugoniot relations, and an isentropic vortex is pushed through the shock.

This problem has no closed-form solution, so no order can be extracted from it.

What can be checked is the jump. At $M_1 = 1.2$ and $\gamma = 1.4$ the ratios are $\rho_2/\rho_1 = 1.3416$, $p_2/p_1 = 1.5133$, $M_2 = 0.8422$. If the code cannot reproduce these, every picture downstream of the shock is meaningless. Mass and momentum integrals along the periodic direction should stay fixed to machine precision. A functional such as the acoustic amplitude or an integral of vorticity gets measured on three grids, with the finest as the reference. That is grid convergence, not an order of accuracy.

What you pick when you write one number in a report#

The sentence "second-order accuracy confirmed" hides three choices: which solution, which norm, and whether the limiter was on.

If Sod reports 0.79, the thing to fix is the sentence in the report, not the scheme. If the vortex reports 2.00, the limiter setting belongs next to it. Measuring order against a smooth exact solution is the same job as checking a pipe flow against the Poiseuille profile. Choosing where to measure is half of the work.