The Shock Never Left the Starting Line — Where Conservative and Primitive Forms Split

I once wrote two versions of a 1D Burgers solver and ran them side by side. One differenced fluxes; the other multiplied a velocity by a gradient. On paper the two forms convert into each other with a single application of the chain rule. Then I fed them a Riemann problem, and one shock refused to move at all. This post traces that stall back to the control-volume derivation and shows, in Python, that refining the grid eight-fold does not fix it.

The shock never left the starting line#

The same equation admits two writings. The conservative form pairs a rate of change with a flux divergence.

The primitive (non-conservative) form expands the derivative into a velocity times a gradient.

Where $u$ is smooth, $\partial_x(u^2/2) = u\,\partial_x u$, so the two are the same equation. The product rule says so.

Now hand them a Riemann problem: $u_L = 1$ on the left, $u_R = 0$ on the right. The exact answer is a shock traveling right at speed $0.5$. The conservative Godunov scheme returned $0.5006$. The primitive upwind difference returned $0.0000$. The shock stayed exactly where it started.

Drop u_R toward 0 and only the pink front stops; raising grid N does not start it moving again.

Green on top is the conservative track, pink below is the primitive one, and the white dashed line marks the exact shock position. Watch how far the pink front lags that line, and whether N = 320 closes the gap at all.

The equation came out of a control volume already in flux form#

Why is the flux form the original? Retrace the derivation and it shows.

Take a small box $dx\,dy\,dz$ and count the mass crossing each face. What passes through one face is the density at the face center, the velocity normal to it, and the area, multiplied together: $\dot m = \rho V_n A$. Face-center values come from a Taylor expansion about the cell center with second order and beyond dropped. Sum all six faces, divide by $dx\,dy\,dz$, and continuity falls out.

Momentum follows the same recipe: momentum carried through the faces, plus body forces, plus surface forces.

Here $\rho$ is density, $u_i$ the velocity components, $p$ the pressure, and $\tau_{ij}$ the viscous stress tensor — the deviatoric stress, linear in the velocity gradients under the Newtonian-fluid assumption.

What matters is not the shape of these equations but their origin. Every term is defined as something that crossed a face. The divergence form is not a stylistic choice; it is what the derivation produced.

The primitive form takes one further step. Expand the product derivative, subtract the continuity equation times $u_i$, and divide through by $\rho$. Assume $\rho = \text{const}$ and you get $\partial u_i / \partial x_i = 0$ along with

Every one of those manipulations assumes differentiability. Across a discontinuity there is none to assume.

One division erased the telescoping sum#

At the discrete level the loss is easier to see. The finite-volume update in conservative form reads

Sum this over all cells. An interior face flux $f_{i+1/2}$ is subtracted in cell $i$ and added in cell $i+1$. The signs are opposite, so it cancels exactly. That is the telescoping sum, and only the fluxes at the two domain ends survive it.

The change in the total equals what crossed the boundaries, with no exceptions beyond round-off.

The primitive form breaks here. The term $u_i \cdot (u_i - u_{i-1})/\Delta x$ carries a different coefficient $u_i$ in front of every cell. Neighboring contributions have different magnitudes and no longer cancel. The leftover accumulates every step.

The code below measures it. In the second Riemann problem ($u_L = 1$, $u_R = 0.4$) the amount that should have entered the domain is $+0.168000$. The conservative scheme reproduces that to six decimals. The primitive one gives $+0.152765$, losing roughly 9%.

The same leak is already visible before a front exists. Where characteristics collide and a Burgers shock is born shows the two forms agreeing until that crossing. Once the ledger of the front no longer balances, the total drifts.

Rankine–Hugoniot answers only to the flux#

Where does the shock speed come from? Apply the conservation law to a thin control volume wrapped around the discontinuity.

Here $s$ is the propagation speed of the discontinuity and $f$ is the flux. For Burgers, $f = u^2/2$, which gives

Only $f$ appears in that relation. The expression $u\,\partial_x u$ never shows up, and it could not. At a discontinuity $\partial_x u$ is a delta function, and multiplying it by a $u$ that jumps is undefined in distribution theory. Such a term is called a non-conservative product.

The Lax–Wendroff theorem guards exactly this boundary. If the numerical solution of a conservative scheme converges, the limit is necessarily a weak solution of the conservation law, and therefore satisfies Rankine–Hugoniot. Non-conservative schemes carry no such guarantee. What Hou and LeFloch showed is worse: they do converge, but to the wrong speed.

Measuring the speed and the total in Python#

Same grid, same CFL, same initial data, two schemes. Standard library only.

def riemann_setup(nx, ul, ur, xs=0.3):
    dx = 1.0 / nx
    return dx, [ul if (i + 0.5) * dx < xs else ur for i in range(nx)]
 
def godunov_flux(a, b):
    if a > b:                                  # shock: pick the upwind side
        return 0.5 * a * a if a + b >= 0 else 0.5 * b * b
    if a >= 0:
        return 0.5 * a * a
    return 0.5 * b * b if b <= 0 else 0.0      # transonic rarefaction
 
def step_conservative(u, dx, dt):              # u_t + (u^2/2)_x = 0
    n = len(u)
    f = [0.5 * u[0] ** 2] + [godunov_flux(u[i], u[i + 1]) for i in range(n - 1)] \
        + [0.5 * u[-1] ** 2]
    return [u[i] - dt / dx * (f[i + 1] - f[i]) for i in range(n)]
 
def step_primitive(u, dx, dt):                 # u_t + u u_x = 0
    n, out = len(u), []
    for i in range(n):
        im, ip = max(i - 1, 0), min(i + 1, n - 1)
        g = (u[i] - u[im]) / dx if u[i] >= 0 else (u[ip] - u[i]) / dx
        out.append(u[i] - dt * u[i] * g)
    return out
 
def shock_locate(u, dx, level):
    for i in range(1, len(u)):
        if u[i] < level <= u[i - 1]:
            return (i - 0.5) * dx + dx * (u[i - 1] - level) / (u[i - 1] - u[i])
    return float("nan")
 
def march_burgers(nx, ul, ur, tend, step):
    dx, u = riemann_setup(nx, ul, ur)
    t = 0.0
    while t < tend - 1e-12:
        dt = min(0.4 * dx / max(max(abs(v) for v in u), 1e-12), tend - t)
        u = step(u, dx, dt)
        t += dt
    return dx, u
 
T, XS = 0.4, 0.3
for ul, ur in ((1.0, 0.0), (1.0, 0.4)):
    s = 0.5 * (ul + ur)
    influx = (0.5 * ul ** 2 - 0.5 * ur ** 2) * T        # exact net flux into the domain
    print("uL=%.1f uR=%.1f | Rankine-Hugoniot speed = %.3f" % (ul, ur, s))
    print("    N   conservative   primitive")
    for nx in (100, 200, 400, 800):
        v = []
        for step in (step_conservative, step_primitive):
            dx, u = march_burgers(nx, ul, ur, T, step)
            v.append((shock_locate(u, dx, s) - XS) / T)
        print("%5d      %7.4f     %7.4f" % (nx, v[0], v[1]))
    for name, step in (("conservative", step_conservative), ("primitive  ", step_primitive)):
        dx, u = march_burgers(400, ul, ur, T, step)
        dx0, u0 = riemann_setup(400, ul, ur)
        print("  N=400 %s : d(int u dx) = %+.6f  (exact %+.6f)"
              % (name, sum(u) * dx - sum(u0) * dx0, influx))
    print()

uL=1.0 uR=0.0 | Rankine-Hugoniot speed = 0.500
    N   conservative   primitive
  100       0.5006      0.0000
  200       0.5003      0.0000
  400       0.5002      0.0000
  800       0.5001      0.0000
  N=400 conservative : d(int u dx) = +0.200000  (exact +0.200000)
  N=400 primitive   : d(int u dx) = +0.000000  (exact +0.200000)
 
uL=1.0 uR=0.4 | Rankine-Hugoniot speed = 0.700
    N   conservative   primitive
  100       0.7009      0.6263
  200       0.7005      0.6330
  400       0.7002      0.6363
  800       0.7001      0.6379
  N=400 conservative : d(int u dx) = +0.168000  (exact +0.168000)
  N=400 primitive   : d(int u dx) = +0.152765  (exact +0.168000)

The first case is the extreme one. With $u_R = 0$, the term $u\,\partial_x u$ vanishes identically in every cell to the right of the jump. There is nothing to update, so the front never starts moving. The total change is exactly zero as well: the $0.2$ that entered through the left boundary shows up nowhere.

The stall and the missing mass are the same arithmetic. With $u_L=1$ and $u_R=0$, $s=0.5$. Over $T=0.4$ the front must move $0.2$, from $x_s=0.3$ to $0.5$. Speed 0.5006 at N=100 is a position error of $0.00024$. Speed 0.0000 leaves the front at $x=0.3$.

Boundary flux gives that same 0.2: $u_L^2/2 = 0.5$ on the left, zero on the right, times $T$. That is the printed $+0.200000$. The primitive line $+0.000000$ drops the 0.2 out of the cell sum.

For $u_R=0.4$, $s=0.7$ and the net influx is $(0.5-0.08)\times 0.4 = 0.168$. The primitive total $+0.152765$ is short by $0.015235$, about 9% of 0.168. At N=800 the speed is still 0.6379, and $(0.700-0.6379)/0.700 = 0.0887$. The primitive steps are 0.0067, then 0.0033, then 0.0016, so the column settles near 0.639 rather than 0.700. On the same four grids the conservative column reads 0.7009, 0.7005, 0.7002, 0.7001. Its gap to 0.700 shrinks from 0.0009 to 0.0001. The primitive gap, about 0.062, does not follow that ratio.

Does refining the grid rescue it?#

The second case is the dangerous one in practice. The primitive speed walks through $0.6263 \to 0.6330 \to 0.6363 \to 0.6379$. Refine eight-fold and the value settles down. It looks like convergence.

The problem is where it converges. The right answer is $0.700$, and this sequence is heading for roughly $0.639$ — about 8.7% low. An honest grid convergence study will not catch that. You confirm the values on three grids are approaching each other, write "converged," and move on.

The conservative scheme runs $0.7009 \to 0.7001$, closing on the exact value with an error proportional to $\Delta x$. The gap between the two sequences is not a difference in accuracy but a difference in which equation is being solved.

None of this shows up when the solution stays smooth. A code validated only on cases like Taylor–Green passes clean. The moment the first discontinuity forms, a code that had been right until then quietly begins solving different physics.

Where the primitive form still belongs — the shelf life of $\rho=\text{const}$#

None of this makes the primitive form unusable. Nearly all incompressible solvers use it, for good reasons.

The unknown count drops. Two-dimensional compressible flow carries $\rho, u, v, p, T$ — five unknowns needing mass, two momentum components, energy, and an equation of state. Incompressible flow freezes $\rho$ and drops the energy equation and the state relation with it. Only $u, v, p$ remain. Pressure stops being a thermodynamic variable and becomes the Lagrange multiplier of the divergence constraint, so it is solved through a pressure Poisson equation.

Incompressible flow has no shock. There is no discontinuity for Rankine–Hugoniot to govern, so the failure above never arises.

The shelf life is set by the Mach number. The isentropic relation gives the density variation as

where $\rho_0$ is the stagnation density, $\gamma$ the ratio of specific heats, and $M$ the Mach number. Expanded for small $M$, the density change goes as $M^2/2$: about 2% at $M = 0.2$ and about 4.5% at $M = 0.3$. The familiar $M < 0.2$ rule of thumb is that number.

Drag exit Mach up from 0.05 and the spacing between the green markers (density allowed to vary) and the pink ones (density frozen) opens up. Below 0.2 the two rows track each other; past 0.3 the yellow dot pulls away from the $M^2/2$ dashed curve.

When the shock arrives late, look here first#

If the solver parks the shock in the wrong place, check whether the time update is a difference of face fluxes. The change in $\sum_i u_i \Delta x$ has to match the boundary flux digit for digit. At N=400 and $u_R=0$ the conservative line is $+0.200000$ and the primitive line is $+0.000000$. Those two cells already separate the forms. Curvilinear or axisymmetric terms that were slid out of the divergence and into the source produce the same symptom. Smooth data stay quiet. The speed goes wrong at the first jump.

A product such as $\alpha\,\partial_x p$ in a multiphase model is non-conservative in principle and needs its own path-integral reading. With a term like that, taking the grid from 100 to 800 still leaves the speed near 0.6379. The ledger that solves a Riemann problem in every cell is the same book as Godunov's shock tube inside one grid cell. If refinement leaves the shock sitting still, doubt the update, not the accuracy.