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

Los Alamos, 1950. John von Neumann and Robert Richtmyer published a short paper in the Journal of Applied Physics titled "A Method for the Numerical Calculation of Hydrodynamic Shocks". Its central trick fit in one line — if a shock won't sit cleanly on the grid, then deliberately smear it ourselves before the grid does worse. This post walks through why that single line of artificial viscosity works, what it costs, and what new problem — odd-even decoupling — it left behind on collocated grids. Fifty lines of NumPy will catch a Mach-2 shock by the end.

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.

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.

Three clever choices hide in that line. $Q$ scales with $(\Delta u)^2$, so it is essentially zero in smooth flow and large only near shocks. The condition $u_{i+1} \le u_{i-1}$ — turn on only where the velocity is locally converging — leaves expansion fans and acoustic waves untouched. Multiplying by $\rho_i$ gives heavier media a stronger pressure boost.

Try It#

Below is a classical Sod shock tube: $\rho = 1, p = 1$ on the left, $\rho = 0.125, p = 0.1$ on the right, diaphragm pulled at $t = 0$. The ξ slider dials the artificial viscosity from 0 to 5.

At ξ = 0 a small sawtooth lives behind the shock front. Near ξ = 2 the wiggles die but the shock spreads to two or three cells. Push ξ = 5 and the shock itself gets so blunt that the plateau goes mushy. Von Neumann's recommended $\xi \approx 2$-$3$ is visibly the sweet spot of the trade-off.

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. Consider 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.

One thing worth carrying home: a sawtooth behind a shock is not always a bug. It has been a known grid pathology since 1944.

Things to Remember#

1. Artificial viscosity $Q = \xi^2 (\Delta u)^2 \rho$ activates only under compression ($\Delta u < 0$). Smooth flow stays untouched, shocks spread over ~3 cells, wiggles die.
2. The trade-off is unambiguous. Small ξ leaves oscillations, large ξ blunts the shock. $\xi \approx 2$-$3$ has been the recommended range since the 1950s.
3. Odd-even decoupling is a structural defect of collocated grids. Pick a way around it — staggered (ZEUS-style) or Rhie-Chow (OpenFOAM-style). Opening the source to see which one you are looking at is a habit that pays off.