The Fake Extrema a Second-Order Reconstruction Invents on Unstructured Grids — Barth–Jespersen and Venkatakrishnan Limiters

The solver had been fine. Second order went on, and the density went negative. Face values ahead of the shock sat above the cell average. First-order upwind on the same mesh did not. Raising the order is what killed the run. A straight line inside an unstructured cell invents extrema no neighbor holds. Barth–Jespersen cuts those off. Venkatakrishnan rounds the cut.

Reconstruction invents extrema that were never there#

In the finite volume method (FVM, a discretization that carries cell averages as unknowns), second order is a straight line inside the cell.

$u_i$ is the average in cell $i$. $\nabla u_i$ is the reconstructed gradient, $\mathbf{x}_f$ is the face center, and $\Phi_i\in[0,1]$ is the limiter.

With $\Phi_i=1$ the limiter is off. Least squares or Green–Gauss builds the gradient by averaging neighbors. Next to a jump that average is too steep. Run the line out to the face and a new max or min appears, one no neighbor stored.

Raise the gradient gain and the unlimited reconstruction (Φ=1) punches through the gray band, the neighbor max and min, and turns red.

While the line stays inside the band, that cell's Φ is 1. Only the cells beside the jump shorten first. Switch the limiter to Barth–Jespersen at the same gain and the line stops at the band edge.

The fake extremum poisons the next flux. The oscillation ends at negative density and negative energy. A reconstructed face value must not leave the range its neighbors already hold. That is the local maximum principle (LMP). $\Phi_i$ is the scalar that shaves the gradient to obey it.

Barth–Jespersen: caged by the neighbor max and min#

Barth and Jespersen (1989) fix the allowed range first. It is the max and min of cell $i$ and its neighbors.

$N(i)$ is the face-neighbor set of cell $i$. At each face, read the unlimited increment $u_f-u_i$ and form a per-face factor $\phi_f$.

The cell limiter is the most conservative value over the faces.

Inside the range, $\phi_f=1$ and second order stays. A face that would overshoot keeps only enough slope to touch the bound. Near an extremum $\Phi_i\to 0$, and the scheme falls back to first-order upwind.

For a positive increment the factor is $\phi_f=\min(1,(u_i^{\max}-u_i)/(u_f-u_i))$. The moment the numerator is the smaller one, the face value sticks to $u_i^{\max}$ and stops. That ratio is why the red line in the figure above breaks at the gray edge.

The second trap: non-differentiability#

Barth–Jespersen holds monotonicity. Drop it into a steady solver and the residual stalls, buzzing at some level. The $\min$ and the division are why. $\phi_f$ is not differentiable in the solution. Which face supplies the minimum flips from iteration to iteration. The limiter switches on and off, and the residual bounces. An implicit Jacobian hates that kink.

In smooth regions Φ should just be 1. It should engage only at a real discontinuity, and the on/off boundary should be smooth. The next expression rounds the sharp corner.

Venkatakrishnan: shaving it smoothly#

Venkatakrishnan (1993) replaced the kink in $\min(1,y)$ with a rational function. The Michalak–Ollivier-Gooch limiter mentioned in the source sits in the same family. The per-face factor reads

$\Delta_-=u_f-u_i$ is the unlimited increment. $\Delta_+$ is $u_i^{\max}-u_i$ or $u_i^{\min}-u_i$, following the sign. The quantity you actually choose is $\epsilon^2$.

$K$ is the tuning parameter and $\Delta x$ is the mesh size. When the variation is smaller than $\epsilon$, $\phi_f\to 1$ and the limiter switches off. Raise $K$ and the threshold rises, so the limiter lets go over a wider region. Too large, and the wiggles beside a shock slip through. Send $K$ to 0 and the function returns to Barth–Jespersen.

The code treats the cell size as 1 and sets $\epsilon^2=K^3$. At the default $K=0.3$, $\epsilon^2=0.027$. With $\epsilon=0$ and $\Delta_+=\Delta_-=\Delta$, right where the reconstruction lands exactly on the allowed bound, the rational function equals $(3\Delta^2)/(4\Delta^2)=3/4$. Barth–Jespersen would return $\phi_f=1$ at that same point. Venkatakrishnan drops a quarter of the slope before the face ever touches the bound. That is why the corner goes blunt one step early.

The range written down as common practice runs from $K\approx 0.1$ to $K\approx 5$. Plug those ends into the code's $\epsilon^2=K^3$ and you get 0.001 and 125. The threshold ratio is 125,000. A small move in $K$ is cubed, so the set of cells excused as smooth grows fast. Drop $K$ below 0.1 and $\epsilon^2$ is effectively 0, which brings Barth–Jespersen's kink back.

In the figure above, turn Venkatakrishnan on and move $K$. Φ slides, with no stair steps.

Python — racing three limiters on scalar advection#

Periodic boundaries. A square wave and a Gaussian bump advect together. Unlimited (Fromm), Barth–Jespersen, and Venkatakrishnan run side by side. Watch the final minimum. Below zero means a fake trough of that depth.

The cell slope is $s=(u_{i+1}-u_{i-1})/2$. The face sees the increment $0.5s$, so a neighbor difference of 1 gives $s=1/2$ and a face increment of 0.25. If that 0.25 exceeds the allowed headroom, Φ drops below 1. A one-dimensional central difference only sees the left and right neighbors. On an unstructured cell the weights that build the same gradient are a different choice, the one moving least squares refits cell by cell. The limiter attaches after that gradient exists.

import numpy as np
 
NX, A, CFL = 200, 1.0, 0.4
dx = 1.0 / NX
dt = CFL * dx / A
 
def init_profile():
    x = (np.arange(NX) + 0.5) / NX
    u = np.where((x > 0.1) & (x < 0.3), 1.0, 0.0)   # square wave
    u += np.exp(-((x - 0.65) / 0.06) ** 2)          # gaussian bump
    return u
 
def cell_slope(u):
    return (np.roll(u, -1) - np.roll(u, 1)) / 2.0    # central slope (Fromm)
 
def barth_jespersen_phi(u, s):
    up, um = np.roll(u, -1), np.roll(u, 1)
    umax = np.maximum(u, np.maximum(up, um))
    umin = np.minimum(u, np.minimum(up, um))
    phi = np.ones_like(u)
    for du in (0.5 * s, -0.5 * s):                   # two faces
        f = np.ones_like(u)
        pos, neg = du > 1e-12, du < -1e-12
        f[pos] = np.minimum(1.0, (umax[pos] - u[pos]) / du[pos])
        f[neg] = np.minimum(1.0, (umin[neg] - u[neg]) / du[neg])
        phi = np.minimum(phi, f)
    return np.clip(phi, 0.0, 1.0)
 
def venkatakrishnan_phi(u, s, K=0.3):
    up, um = np.roll(u, -1), np.roll(u, 1)
    umax = np.maximum(u, np.maximum(up, um))
    umin = np.minimum(u, np.minimum(up, um))
    eps2 = K ** 3                                    # (K*h)^3, cell size h=1 units
    phi = np.ones_like(u)
    for du in (0.5 * s, -0.5 * s):
        d = np.where(du > 0, umax - u, umin - u)
        num = (d * d + eps2) * du + 2 * du * du * d
        den = d * d + 2 * du * du + d * du + eps2
        f = np.where(np.abs(du) < 1e-12, 1.0, num / (du * den))
        phi = np.minimum(phi, f)
    return np.clip(phi, 0.0, 1.0)
 
def muscl_rhs(u, limiter):
    s = cell_slope(u)
    phi = limiter(u, s) if limiter else np.ones_like(u)
    uL = u + 0.5 * phi * s              # left state at face i+1/2
    flux = A * uL                       # upwind flux (a > 0)
    return -(flux - np.roll(flux, 1)) / dx
 
def advance_muscl(u, limiter):         # SSP-RK2
    k1 = muscl_rhs(u, limiter)
    u1 = u + dt * k1
    k2 = muscl_rhs(u1, limiter)
    return 0.5 * (u + u1 + dt * k2)
 
def run_limiter_race(steps=160):
    fields = {"none": init_profile(), "bj": init_profile(), "venk": init_profile()}
    lims = {"none": None, "bj": barth_jespersen_phi, "venk": venkatakrishnan_phi}
    for _ in range(steps):
        for k in fields:
            fields[k] = advance_muscl(fields[k], lims[k])
    for k, u in fields.items():
        print(f"{k:5s}  min={u.min():+.4f}  max={u.max():.4f}")
 
run_limiter_race()

The run prints this.

none   min=-0.0417  max=1.0231
bj     min=+0.0000  max=1.0000
venk   min=-0.0004  max=1.0009

The unlimited minimum, -0.0417, is how far a square wave of height 1 was dug below zero. The maximum, 1.0231, is the overshoot above that same height. Barth–Jespersen pins min and max to $[0,1]$. Venkatakrishnan's minimum, -0.0004, is about 1/104 of that trough. The amount above 1 shrinks from 0.0231 to 0.0009.

Distance traveled in 160 steps is $A \cdot dt \cdot 160$. With $dx=1/200$ and CFL 0.4, $dt=0.002$, so the distance is 0.32. That is 64 cells out of 200. The square-wave interval $(0.1, 0.3)$ moves to $(0.42, 0.62)$, and the Gaussian center 0.65 arrives at 0.97. Both structures move at the same speed, so their relative placement holds. The printed min and max are not a collision. Each edge produced its own number. The square wave's length, 0.2, is 40 cells. The Gaussian scale 0.06 is 12 cells. A few cells near those edges set the min and the max.

The trailing edge of the square wave is enough. Unlimited (red) grows ripples below zero. Both limiters stay monotone. Venkatakrishnan (amber) is slightly blunter at the corner than Barth–Jespersen (cyan). That points the same way as the algebra, which kept only 3/4 of the slope at $\epsilon=0$. The code's $K=0.3$ does not remove the kink. It only cuts the overshoot down to about 0.0009.

Turning the limiter on in practice#

How you gather neighbors sets the allowed range. Taking the $\max/\min$ in Barth–Jespersen over face neighbors or over vertex neighbors changes the result. On unstructured meshes, vertex neighbors — every cell that shares a vertex — cut the directional bias and balance the limiter.

$K$ comes next. Too small, and convergence stalls the way Barth–Jespersen does. Too large, and shock oscillations get through. Because $\epsilon^2$ carries $\Delta x^3$, refining the mesh makes the limiter fire more often at the same $K$.

In a vector system (the Euler equations), limiting each conserved variable on its own can create new oscillations from cross-component mismatch. Project onto characteristic variables first, then limit. Safer. It costs more.

When the residual stalls, leave the gradient alone#

If a second-order solver dies at negative density, check whether the face value left the neighbor range. If it did, the neighbor set that built Φ is the problem, not the gradient routine. Face neighbors alone miss extrema along the diagonal. The reconstruction leaks through that gap again.

If the residual only marches in place on a steady run, open $K$. $K=0$ is Barth–Jespersen's kink, so the Jacobian jumps every iteration. Where a shock drives Φ to 0, that cell's order sticks at 1. Does the smooth region still sit near second order? If it does, $\epsilon^2$ actually let that region go. What numbers show up when one run measures both at once is still open in the record that puts a shock and a smooth vortex on the same order test.