Where Negative Density Leaks In — A Positivity-Preserving Flux Limiter for Blast Waves

2:27 AM. A Le Blanc shock tube spat out NaN after 0.003 s. CFL at 0.1, grid doubled — both curves still collapsed. The cause was small. At the head of a strong rarefaction, one cell's density crossed below zero. Two steps later the sound speed was NaN. The step after that killed the domain. The θ-blend of Hu, Adams, and Shu (2013) backs the high-order flux off by one notch, and only there.

Drop the right-state density ρ_R below 10⁻³ and the allowed curve enters the red region, while θ on that face falls under 0.2.

The red region is where the next density breaks the floor. If θ stops on that boundary, the face is effectively first-order Lax–Friedrichs.

Where negative density leaks in#

A high-order scheme wobbles between stability and accuracy. First-order upwind smears a shock across about five cells. Second-order central difference plants a sawtooth on both sides. Light ringing only shakes the pressure. Once the density ratio of neighboring cells passes 100:1 — Le Blanc, a Sedov blast, a supersonic nozzle wake — the failure changes character.

Density conservation on cell $i$ is

$F^\rho = \rho u$ is the density flux. In a cell with $F_{i+1/2}^\rho \gg F_{i-1/2}^\rho$, mass leaving faster than it arrives, $\rho_i^{n+1}$ can go negative. In 1D Euler that cell is the first one in the expansion just behind a strong shock.

Negative density kills the sound speed first. With $a = \sqrt{\gamma p / \rho}$, $\rho \to 0^-$ makes $a^2 < 0$. The next flux evaluation returns NaN from $\sqrt{\cdot}$. One cell's NaN reaches both neighbors. Two steps later the domain is gone.

Lax–Friedrichs — blurry but never broken#

The LF flux Lax wrote down in 1954 is outside this trap.

$\alpha$ is the largest signal speed of the two cells, and $U = (\rho, \rho u, E)^T$ is the conservative state. The subtracted term is artificial viscosity. It averages the neighbors and kills the sawtooth. Same move as the viscosity von Neumann deliberately smeared across a shock in 1950. Under CFL$\leq$1, Lax showed $\rho^{n+1} \geq 0$, a stripped-down form of Hu et al., Theorem 1. Fast shock or weak shock, LF stays up.

Accuracy is the cost. A two-cell shock becomes five cells. A contact smears by another cell every step. Nobody ships LF alone. The safety margin is still where the limiter starts.

The θ-blend — back off one notch at a time#

A high-order flux $F_{i+1/2}^{HO}$ (WENO, MUSCL, central) is sharper than LF and can still make negative density. Hu, Adams, and Shu tie the two fluxes in one line.

$\theta \in [0, 1]$ is a weight per face. $\theta=1$ is pure high order. $\theta=0$ is pure LF. The limiter does one thing: pick the largest $\theta$ that leaves every next density positive. It does not retreat across the domain. One notch.

The LF state already has $\rho^{LF}_i \geq \rho_{\min}$, so each face solves for the largest $\theta$ that keeps the update positive after adding $\Delta F = F^{HO} - F^{LF}$. A face touches two cells. Both must clear the floor, so the smaller bound wins.

That is the cap that dropped θ under 0.2 when ρ_R came down. The face is then solved as first-order LF.

Two-step limiting — density first, pressure next#

Positive density is not enough. Negative pressure kills the solver the same way, because the sound speed goes imaginary. Hu et al. split the limiter into two passes.

Density pass. The formula above sets $\theta^\rho$ on each face so every next density stays at or above $\rho_{\min}$.

Pressure pass. With the density flux frozen by $\theta^\rho$, the momentum and energy fluxes get a $\theta^p \leq \theta^\rho$. The pressure constraint

is nonlinear, a quadratic in general. A conservative choice, taking the smaller root for $\theta^p$, fits on one line.

One weight. $\theta = \min(\theta^\rho, \theta^p)$ blends every conservative flux. One face, one $\theta$, so the conserved variables stay consistent with each other.

50 lines of NumPy — auto-θ on a Le Blanc tube#

Le Blanc here is a strong tube, $\rho_L=1$, $p_L=100$ against $\rho_R=0.02$, $p_R=0.05$. A strong rarefaction runs left and the first cell's density goes negative if nothing guards it.

Divide the two states. The density ratio is $1/0.02 = 50$. The pressure ratio is $100/0.05 = 2000$. The 100:1 line above was a rule of thumb, and this pair's density jump sits under it. The violence is in the pressure. With $\gamma = 1.4$, $a = \sqrt{\gamma p / \rho}$ gives $\sqrt{140} \approx 11.83$ on the left and $\sqrt{3.5} \approx 1.871$ on the right. The ratio is exactly $\sqrt{40} \approx 6.32$.

The grid numbers say the timestep was not the killer. $N = 200$ so $\Delta x = 1/200 = 0.005$, and $\Delta t = 9.0 \times 10^{-5}$. Then $\Delta t/\Delta x = 0.018$. An acoustic Courant number on the left state alone is $a_L \Delta t/\Delta x \approx 0.213$. That matches a run that already died at CFL 0.1. From $\rho^{n+1} = \rho^n - (\Delta t/\Delta x)\Delta F$, once the flux difference exceeds the density slack, a smaller step still punches through the floor.

At the right state $\rho = 0.02$, the slack down to the floor $10^{-3}$ is $0.019$. For $\theta$ to stay at 1, the density correction $|d|$ cannot exceed $0.019/0.018 \approx 1.06$. At $|d| = 2$ the cap is $0.019/(0.018 \times 2) \approx 0.53$. If the slack shrinks to $5 \times 10^{-4}$, the same $|d| = 1$ drops the cap to $0.0005/0.018 \approx 0.028$. That face is almost pure LF.

By weight the blend is $F = (1-\theta)F^{LF} + \theta F^{HO}$. At $\theta = 0.53$ the high-order share is 53%. At $\theta = 0.028$ it is 2.8% and the rest is LF. The cap is the update bound itself: the correction moves density by $(\Delta t/\Delta x)\,\theta\,|d|$, which must stay inside the slack $0.019$, so $\theta \le 0.019/(0.018\,|d|)$. Equality at $\theta = 1$ sits near $|d| \approx 1.06$.

The pressure formula reads the same way. $p = (\gamma - 1)(E - \tfrac{1}{2}\rho u^2)$ and $\gamma - 1 = 0.4$. Hold $E$ and $\rho$ fixed and one unit of kinetic energy cuts pressure by 0.4. The right state $p = 0.05$ hits zero at a kinetic energy of $0.05/0.4 = 0.125$. The left state $p = 100$ survives up to $250$. Guard only the density floor and the leak moves to a positive density with an imaginary sound speed.

import numpy as np
 
GAMMA = 1.4
RHO_FLOOR = 1e-3
N, NSTEP = 200, 220
dx = 1.0 / N
dt = 9.0e-5
 
def cons_to_prim(U):
    rho, mu, E = U[0], U[1], U[2]   # conservative -> primitive
    u = mu / rho
    p = (GAMMA - 1.0) * (E - 0.5 * rho * u * u)
    return u, p
 
def euler_flux(U):
    u, p = cons_to_prim(U)
    return np.array([U[1], U[1] * u + p, u * (U[2] + p)])
 
def alpha_face(UL, UR):
    uL, pL = cons_to_prim(UL); uR, pR = cons_to_prim(UR)
    aL = np.sqrt(GAMMA * max(pL, 0.0) / UL[0])
    aR = np.sqrt(GAMMA * max(pR, 0.0) / UR[0])
    return max(abs(uL) + aL, abs(uR) + aR)
 
def density_theta(U, dF, dt_dx):
    """Per-face θ: largest weight that keeps ρ above the floor."""
    theta = np.ones(N + 1)
    for f in range(1, N):
        d = dF[0, f]                          # density anti-diffusion
        if abs(d) < 1e-14: continue
        rL_room = U[0, f - 1] - RHO_FLOOR     # left-cell headroom
        rR_room = U[0, f]     - RHO_FLOOR     # right-cell headroom
        cap = (rL_room if d > 0 else rR_room) / (dt_dx * abs(d))
        theta[f] = max(0.0, min(theta[f], cap))
    return theta

The driver is one first-order LF pass, one central anti-diffusion pass, one pass for θ. 220 steps, 200 cells. Limiter off: NaN near step 30. Limiter on: the run finishes.

Switch the limiter off and the red curve drags density below zero near step 30.

That crossing is the collapse. Cyan flattens near $\rho \approx 0.02$ at the same time. The accuracy gap is under 5%. One run lives. The other dies.

What the limiter costs#

A face with $\theta = 0$ is solved by first-order LF. Shock fronts and rarefaction heads go blunt there. That region is under 5% of the domain, so the global error barely moves. Mean θ sits between 0.95 and 1.0.

Set $\rho_{\min}$ to something like $10^{-12}$ and the limiter wakes up late. The same NaN comes back. In practice the floor is near $10^{-4}$ of the mean density. Multi-stage Runge–Kutta has to rerun the limiter on every stage. Pass only the first stage and the next stage diverges anyway.

A floor of 1e-12 wakes the limiter after the NaN#

$\rho_{\min} = 10^{-12}$ is decoration, not a floor. Against a right-state density of 0.02 almost the whole slack is still unused, and the cap above stays near 1. The bound tightens only after density has already fallen below $10^{-8}$. By then $\gamma p / \rho$ has changed sign. The limiter retreats to LF after the cell is dead.

Face θ guards the cell average. If that average stays positive while a reconstructed point inside the cell crosses zero, this one line never sees the point. Whether a reconstruction that stays positive in the mean and goes negative at a point still needs its own fix is the next cut.

“

Reference: X. Y. Hu, N. A. Adams, C.-W. Shu. Positivity-preserving method for high-order conservative schemes solving compressible Euler equations. J. Comput. Phys. 242 (2013) 169–180.