1950 — When von Neumann Deliberately Smeared the Shock: Artificial Viscosity and Odd-Even Decoupling

Los Alamos, 1950. John von Neumann and Robert Richtmyer put a short paper in the Journal of Applied Physics: "A Method for the Numerical Calculation of Hydrodynamic Shocks". If the shock will not sit in one cell, smear it on purpose. That was the whole line.

The tube below is Sod's, diaphragm gone at $t = 0$: $\rho = 1, p = 1$ on the left, $\rho = 0.125, p = 0.1$ on the right. Raise the ξ slider from 0 to 5 and the post-shock sawtooth dies while the jump spreads over a few cells.

At ξ = 0 the backside of the shock is a saw, near ξ = 2 the wiggles die and the jump is two or three cells wide, and at ξ = 5 the plateau itself goes mushy. Von Neumann's recommended $\xi \approx 2$-$3$ is visibly the sweet spot of the trade-off.

Why the fake viscosity works, how thick it leaves the shock, and the odd-even mode it does not cure: that is the note. The NumPy listing at the end catches one Mach 2 tube.

Shocks Oscillate on a Grid#

Ideal Euler equations carry no viscosity, so a real shock compresses to a few mean free paths. A CFD cell is roughly a million times wider than that. Density and pressure are supposed to double inside one cell.

Pass that discontinuity through almost any conservative scheme and a single phenomenon shows up: a sawtooth oscillation behind the shock, the wiggles. Von Neumann had already met them in 1944 on the Manhattan Project. When the post-shock region rang, the TNT yield calculation came out wrong.

The cause is plain. First-order upwind has strong numerical diffusion that traps shocks but smears everything else. Second-order schemes are accurate but ring at discontinuities, a numerical Gibbs phenomenon. Nature manufactures the missing entropy with molecular viscosity inside the shock. The grid has no such channel. Solving the jump inside one cell as a Riemann problem is Godunov's shock tube in every cell. Artificial viscosity does the opposite: it stretches the jump over about three cells so the wiggles die.

The Single Line from 1950#

The von Neumann-Richtmyer answer was to imitate nature. We can't resolve molecular viscosity, so add a fake bulk pressure on the cell scale.

$Q_i$ is an artificial bulk pressure at cell $i$. $\xi$ is a tuning parameter that sets how many cells the shock will spread over (typically 2-3). $u_{i+1} - u_{i-1}$ is the velocity divergence stencil. $\rho_i$ scales with local density. Replace every $P$ in the momentum and energy equations with $P + Q$. That's it.

$Q$ scales with $(\Delta u)^2$, so it is essentially zero in smooth flow and large only near shocks. The switch is $u_{i+1} \le u_{i-1}$: on only where the velocity is locally converging. Expansion fans and acoustic waves stay untouched. Multiplying by $\rho_i$ gives heavier media a stronger pressure boost.

The written ξ and the code's 0.25 are not the same thickness#

The short form is $Q = \xi^2 (\Delta u)^2 \rho$. It turns on only for $\Delta u < 0$. The listing multiplies by 0.25 * xi**2 * du**2 * rho. That du skips a cell, $u_{i+1} - u_{i-1}$. A one-cell slope puts $2\Delta x$ in the denominator, and the square contributes $1/4$. The 0.25 is that centered difference. The displayed line keeps $\xi^2$ and drops the factor.

At $\xi = 2.5$, $\xi^2 = 6.25$. The code actually multiplies by $0.25 \times 6.25 = 1.5625$. Drop the 0.25 and the same $\xi$ makes $Q$ four times larger. Because $Q$ scales with $\xi^2$, a factor of four is the same as doubling $\xi$. A shock that goes quiet at $\xi = 2$ on the slider matches an implementation of the bare formula at $\xi = 4$. At $\xi = 5$ the plateau is already mushy. Feeding $\xi = 2.5$ into the formula without the 0.25 lands on that blunt side.

The listing fixes 200 cells, a time of 0.2, and a CFL of 0.45. The sentence under the code only says the density peak stays within a few percent of the Mach 2 Rankine–Hugoniot value 2.667. Check the formula with the 0.25 removed against that sentence and the shock thickness being measured is no longer the same one.

Following the Algorithm in NumPy#

Fifty lines hold the core. Donor-cell advection plus the artificial viscosity term, nothing more.

import numpy as np
 
def vnr_shock_solver(N=200, T=0.2, xi=2.5, gamma=1.4):
    dx = 1.0 / N
    rho = np.where(np.arange(N) < N // 2, 1.0, 0.125)
    mom = np.zeros(N)
    p   = np.where(np.arange(N) < N // 2, 1.0, 0.1)
    E   = p / (gamma - 1)
    t = 0.0
    while t < T:
        u = mom / rho
        c = np.sqrt(gamma * p / np.maximum(rho, 1e-9))
        dt = 0.45 * dx / np.max(np.abs(u) + c)
        # von Neumann-Richtmyer Q
        du = np.zeros(N)
        du[1:-1] = u[2:] - u[:-2]
        Q = np.where(du < 0, 0.25 * xi**2 * du**2 * rho, 0.0)
        # donor-cell advection of rho, mom, E
        uf = 0.5 * (u[:-1] + u[1:])
        for q_arr, q_new in [(rho, 'rho'), (mom, 'mom'), (E, 'E')]:
            flux = np.where(uf > 0, q_arr[:-1] * uf, q_arr[1:] * uf)
            q_arr[1:-1] -= dt / dx * (flux[1:] - flux[:-1])
        # pressure update + body force from P+Q
        p = (gamma - 1) * (E - 0.5 * rho * (mom / rho) ** 2)
        Ptot = p + Q
        mom[1:-1] -= dt / (2 * dx) * (Ptot[2:] - Ptot[:-2])
        E[1:-1]   -= dt / (2 * dx) * (Ptot[2:] * u[2:] - Ptot[:-2] * u[:-2])
        t += dt
    return rho, mom, E
 
rho, mom, E = vnr_shock_solver(xi=2.5)
print(f"shock peak rho ≈ {rho.max():.3f}   (Rankine-Hugoniot exact ≈ 2.667 for Mach 2)")

Run it and the peak density lands within a percent of the Rankine-Hugoniot value. Run the same code with ξ = 0 and the oscillation kicks max and min asymmetrically away from the true plateau.

Odd-Even Decoupling: The Other Trap#

Catching the shock isn't the end of the story. The same collocated grid (all variables at cell centers) hides another pathology von Neumann noticed. Take a fluid at rest with constant thermal energy 1 and a sawtooth density.

The pressure $p_i = \rho_i$ in the isothermal case, so cell $i = 2$ sits at twice its neighbors' pressure. It should relax instantly. But the momentum update at cell $i$ reads

because $p_{i+1}$ and $p_{i-1}$ come from the same parity sub-grid. The grid has split into two independent sub-grids that don't see each other. The sawtooth mode neither grows nor decays.

Switch the mode above from collocated to staggered. A staggered grid puts momentum on the cell faces. The pressure difference at face $i+1/2$ is $p_{i+1} - p_i$ — even and odd parities are forced to meet across a single cell. The sawtooth lands inside the first-derivative stencil and is damped at once.

What Was Bought and What Was Lost#

The 1950 trick showed two things at once. Grid limitations could be patched by imitation. And every imitation creates its own grid pathology. Both lessons are alive in modern CFD.

Astrophysics codes like ZEUS, Athena, and PLUTO still use staggered grids with artificial viscosity. OpenFOAM and ANSYS Fluent take collocated grids and bolt Rhie-Chow interpolation on top to dodge the same decoupling. Two paths around the same problem.

A sawtooth behind a shock is not always a bug. It has been a known grid pathology since 1944.

Dropping the 0.25 at ξ = 2.5 lands on the ξ = 5 side of the slider#

Strip the 0.25 and insert $\xi = 2.5$ into the written line, and $Q$ matches what the listing produces at $\xi = 5$. On the slider that setting is already a blunt shock. The range $\xi \approx 2$-$3$ belongs to the code that still carries the 0.25. Move the number across without the factor and the thickness story changes.

Even and odd cells that cannot see each other are a separate illness from the artificial viscosity. Putting velocity on the face lets the pressure jump look at a one-cell neighbor. That layout failing for a reason other than the post-shock sawtooth is recorded in the curl constraint breaking on a staggered grid.