WENO schemes — high-order shock capturing for hyperbolic equations

1. Hyperbolic conservation laws and the Gibbs problem

A scalar hyperbolic conservation law has the form ∂u/∂t + ∂f(u)/∂x = 0. Even when the flux f(u) and the initial data are smooth, characteristics can cross in finite time and produce a genuine discontinuity — a shock. Systems such as the compressible Euler equations behave the same way: a smooth initial pulse steepens into a shock front carrying a jump in density, velocity and pressure.

Standard high-order finite-difference or finite-volume reconstructions (centred, or fixed-stencil polynomial interpolation) are designed assuming local smoothness. Near a discontinuity that assumption fails, and the reconstructed polynomial overshoots and undershoots the true solution — the numerical analogue of the Gibbs phenomenon in Fourier series. Low order but robust schemes (first-order upwind) avoid oscillation by smearing the shock over many cells, destroying resolution. WENO schemes target the middle ground: high order where the solution is smooth, first-order-like robustness where it is not — without a human-tuned shock detector.

2. From ENO to WENO — the reconstruction idea

The predecessor of WENO is ENO (Essentially Non-Oscillatory, Harten & Osher, 1987). For each cell, ENO builds several candidate polynomial stencils of the same order and picks the single smoothest one — measured by a divided-difference smoothness test — discarding the others entirely. This avoids oscillation but the reconstruction is not smooth in coefficients: a tiny change in the data can flip which stencil is "selected", hurting accuracy and steady-state convergence.

WENO (Liu, Osher & Chan, 1994; Jiang & Shu, 1996) replaces the binary choice with a convex combination of all candidate stencils, each weighted by how smooth it is. In smooth regions, all stencils contribute close to their "ideal" (linear) weights and the combination recovers the maximum possible order (2r−1 for r candidate stencils of r points each). Near a discontinuity, the weight of any stencil that crosses the jump collapses toward zero, and the scheme degrades gracefully toward the smooth, low-order stencils that avoid the jump.

For WENO5 (fifth-order, the most widely used variant) the 5-point stencil above is split into three overlapping 3-point candidate stencils: S₀ = {i−2, i−1, i}, S₁ = {i−1, i, i+1}, S₂ = {i, i+1, i+2}. Each gives a third-order polynomial reconstruction of the value at the cell interface i+½; WENO blends the three.

3. Smoothness indicators

The Jiang-Shu smoothness indicator β_k for stencil S_k measures the total variation of the derivatives of the local reconstruction polynomial p_k over the cell:

β_k is large where the candidate stencil straddles a jump or a steep gradient, and small (ideally zero for a linear function) where the data is smooth. β_k is the single quantity that lets WENO tell the difference between "smooth" and "discontinuous" without an explicit shock-detection flag.

4. Nonlinear weights — WENO5

The three candidate reconstructions of u at the interface i+½ are:

Combined with linear ("ideal") weights d₀ = 1/10, d₁ = 6/10, d₂ = 3/10 they give the optimal fifth-order centred scheme in smooth regions. WENO replaces d_k with nonlinear weights ω_k that depend on β_k:

Where the solution is smooth, β₀ ≈ β₁ ≈ β₂ and ω_k → d_k, recovering the full fifth-order accuracy. Near a shock, the β_k of any stencil crossing the discontinuity blows up, its ω_k → 0, and the reconstruction is dominated by the remaining smooth stencils — automatically and continuously, with no tunable threshold.

5. Flux reconstruction and Riemann solvers

For a system of conservation laws (e.g. the 1D Euler equations), WENO reconstructs not the raw field but the flux itself, after splitting it into left-going and right-going parts (flux splitting, e.g. Lax-Friedrichs: f± = ½(f(u) ± α u) with α = max|characteristic speed|). WENO reconstructs f⁺ using an upwind-biased (leftward) stencil and f⁻ using a mirrored stencil, and the numerical flux at the interface is F_{i+1/2} = f⁺_{i+1/2} + f⁻_{i+1/2}.

For genuinely nonlinear systems, reconstructing in characteristic variables (projecting onto local eigenvectors of the flux Jacobian) rather than conservative variables greatly reduces oscillations near strong shocks and contact discontinuities, at the cost of an extra matrix diagonalisation per cell interface.

6. Time integration — SSP-RK3

A high-order spatial reconstruction is wasted if paired with a low-order or unstable time integrator. WENO is almost always combined with a Strong Stability Preserving Runge-Kutta scheme (Shu-Osher SSP-RK3), which is a convex combination of forward-Euler substeps and therefore inherits the TVD (total variation diminishing) property of the spatial operator:

The time step is limited by the usual CFL condition Δt ≤ CFL · Δx / max|λ|, with CFL typically 0.3–0.6 for WENO5 combined with SSP-RK3.

7. JavaScript implementation — 1D advection and Euler

// ── WENO5 interface reconstruction (left-biased) ─────────────────
function weno5(um2, um1, u0, up1, up2) {
  const EPS = 1e-6;

  // Candidate stencil reconstructions
  const u0r = (1/3)*um2 - (7/6)*um1 + (11/6)*u0;
  const u1r = -(1/6)*um1 + (5/6)*u0  + (1/3)*up1;
  const u2r = (1/3)*u0  + (5/6)*up1 - (1/6)*up2;

  // Smoothness indicators (Jiang-Shu)
  const b0 = (13/12)*(um2-2*um1+u0)**2 + 0.25*(um2-4*um1+3*u0)**2;
  const b1 = (13/12)*(um1-2*u0+up1)**2  + 0.25*(um1-up1)**2;
  const b2 = (13/12)*(u0-2*up1+up2)**2  + 0.25*(3*u0-4*up1+up2)**2;

  // Nonlinear weights (linear weights d = 0.1, 0.6, 0.3)
  const a0 = 0.1 / (EPS+b0)**2;
  const a1 = 0.6 / (EPS+b1)**2;
  const a2 = 0.3 / (EPS+b2)**2;
  const sum = a0 + a1 + a2;

  return (a0*u0r + a1*u1r + a2*u2r) / sum;
}

// ── Lax-Friedrichs flux splitting + WENO5 + SSP-RK3 ──────────────
function wenoRHS(u, flux, alpha, dx, N) {
  const fp = new Float64Array(N), fm = new Float64Array(N);
  for (let i = 0; i < N; i++) {
    const f = flux(u[i]);
    fp[i] = 0.5 * (f + alpha * u[i]);  // right-going part
    fm[i] = 0.5 * (f - alpha * u[i]);  // left-going part
  }
  const rhs = new Float64Array(N);
  for (let i = 2; i < N - 2; i++) {
    const Fp = weno5(fp[i-2], fp[i-1], fp[i], fp[i+1], fp[i+2]);
    const Fm = weno5(fm[i+2], fm[i+1], fm[i], fm[i-1], fm[i-2]);
    const FpL = weno5(fp[i-3]??fp[i-2], fp[i-2], fp[i-1], fp[i], fp[i+1]);
    const FmL = weno5(fm[i+1], fm[i], fm[i-1], fm[i-2], fm[i-3]??fm[i-2]);
    const F_ip = Fp + Fm;      // numerical flux at i+1/2
    const F_im = FpL + FmL;   // numerical flux at i-1/2
    rhs[i] = -(F_ip - F_im) / dx;
  }
  return rhs;
}

// ── SSP-RK3 time step ─────────────────────────────────────────────
function sspRK3Step(u, dt, rhsFn) {
  const N = u.length;
  const u1 = new Float64Array(N), u2 = new Float64Array(N), out = new Float64Array(N);
  const L0 = rhsFn(u);
  for (let i = 0; i < N; i++) u1[i] = u[i] + dt * L0[i];
  const L1 = rhsFn(u1);
  for (let i = 0; i < N; i++) u2[i] = 0.75*u[i] + 0.25*u1[i] + 0.25*dt*L1[i];
  const L2 = rhsFn(u2);
  for (let i = 0; i < N; i++) out[i] = (1/3)*u[i] + (2/3)*u2[i] + (2/3)*dt*L2[i];
  return out;
}

8. Applications and variants

• Sod shock tube: the canonical 1D Euler benchmark — a diaphragm separating two constant states bursts, producing a shock, a contact discontinuity and a rarefaction fan. WENO5 resolves all three with sharp, non-oscillatory profiles.
• WENO-Z and mapped WENO (WENO-M): variants that reduce excess numerical dissipation at smooth extrema while keeping the shock-capturing robustness of classical WENO5.
• Compact / hybrid schemes: combine a low-cost central compact scheme in smooth regions with WENO reconstruction only near detected discontinuities, cutting computational cost for mostly-smooth flows (e.g. turbulence with occasional shocks).
• Multi-dimensional systems: WENO is applied dimension-by-dimension on structured grids (finite-difference WENO) or via reconstruction on unstructured meshes (finite-volume WENO) for complex geometries such as supersonic aerofoils.
• Relation to LBM and Navier-Stokes solvers on this site: unlike Lattice-Boltzmann or the pressure-projection Navier-Stokes solver, WENO targets the fully compressible, shock-forming regime — the natural next step once flows are no longer safely incompressible.