L06 · Virtual Wind Tunnel

"""L06 Virtual Wind Tunnel: a D2Q9 lattice Boltzmann (BGK) solver in Numba, run through a V&V workflow.

  code verification      Poiseuille channel vs the exact parabola, observed order of accuracy
  solution verification  grid-convergence index (GCI) on the Schafer-Turek 2D-2 cylinder benchmark
  validation             unconfined cylinder at Re = 100: Strouhal number and drag vs Williamson 1996
  showcase               NACA 0012 angle-of-attack sweep (lift/drag polar) + a vorticity movie

Run from the repo root:   uv run python L06_virtual_wind_tunnel/sim.py          (full run, ~6 min)
                          uv run python L06_virtual_wind_tunnel/sim.py --quick  (smoke test, ~1.5 min)
Everything is in lattice units: dx = dt = 1, sound speed c_s = 1/sqrt(3), viscosity nu = (tau - 1/2) / 3.
"""
import os

os.environ.setdefault("NUMBA_NUM_THREADS", "3")  # shared machine: cap threads BEFORE numba loads
import sys
import time
from pathlib import Path

import imageio_ffmpeg
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FFMpegWriter
from matplotlib.colors import AsinhNorm, LinearSegmentedColormap
from matplotlib.patches import Polygon
from matplotlib.path import Path as MplPath
from numba import njit, prange, stencil
from scipy import fft

OUT = Path(__file__).parent / "output"
QUICK = "--quick" in sys.argv
plt.rcParams["animation.ffmpeg_path"] = imageio_ffmpeg.get_ffmpeg_exe()

# D2Q9 lattice: rest, 4 axis neighbours, 4 diagonals. OPP[i] is the direction opposite to i.
CX = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
CY = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
W = np.array([4 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 36, 1 / 36, 1 / 36, 1 / 36])
OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6])

BG = "#0b0e14"  # figure background
VORT = LinearSegmentedColormap.from_list(  # diverging, with a dark centre so still air is the background
    "vort", ["#e8f7ff", "#56cfe1", "#2463a6", "#12213d", BG, "#46142b", "#b8233f", "#f3722c", "#ffe8a3"])


# ----------------------------------------------------------------------------- the solver
def _make_stream_collide(forced):
    @njit(parallel=True, fastmath=True)
    def stream_collide(f, fn, solid, omega, gx):
        """One fused pull-stream + BGK-collide step on every interior fluid node (+ Guo body force gx if forced).

        Pull scheme: node (y, x) gathers population i from its upstream neighbour (y - cy_i, x - cx_i), relaxes
        the nine towards the local equilibrium and writes them to fn. A 3x3 stencil with no write conflicts, so
        rows are split over threads with prange. Measured speed notes: the nine directions written out by hand
        instead of a loop over i with a scratch array: 5x faster (values stay in registers). No branches inside
        the x loop beyond the solid test (an optional "also store rho, u" flag cost another 5x, because LLVM
        then gives up on AVX vectorisation). Solid neighbours hold bounce-back values put there by bounce().
        Fusing the sponge in here (one more multiply-add per population) was 5x slower: same vectoriser cliff.
        So was a 1024-wide grid: 4 kB rows put all 27 input streams in the same cache sets. Avoid powers of 2.
        """
        _, ny, nx = f.shape
        # ponytail: BGK single relaxation time, unstable as tau -> 1/2 (see limits()); TRT / MRT / regularised
        # collision relax the non-hydrodynamic moments separately and reach much higher Re on the same grid
        a = 1.0 - omega
        k = (1.0 - 0.5 * omega) * gx  # Guo et al. 2002 source prefactor
        for y in prange(ny):
            ym = y - 1 if y > 0 else ny - 1  # y is periodic; walls are rows of solid nodes
            yp = y + 1 if y < ny - 1 else 0
            for x in range(1, nx - 1):  # columns 0 and nx-1 are boundary/ghost columns, set by edges()
                if solid[y, x]:
                    continue
                f0 = f[0, y, x]
                f1 = f[1, y, x - 1]
                f2 = f[2, ym, x]
                f3 = f[3, y, x + 1]
                f4 = f[4, yp, x]
                f5 = f[5, ym, x - 1]
                f6 = f[6, ym, x + 1]
                f7 = f[7, yp, x + 1]
                f8 = f[8, yp, x - 1]
                rho = f0 + f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8
                ux = (f1 - f3 + f5 - f6 - f7 + f8) / rho
                uy = (f2 - f4 + f5 + f6 - f7 - f8) / rho
                if forced:
                    ux += 0.5 * gx / rho
                usq = 1.5 * (ux * ux + uy * uy)
                r0 = omega * rho * (4.0 / 9.0)
                r1 = omega * rho * (1.0 / 9.0)
                r2 = omega * rho * (1.0 / 36.0)
                # f_i <- f_i + omega (feq_i - f_i),  feq_i = w_i rho (1 + 3 c.u + 4.5 (c.u)^2 - 1.5 u^2)
                g0 = a * f0 + r0 * (1.0 - usq)
                c = 3.0 * ux
                g1 = a * f1 + r1 * (1.0 + c + 0.5 * c * c - usq)
                g3 = a * f3 + r1 * (1.0 - c + 0.5 * c * c - usq)
                c = 3.0 * uy
                g2 = a * f2 + r1 * (1.0 + c + 0.5 * c * c - usq)
                g4 = a * f4 + r1 * (1.0 - c + 0.5 * c * c - usq)
                c = 3.0 * (ux + uy)
                g5 = a * f5 + r2 * (1.0 + c + 0.5 * c * c - usq)
                g7 = a * f7 + r2 * (1.0 - c + 0.5 * c * c - usq)
                c = 3.0 * (uy - ux)
                g6 = a * f6 + r2 * (1.0 + c + 0.5 * c * c - usq)
                g8 = a * f8 + r2 * (1.0 - c + 0.5 * c * c - usq)
                if forced:  # + (1 - omega/2) S_i, S_i = w_i [3 (c_i - u) + 9 (c_i.u) c_i] . (gx, 0)
                    g0 -= (4.0 / 3.0) * ux * k
                    g1 += (1.0 / 9.0) * (3.0 + 6.0 * ux) * k
                    g3 += (1.0 / 9.0) * (6.0 * ux - 3.0) * k
                    g2 -= (1.0 / 3.0) * ux * k
                    g4 -= (1.0 / 3.0) * ux * k
                    g5 += (1.0 / 36.0) * (3.0 + 6.0 * ux + 9.0 * uy) * k
                    g7 += (1.0 / 36.0) * (6.0 * ux + 9.0 * uy - 3.0) * k
                    g6 += (1.0 / 36.0) * (6.0 * ux - 9.0 * uy - 3.0) * k
                    g8 += (1.0 / 36.0) * (3.0 + 6.0 * ux - 9.0 * uy) * k
                fn[0, y, x] = g0
                fn[1, y, x] = g1
                fn[2, y, x] = g2
                fn[3, y, x] = g3
                fn[4, y, x] = g4
                fn[5, y, x] = g5
                fn[6, y, x] = g6
                fn[7, y, x] = g7
                fn[8, y, x] = g8

    return stream_collide


# "forced" is a closure constant, so Numba compiles the force terms out of the tunnel version entirely
stream_collide = _make_stream_collide(False)  # wind-tunnel runs (float32: ~4x the throughput of float64)
stream_collide_forced = _make_stream_collide(True)  # body-force-driven channel (verification, float64)


@njit
def bounce(f, links, body):
    """Half-way bounce-back + momentum-exchange force (Ladd 1994), before each stream_collide.

    For every fluid->solid link, the population that left fluid node x towards the wall comes straight back:
    f_i(x, t+1) = f_opp(i)(x, t). Writing it into the solid node the pull reads from makes the main loop
    branch-free. Each link turns the momentum c_opp f around, handing 2 c_opp f to the body.
    """
    fx = 0.0
    fy = 0.0
    # ponytail: staircase mask + wall assumed half-way along every link = 1st order on curved walls;
    # Bouzidi et al. (2001) interpolated bounce-back uses the true wall distance for 2nd order
    for k in range(links.shape[0]):
        i, y, x, sy, sx = links[k, 0], links[k, 1], links[k, 2], links[k, 3], links[k, 4]
        v = f[OPP[i], y, x]
        f[i, sy, sx] = v
        if body[k]:
            fx -= 2.0 * CX[i] * v
            fy -= 2.0 * CY[i] * v
    return fx, fy


@njit
def feq(i, rho, ux, uy):
    cu = 3.0 * (CX[i] * ux + CY[i] * uy)
    return W[i] * rho * (1.0 + cu + 0.5 * cu * cu - 1.5 * (ux * ux + uy * uy))


@njit
def edges(f, inflow, solid):
    """Velocity inlet at x = 0 and pressure outlet (rho = 1) at x = nx-1, both by non-equilibrium extrapolation
    (Guo, Zheng & Shi 2002): the edge node gets the equilibrium of the prescribed state, taking the unknown
    quantity (inlet density, outlet velocity) and the non-equilibrium part f - feq from its inner neighbour.
    A first version pinned rho = 1 at the inlet and copied the outlet: nothing anchored the pressure, and the
    tunnel slowly pressurised while the flow rate sagged by 16 %."""
    _, ny, nx = f.shape
    for y in range(ny):
        if solid[y, 0]:
            continue
        for xb, xn in ((0, 1), (nx - 1, nx - 2)):
            rho = 0.0
            jx = 0.0
            jy = 0.0
            for i in range(9):
                rho += f[i, y, xn]
                jx += CX[i] * f[i, y, xn]
                jy += CY[i] * f[i, y, xn]
            ux, uy = jx / rho, jy / rho
            rb, uxb, uyb = (rho, inflow[y], 0.0) if xb == 0 else (1.0, ux, uy)
            for i in range(9):
                f[i, y, xb] = feq(i, rb, uxb, uyb) + f[i, y, xn] - feq(i, rho, ux, uy)


@njit(parallel=True, fastmath=True)
def sponge(f, sigma, x_start, target):
    """Absorbing layer: from column x_start[y] on, relax towards target[:, y] by a fraction sigma[y, x]."""
    _, ny, nx = f.shape
    for y in prange(ny):
        for i in range(9):
            t = target[i, y]
            for x in range(x_start[y], nx - 1):
                f[i, y, x] += sigma[y, x] * (t - f[i, y, x])


def sponge_profile(solid):
    """Absorbing layer strength: relax towards the undisturbed inflow state by a fraction sigma per step,
    0 -> 0.1 (quadratic) over the last 10 % of the length, and over 10 % of the height at the top and bottom
    unless those are walls. An absorbing far field around the test section.

    Why: LBM is weakly compressible, so a closed box is also an organ pipe. The impulsive start and the
    shedding itself launch sound waves that BGK at low viscosity barely damps; without this layer they bounce
    between the inlet, the outlet and the periodic sides and ride on top of the force signal."""
    # ponytail: plain relaxation sponge, still reflects a little and costs 10-20 % of the domain;
    # characteristic (non-reflecting) boundary conditions if the acoustics themselves matter
    ny, nx = solid.shape
    ramp = lambda n, m: 0.1 * np.clip(1 - np.arange(n) / m, 0, 1) ** 2  # 0.1 at index 0, zero beyond m
    sigma = np.repeat(ramp(nx, nx // 10)[None, ::-1], ny, 0)
    if not solid[0].all():
        side = ramp(ny, ny // 10)
        sigma = np.maximum(sigma, np.maximum(side, side[::-1])[:, None])
    return sigma.astype(np.float32), np.argmax(sigma > 0, axis=1)


@stencil
def _curl(ux, uy):
    return 0.5 * (uy[0, 1] - uy[0, -1] - ux[1, 0] + ux[-1, 0])  # [y, x] indexing: dv/dx - du/dy


@njit(parallel=True)
def vorticity(ux, uy):
    """numba.stencil writes the relative-index kernel once; Numba generates the parallel loop and zero borders."""
    return _curl(ux, uy)


def equilibrium(ux, uy, dtype):
    cu = 3 * (CX[:, None, None] * ux + CY[:, None, None] * uy)
    return (W[:, None, None] * (1 + cu + 0.5 * cu**2 - 1.5 * (ux**2 + uy**2))).astype(dtype)


def find_links(solid, body):
    """All (direction, fluid node, solid node) links where the pull of an interior node reads from a solid."""
    ny, nx = solid.shape
    fluid = ~solid
    fluid[:, [0, -1]] = False  # boundary/ghost columns are set by edges(), not streamed
    out = []
    for i in range(1, 9):
        src_solid = np.roll(solid, (CY[i], CX[i]), axis=(0, 1))  # src_solid[y, x] = solid[y - cy, x - cx]
        y, x = np.nonzero(fluid & src_solid)
        out.append(np.stack([np.full_like(y, i), y, x, (y - CY[i]) % ny, (x - CX[i]) % nx], 1))
    links = np.concatenate(out)
    return links, body[links[:, 3], links[:, 4]]


def moments(f, solid, gx=0.0):
    """Density and velocity from (post-collision) populations; zero velocity inside solids.
    With a Guo force the post-collision momentum is rho u + g/2, hence the -g/2."""
    rho = f.sum(0)
    ux = (np.tensordot(CX.astype(f.dtype), f, 1) - 0.5 * gx) / rho
    uy = np.tensordot(CY.astype(f.dtype), f, 1) / rho
    ux[solid] = 0
    uy[solid] = 0
    return rho, ux, uy


def simulate(solid, tau, steps, f, inflow=None, body=None, gx=0.0, every=0, on_frame=None):
    """Advance f by `steps`; returns (f, force history [steps, 2] on the body, lattice units).

    inflow=None: periodic in x (columns 0 and nx-1 are ghost copies), else the inlet u_x profile per row.
    body: which solids count towards the force (default: all). on_frame(t, f, F) is called every `every` steps.
    """
    # ponytail: one uniform grid everywhere; octree refinement (PowerFLOW "VR regions") puts cells where needed
    links, on_body = find_links(solid, solid if body is None else body)
    if inflow is not None:
        sigma, x_start = sponge_profile(solid)
        target = equilibrium(inflow[:, None], np.zeros((len(inflow), 1)), f.dtype)[:, :, 0]
    kernel = stream_collide_forced if gx else stream_collide
    fn = f.copy()
    F = np.empty((steps, 2))
    for t in range(steps):
        F[t] = bounce(f, links, on_body)
        kernel(f, fn, solid, 1.0 / tau, gx)
        if inflow is None:
            fn[:, :, 0], fn[:, :, -1] = fn[:, :, -2], fn[:, :, 1]
        else:
            edges(fn, inflow, solid)
            sponge(fn, sigma, x_start, target)
        f, fn = fn, f
        if every and t % every == 0:
            on_frame(t, f, F)
    return f, F


def spectrum(sig):
    """Hann-windowed, 16x zero-padded power spectrum (scipy.fft). Frequencies in cycles per time step."""
    s = (sig - sig.mean()) * np.hanning(len(sig))
    n = 16 * len(s)  # padding interpolates the spectrum between bins (it does not add resolution)
    return fft.rfftfreq(n), np.abs(fft.rfft(s, n)) ** 2


def zero_cross_period(sig):
    """Mean period from upward zero crossings (linear interpolation), an independent check on the FFT."""
    s = sig - sig.mean()
    i = np.flatnonzero((s[:-1] < 0) & (s[1:] >= 0))
    tc = i + s[i] / (s[i] - s[i + 1])
    return np.diff(tc).mean() if len(tc) > 2 else np.nan


# ----------------------------------------------------------------------------- geometry
def disk(ny, nx, xc, yc, d):
    """Solid mask of a circle. Put the centre between nodes (x.5) so half-way walls sit at exactly r on the axes."""
    yy, xx = np.mgrid[:ny, :nx]
    return (xx - xc) ** 2 + (yy - yc) ** 2 < (d / 2) ** 2


def naca00xx(alpha_deg, c, xp, yp, t=0.12, n=400):
    """Closed polygon of a symmetric NACA 4-digit section, chord c, pitched nose-up by alpha about the
    quarter chord placed at (xp, yp). Cosine spacing clusters points at the leading edge."""
    xs = 0.5 * (1 - np.cos(np.linspace(0, np.pi, n)))
    yt = 5 * t * (0.2969 * np.sqrt(xs) - 0.1260 * xs - 0.3516 * xs**2 + 0.2843 * xs**3 - 0.1036 * xs**4)
    x = np.r_[xs[::-1], xs[1:]] - 0.25
    y = np.r_[yt[::-1], -yt[1:]]
    a = np.radians(alpha_deg)
    return np.c_[xp + c * (x * np.cos(a) + y * np.sin(a)), yp + c * (-x * np.sin(a) + y * np.cos(a))]


def polygon_mask(poly, ny, nx):
    yy, xx = np.mgrid[:ny, :nx]
    return MplPath(poly).contains_points(np.c_[xx.ravel(), yy.ravel()]).reshape(ny, nx)


def tunnel_start(solid, inflow, kick, dtype):
    """Impulsive start: every fluid node at the inflow velocity, plus a cross-flow blob kick = (xc, yc, r, amp)
    that breaks the symmetry, so vortex shedding starts in a few periods instead of growing from round-off."""
    yy, xx = np.mgrid[:solid.shape[0], :solid.shape[1]]
    xc, yc, r, amp = kick
    uy = amp * np.exp(-((xx - xc) ** 2 + (yy - yc) ** 2) / r**2)
    return equilibrium(inflow[:, None] * ~solid, uy * ~solid, dtype)


# ----------------------------------------------------------------------------- 1. code verification
def poiseuille(N, tau, u_max=0.05, tol=1e-12):
    """Body-force driven channel, N fluid rows between half-way bounce-back walls, run to steady state."""
    solid = np.zeros((N + 2, 3), bool)  # one real column + 2 periodic ghosts: exactly unidirectional flow
    solid[[0, -1]] = True
    nu = (tau - 0.5) / 3
    gx = 8 * nu * u_max / N**2
    s = np.arange(N + 2) - 0.5  # distance from the lower wall, which sits half-way between rows 0 and 1
    exact = (gx / (2 * nu) * s * (N - s))[1:-1]
    f = equilibrium(np.zeros(solid.shape), np.zeros(solid.shape), np.float64)
    u_old = np.zeros(N)
    while True:
        f, _ = simulate(solid, tau, 2000, f, gx=gx)
        u = moments(f, solid, gx)[1][1:-1, 1]
        if np.abs(u - u_old).max() < tol * u_max:
            return s[1:-1], u, exact
        u_old = u


def verify_poiseuille():
    tau = 0.8
    Ns = [8, 16, 32, 64] if QUICK else [8, 16, 32, 64, 128]
    errs = []
    for N in Ns:
        s, u, ex = poiseuille(N, tau)
        errs.append(np.sqrt(np.sum((u - ex) ** 2) / np.sum(ex**2)))
    errs = np.array(errs)
    p = np.log2(errs[:-1] / errs[1:])
    print("[verify] Poiseuille, BGK tau = 0.8, relative L2 error of u(y):")
    for N, e in zip(Ns, errs):
        print(f"         N = {N:4d}   error = {e:.3e}")
    print(f"         observed order p = {', '.join(f'{q:.3f}' for q in p)}  (expected 2)")

    taus = np.array([0.55, 0.6, 0.7, 0.8, 0.9, 0.95, 1.0, 1.2, 1.5, 2.0])
    slip = []  # signed error at the centreline, N = 16
    for t in taus:
        s, u, ex = poiseuille(16, t)
        slip.append((u[len(u) // 2] - ex[len(ex) // 2]) / ex[len(ex) // 2])
    slip = np.array(slip)
    tau_magic = 0.5 + np.sqrt(3 / 16)
    fit = np.polyfit((taus - 0.5) ** 2, slip, 1)  # the BGK wall error is linear in (tau - 1/2)^2 ...
    tau_zero = 0.5 + np.sqrt(-fit[1] / fit[0])  # ... and vanishes at the "magic" tau
    print(f"[verify] wall error vanishes at tau = {tau_zero:.4f} (fit), theory 1/2 + sqrt(3/16) = {tau_magic:.4f}")

    fig, ax = plt.subplots(1, 3, figsize=(15, 4.4))
    s8, u8, e8 = poiseuille(8, tau)
    ss = np.linspace(0, 8, 200)
    ax[0].plot(ss, 4 * 0.05 * ss * (8 - ss) / 64, color="#adb5bd", lw=2, label="exact parabola")
    ax[0].plot(s8, u8, "o", color="#f3722c", ms=8, label="LBM, N = 8 fluid rows")
    ax[0].axvline(0, color="#56cfe1", lw=3)
    ax[0].axvline(8, color="#56cfe1", lw=3)
    ax[0].set(xlabel="y (lattice units, walls at 0 and N)", ylabel="u_x (lattice units)",
              title="Poiseuille profile: 8 cells already match")
    ax[0].legend(loc="lower center")
    ax[1].loglog(Ns, errs, "o-", color="#f3722c", lw=2, ms=8, label="BGK, tau = 0.8")
    ax[1].loglog(Ns, errs[0] * (Ns[0] / np.array(Ns)) ** 2, "--", color="#adb5bd", label="slope -2")
    for N, e, q in zip(Ns[1:], errs[1:], p):
        ax[1].annotate(f"p = {q:.2f}", (N, e), textcoords="offset points", xytext=(6, 6), fontsize=9)
    ax[1].set_xticks(Ns, [str(n) for n in Ns], minor=False)
    ax[1].minorticks_off()
    ax[1].set(xlabel="fluid rows N", ylabel="relative L2 error", title="Grid convergence: observed order p",
              xlim=(Ns[0] / 1.3, Ns[-1] * 1.5))
    ax[1].legend()
    ax[2].axhline(0, color="#adb5bd", lw=0.8)
    ax[2].plot(taus, slip * 100, "o-", color="#f3722c", lw=2, ms=7, label="LBM, N = 16")
    tt = np.linspace(0.52, 2.0, 200)
    ax[2].plot(tt, np.polyval(fit, (tt - 0.5) ** 2) * 100, ":", color="#adb5bd",
               label="fit  a (tau - 1/2)^2 + b")
    ax[2].axvline(tau_magic, color="#56cfe1", ls="--", label="theory zero, tau = 1/2 + sqrt(3/16)")
    ax[2].set(xlabel="relaxation time tau", ylabel="centre-line velocity error (%)",
              title="Bounce-back wall error depends on viscosity")
    ax[2].legend(fontsize=8)
    for a in ax:
        a.grid(alpha=0.2)
    fig.tight_layout()
    fig.savefig(OUT / "verify_poiseuille.png", dpi=130)
    plt.close(fig)
    return dict(Ns=Ns, errs=errs, p=p, tau_zero=tau_zero, tau_magic=tau_magic)


# ----------------------------------------------------------------------------- 2. cylinder runs
def cylinder_run(D, U, Re, ny, nx, xc, yc, inflow, walls, n_conv, n_meas, dtype=np.float32):
    """Cylinder of diameter D (cells) at Re = U D / nu. Runs n_conv + n_meas convective times D/U and
    returns force coefficients, sampled every step, over the last n_meas of them."""
    body = disk(ny, nx, xc, yc, D)
    solid = body.copy()
    if walls:
        solid[[0, -1]] = True
    nu = U * D / Re
    tau = 3 * nu + 0.5
    steps = int((n_conv + n_meas) * D / U)
    f = tunnel_start(solid, inflow, (xc + D, yc, D, 0.5 * U), dtype)
    t0 = time.perf_counter()
    f, F = simulate(solid, tau, steps, f, inflow=inflow, body=body)
    wall = time.perf_counter() - t0
    m = int(n_meas * D / U)
    cd, cl = (2 * F[-m:] / (U**2 * D)).T  # rho0 = 1
    fr, P = spectrum(cl)
    j = np.argmax(P)
    j += 0.5 * (P[j - 1] - P[j + 1]) / (P[j - 1] - 2 * P[j] + P[j + 1])  # parabolic peak interpolation
    st = j * (fr[1] - fr[0]) * D / U
    return dict(D=D, U=U, tau=tau, cd=cd, cl=cl, fr=fr * D / U, P=P, St=st, St_zc=D / U / zero_cross_period(cl),
                cd_hist=2 * F[:, 0] / (U**2 * D), cl_hist=2 * F[:, 1] / (U**2 * D), f=f, solid=solid,
                mlups=solid.size * steps / wall / 1e6, blowup=not np.isfinite(F).all())


def schafer_turek(D, n_conv=35, n_meas=15):
    """Schafer & Turek (1996) test 2D-2: channel H = 4.1 D, cylinder centre 2D from the inlet and 2D above the
    lower wall (slightly off-centre), parabolic inflow, Re = U_mean D / nu = 100, blockage D/H = 24 %."""
    H = round(4.1 * D)
    ny, nx = H + 2, 22 * D + 1
    s = np.arange(ny) - 0.5
    Um = 0.1  # peak inflow velocity, Ma = 0.17
    inflow = np.clip(4 * Um * s * (H - s) / H**2, 0, None)
    r = cylinder_run(D, 2 * Um / 3, 100, ny, nx, 2 * D + 0.5, 2 * D + 0.5, inflow, True, n_conv, n_meas)
    r.update(cd_max=r["cd"].max(), cl_max=r["cl"].max(), cd_mean=r["cd"].mean())
    return r


def open_cylinder(up=20, H=20, D=20, U=0.08, n_conv=15, n_meas=25):
    """Unconfined cylinder at Re = 100 for validation against experiment: uniform inflow `up` diameters
    upstream, 20 D of wake behind it, periodic sides H diameters apart (an infinite column of cylinders,
    blockage D/H) with absorbing bands along them."""
    ny, nx = H * D, (up + 20) * D
    r = cylinder_run(D, U, 100, ny, nx, up * D + 0.5, ny / 2 + 0.5, np.full(ny, U), False, n_conv, n_meas)
    r.update(xc=up * D, yc=ny // 2, cl_amp=0.5 * (r["cl"].max() - r["cl"].min()))
    return r


# ----------------------------------------------------------------------------- 3. showcase: NACA 0012
def airfoil_tunnel(alpha, c, Re, U, ny, nx, steps, every=0, on_frame=None, f0=None, xp=2.0):
    """NACA 0012 at angle of attack alpha, chord c cells, quarter chord xp chords behind the inlet (continue
    from f0 if given). Returns polygon, solid mask, final f and C_L, C_D histories (lift +y, drag +x)."""
    poly = naca00xx(alpha, c, xp * c, ny / 2 - 0.5)
    solid = polygon_mask(poly, ny, nx)
    tau = 3 * U * c / Re + 0.5
    inflow = np.full(ny, U)
    f = tunnel_start(solid, inflow, (0, 0, 1, 0.0), np.float32) if f0 is None else f0  # no kick: alpha breaks symmetry
    t0 = time.perf_counter()
    f, F = simulate(solid, tau, steps, f, inflow=inflow, every=every, on_frame=on_frame)
    coef = 2 * F / (U**2 * c)
    return dict(poly=poly, solid=solid, tau=tau, cd=coef[:, 0], cl=coef[:, 1], f=f,
                mlups=solid.size * steps / (time.perf_counter() - t0) / 1e6, blowup=not np.isfinite(F).all())


def sweep_alphas():
    return np.array([0, 10, 20.0]) if QUICK else np.arange(0, 20.1, 2.5)


def airfoil_sweep(c=64, Re=1000, U=0.08, n_conv=10, n_meas=6):
    """Angle-of-attack sweep at Re = 1000 (insect / micro-air-vehicle regime, the laptop-LBM ceiling).
    Tunnel 10c x 5c, quarter chord 3.5c behind the inlet. Averages over the last n_meas convective times c/U."""
    ny, nx = 5 * c, 10 * c
    steps, m = int((n_conv + n_meas) * c / U), int(n_meas * c / U)
    rows, snaps = [], []
    for a in sweep_alphas():
        r = airfoil_tunnel(a, c, Re, U, ny, nx, steps, xp=3.5)
        cl, cd = r["cl"][-m:], r["cd"][-m:]
        st = np.nan
        if cl.std() > 0.02:  # unsteady: shedding frequency from the lift spectrum
            fr, P = spectrum(cl)
            st = fr[np.argmax(P)] * c / U
        rows.append((a, cl.mean(), cl.std(), cd.mean(), cd.std(), st, r["tau"], r["mlups"], r["blowup"]))
        snaps.append((a, r["f"], r["solid"], r["poly"]))
        print(f"[airfoil] alpha = {a:4.1f} deg: C_L = {cl.mean():.3f} +- {cl.std():.3f}, C_D = {cd.mean():.3f}, "
              f"L/D = {cl.mean() / cd.mean():5.2f}" + (f", St_c = {st:.3f}" if st == st else ", steady")
              + f"  ({r['mlups']:.0f} MLUPS{', BLOW-UP' if r['blowup'] else ''})")
    rows = np.array(rows, dtype=float)
    # Domain check (the cylinder taught us to distrust a nearby inlet): repeat two angles in a tunnel 1.5x as
    # long and as tall. A first sweep in a 7c x 3.5c tunnel was 6-8 % high in C_L and C_D by this test.
    check = []
    for a in ([10.0] if QUICK else [5.0, 15.0]):
        r = airfoil_tunnel(a, c, Re, U, 3 * ny // 2, 3 * nx // 2, steps, xp=5.25)
        cl, cd = r["cl"][-m:].mean(), r["cd"][-m:].mean()
        cl0, cd0 = rows[rows[:, 0] == a][0, [1, 3]]
        check.append((a, cl, cd))
        print(f"[airfoil] domain check alpha = {a:g}: 10c x 5c tunnel C_L = {cl0:.3f}, C_D = {cd0:.3f};  "
              f"15c x 7.5c tunnel C_L = {cl:.3f} ({100 * (cl / cl0 - 1):+.1f} %), C_D = {cd:.3f} ({100 * (cd / cd0 - 1):+.1f} %)")
    return dict(rows=rows, snaps=snaps, c=c, U=U, Re=Re, check=np.array(check))


# ----------------------------------------------------------------------------- limits of the method
def limits(D=20):
    """Two honest ceilings. (1) BGK stability: cylinder at D = 20 cells, U = 0.1, raise Re (tau -> 1/2) until
    it blows up. (2) Compressibility: Schafer-Turek at D = 20 for four Mach numbers, C_D against Ma^2."""
    probe = []
    U = 0.1
    ny, nx = 8 * D, 16 * D
    for Re in ([400, 3200] if QUICK else [200, 400, 800, 1600, 3200, 6400]):
        r = cylinder_run(D, U, Re, ny, nx, 4 * D + 0.5, ny / 2 + 0.5, np.full(ny, U), False, 25, 5)
        ok = not r["blowup"] and np.abs(r["cd"]).max() < 50
        probe.append((Re, r["tau"], ok))
        print(f"[limits] D = {D}, U = {U}, Re = {Re:5d}: tau = {r['tau']:.5f} -> {'stable' if ok else 'BLOWS UP'}")
    mach = []
    for Um in ([0.1, 0.15] if QUICK else [0.05, 0.075, 0.1, 0.15]):
        H = round(4.1 * D)
        s = np.arange(H + 2) - 0.5
        r = cylinder_run(D, 2 * Um / 3, 100, H + 2, 22 * D + 1, 2 * D + 0.5, 2 * D + 0.5,
                         np.clip(4 * Um * s * (H - s) / H**2, 0, None), True, 35, 15)
        ma = Um / np.sqrt(1 / 3)
        mach.append((ma, r["cd"].mean(), r["St"], r["tau"]))
        print(f"[limits] Schafer-Turek D = {D}: peak Ma = {ma:.3f}, tau = {r['tau']:.4f}: "
              f"C_D mean = {r['cd'].mean():.4f}, St = {r['St']:.4f}")
    mach = np.array(mach)
    fit = np.polyfit(mach[:, 0] ** 2, mach[:, 1], 1)
    fit_st = np.polyfit(mach[:, 0] ** 2, mach[:, 2], 1)
    print(f"[limits] C_D(Ma) ~ {fit[1]:.4f} + {fit[0]:.3f} Ma^2  ->  at Ma = 0.17 the compressibility error is "
          f"{fit[0] * 0.173**2 / fit[1] * 100:+.2f} %")

    fig, ax = plt.subplots(1, 2, figsize=(12, 4.4))
    for Re, tau, ok in probe:
        ax[0].plot(Re, tau - 0.5, "o" if ok else "X", ms=11, color="#56cfe1" if ok else "#f3722c")
    ax[0].set(xscale="log", yscale="log", xlabel="Reynolds number U D / nu", ylabel="tau - 1/2 = 3 nu",
              title=f"BGK stability, cylinder D = {D} cells, U = {U} (o stable, x blew up)")
    ax[0].set_xticks([p[0] for p in probe], [str(p[0]) for p in probe])
    ax[0].minorticks_off()
    mm = np.linspace(0, mach[:, 0].max() * 1.05, 50)
    ax[1].plot(mach[:, 0] ** 2, mach[:, 1], "o", color="#f3722c", ms=9, label="LBM, D = 20")
    ax[1].plot(mm**2, np.polyval(fit, mm**2), "--", color="#adb5bd", label=f"fit {fit[1]:.3f} + {fit[0]:.2f} Ma^2")
    ax[1].set(xlabel="Ma^2 (peak inflow)", ylabel="mean C_D", title="Compressibility error grows like Ma^2")
    ax[1].legend()
    for a in ax:
        a.grid(alpha=0.2, which="both")
    fig.tight_layout()
    fig.savefig(OUT / "limits.png", dpi=130)
    plt.close(fig)
    return dict(probe=probe, mach=mach, fit=fit, fit_st=fit_st)


def error_budget(val, lim, U=0.08):
    """Validation error budget of the largest open-cylinder tunnel: remove the discretisation error (the
    D = 20 vs Richardson ratio of the Schafer-Turek grid study: same D, same Re, a transferred estimate) and
    the compressibility error (the Ma^2 fit at this run's Ma), then compare what is left with the references."""
    if not val["study"]:
        return
    ma2 = (U / np.sqrt(1 / 3)) ** 2
    d20 = val["runs"][0]
    grid_cd = d20["cd_mean"] / val["study"]["cd_mean"][1]
    grid_st = d20["St"] / val["study"]["St"][1]
    mach_cd = np.polyval(lim["fit"], ma2) / lim["fit"][1]
    mach_st = np.polyval(lim["fit_st"], ma2) / lim["fit_st"][1]
    big = val["big"]
    cd = big["cd"].mean() / grid_cd / mach_cd
    st = big["St"] / grid_st / mach_st
    print(f"[budget] C_D: {big['cd'].mean():.4f} (largest tunnel) -> /{grid_cd:.4f} (grid) -> /{mach_cd:.4f} "
          f"(Ma = {np.sqrt(ma2):.3f}) = {cd:.4f}   refs {CD_REF[0]}-{CD_REF[1]}")
    print(f"[budget] St : {big['St']:.4f} (largest tunnel) -> /{grid_st:.4f} (grid) -> /{mach_st:.4f} (Ma) = {st:.4f}"
          f"   Williamson {WILLIAMSON_ST} ({100 * (st / WILLIAMSON_ST - 1):+.1f} %)")


# ----------------------------------------------------------------------------- figures and movie
def vort_field(f, solid, U, L):
    """Non-dimensional vorticity omega L / U, via the numba.stencil curl."""
    _, ux, uy = moments(f, solid)
    return vorticity(ux, uy) * (L / U)


def draw_field(ax, w, patch, lim=12.0, lw=0.03):
    """Vorticity on a dark background, asinh-scaled so the weak far wake and the strong boundary layer both
    show; the body is a crisp vector patch drawn over the staircase mask."""
    ny, nx = w.shape
    im = ax.imshow(w, origin="lower", cmap=VORT, norm=AsinhNorm(lw * lim, -lim, lim),
                   extent=(-0.5, nx - 0.5, -0.5, ny - 0.5), interpolation="bilinear")
    ax.add_patch(patch)
    ax.set_axis_off()
    return im


def body_patch(poly=None, circle=None):
    kw = dict(facecolor="#d8dee9", edgecolor="white", lw=0.8, zorder=3, antialiased=True)
    return Polygon(poly, closed=True, **kw) if poly is not None else plt.Circle(circle[:2], circle[2], **kw)


def gci(h, phi):
    """Grid Convergence Index, Celik et al. (2008, J. Fluids Eng. 130, 078001). Three grids, finest first.
    Returns observed order p, Richardson-extrapolated value, GCI_fine (fraction), convergence type."""
    (h1, h2, h3), (f1, f2, f3) = h, phi
    r21, r32 = h2 / h1, h3 / h2
    e21, e32 = f2 - f1, f3 - f2
    s = np.sign(e32 / e21)  # +1 monotone, -1 oscillatory convergence
    p = 1.0
    for _ in range(200):  # fixed-point iteration, needed because r21 != r32
        p = abs(np.log(abs(e32 / e21)) + np.log((r21**p - s) / (r32**p - s))) / np.log(r21)
    ext = (r21**p * f1 - f2) / (r21**p - 1)
    return p, ext, 1.25 * abs(e21 / f1) / (r21**p - 1), "monotone" if s > 0 else "oscillatory"


ST_REF = {"cd_max": (3.22, 3.24), "cl_max": (0.99, 1.01), "St": (0.295, 0.305)}  # Schafer & Turek 1996, 2D-2
WILLIAMSON_ST = 0.164  # Williamson 1996 (parallel shedding, Re = 100, experiment)
CD_REF = (1.33, 1.35)  # unconfined numerics: Park, Kwon & Choi 1998 (1.33); Liu, Zheng & Sung 1998 (1.35)
CL_AMP_REF = (0.33, 0.34)  # same two papers: lift amplitude 0.332 / 0.339


def validate_cylinders():
    Ds = [20, 30] if QUICK else [20, 30, 40]
    runs = []
    for D in Ds:
        r = schafer_turek(D)
        runs.append(r)
        print(f"[ST 2D-2] D = {D:2d} (tau = {r['tau']:.3f}, {r['mlups']:.0f} MLUPS): C_D,max = {r['cd_max']:.4f}, "
              f"C_L,max = {r['cl_max']:.4f}, St = {r['St']:.4f} (zero-crossing {r['St_zc']:.4f}), "
              f"C_D,mean = {r['cd_mean']:.4f}")
    study = {}
    if len(Ds) == 3:
        h = [1 / D for D in Ds[::-1]]
        for q in ("cd_max", "cl_max", "St", "cd_mean"):
            p, ext, g, kind = gci(h, [r[q] for r in runs[::-1]])
            study[q] = (p, ext, g, kind)
            ref = ST_REF.get(q)
            print(f"[GCI] {q:7s}: p = {p:.2f} ({kind}), Richardson extrapolation {ext:.4f}, GCI_fine = {g * 100:.2f} %"
                  + (f"   benchmark {ref[0]}-{ref[1]}" if ref else ""))
    # Domain-size study: a free-stream experiment has no inlet and no walls; the tunnel has both.
    # Move one boundary at a time until the answer stops changing. (A first version moved only the side
    # walls and blamed them for a +5 % Strouhal error; the inlet, 6 D upstream, turned out to be the culprit.)
    dom = {}
    for up, H in ([(6, 20), (20, 20)] if QUICK else [(6, 20), (12, 20), (20, 20), (20, 10), (20, 40)]):
        r = open_cylinder(up, H)
        dom[up, H] = r
        print(f"[validate] open cylinder Re = 100, D = 20, inlet {up:2d} D upstream, sides {H:2d} D apart "
              f"({r['mlups']:.0f} MLUPS): St = {r['St']:.4f}, C_D = {r['cd'].mean():.4f}, C_L amplitude = {r['cl_amp']:.3f}")
    oc, big = dom[20, 20], dom[max(dom, key=lambda k: k[0] * 100 + k[1])]
    print(f"[validate] largest tunnel vs literature: St = {big['St']:.4f} (Williamson 1996: {WILLIAMSON_ST}, "
          f"{100 * (big['St'] / WILLIAMSON_ST - 1):+.1f} %), C_D = {big['cd'].mean():.4f} (refs {CD_REF[0]}-{CD_REF[1]}), "
          f"C_L amplitude = {big['cl_amp']:.3f} (refs {CL_AMP_REF[0]}-{CL_AMP_REF[1]})")

    # --- figure: Schafer-Turek flow + grid convergence
    fine = runs[-1]
    D = fine["D"]
    fig = plt.figure(figsize=(15, 7.6), facecolor=BG)
    ax = fig.add_axes((0.02, 0.56, 0.96, 0.4))
    w = vort_field(fine["f"], fine["solid"], fine["U"], D)
    ny, nx = w.shape
    draw_field(ax, w, body_patch(circle=(2 * D + 0.5, 2 * D + 0.5, D / 2)))
    for yw in (0, ny - 1):
        ax.add_patch(plt.Rectangle((-0.5, yw - 0.5), nx, 1, color="#d8dee9", zorder=3))
    ax.set_title(f"Schafer-Turek 2D-2: Re = 100 in a channel, blockage 24 %, D = {D} cells "
                 f"(the sponge zone is the last 10 %)", color="w")
    labels = {"cd_max": "C_D,max", "cl_max": "C_L,max", "St": "Strouhal St"}
    for k, q in enumerate(labels):
        a = fig.add_axes((0.06 + 0.32 * k, 0.08, 0.26, 0.38))
        hh = np.array([1 / r["D"] for r in runs])
        a.axhspan(*ST_REF[q], color="#56cfe1", alpha=0.25, label="benchmark range")
        a.plot(hh, [r[q] for r in runs], "o-", color="#f3722c", ms=8, lw=2, label="LBM, D = 20, 30, 40")
        if q in study:
            p, ext, g, kind = study[q]
            a.plot([0, hh[-1]], [ext, runs[-1][q]], ":", color="#f3722c")
            a.plot(0, ext, "*", ms=14, color="#ffe8a3", label=f"Richardson, p = {p:.2f}")
        a.set(xlabel="grid spacing h = 1/D", ylabel=labels[q], xlim=(-0.002, 0.055))
        a.grid(alpha=0.2)
        a.legend(fontsize=8)
    fig.savefig(OUT / "validate_schafer_turek.png", dpi=120, facecolor=BG)
    plt.close(fig)

    # --- figure: open cylinder, wake + forces + spectrum
    fig = plt.figure(figsize=(15, 8.4), facecolor=BG)
    ax = fig.add_axes((0.01, 0.47, 0.98, 0.47))
    w = vort_field(oc["f"], oc["solid"], oc["U"], oc["D"])
    D, xc, yc = oc["D"], oc["xc"], oc["yc"]
    draw_field(ax, w[yc - 3 * D: yc + 3 * D, xc - 2 * D: xc + 15 * D],
               body_patch(circle=(2 * D + 0.5, 3 * D + 0.5, D / 2)), lim=6)
    ax.set_title("Karman vortex street behind a cylinder, Re = 100 (crop of the 40D x 20D tunnel)", color="w")
    a1 = fig.add_axes((0.05, 0.07, 0.27, 0.3))
    t = np.arange(len(oc["cl_hist"])) * oc["U"] / D
    a1.plot(t, oc["cd_hist"], color="#f3722c", lw=1.2, label="C_D")
    a1.plot(t, oc["cl_hist"], color="#56cfe1", lw=1.2, label="C_L")
    a1.axvspan(t[-len(oc["cl"])], t[-1], color="w", alpha=0.07, label="averaging window")
    a1.set(xlabel="t U / D", ylabel="force coefficient", ylim=(-1, 2.5))
    a1.legend(ncol=3, fontsize=8, loc="upper right")
    a2 = fig.add_axes((0.38, 0.07, 0.27, 0.3))
    a2.semilogy(oc["fr"], oc["P"] / oc["P"].max(), color="#56cfe1")
    a2.axvline(WILLIAMSON_ST, color="#ffe8a3", ls="--", label=f"Williamson 1996: St = {WILLIAMSON_ST}")
    a2.axvline(oc["St"], color="#f3722c", ls=":", label=f"LBM: St = {oc['St']:.3f}")
    a2.set(xlim=(0, 1), ylim=(1e-6, 2), xlabel="frequency f D / U", ylabel="lift power (norm.)")
    a2.legend(fontsize=8)
    a3 = fig.add_axes((0.71, 0.07, 0.27, 0.3))
    for keys, dist, ls, name in (([k for k in dom if k[1] == 20], lambda k: k[0], "-", "inlet"),
                                 ([k for k in dom if k[0] == 20], lambda k: k[1] / 2, "--", "side")):
        keys = sorted(keys, key=dist)
        x = [dist(k) for k in keys]
        a3.plot(x, [dom[k]["St"] / WILLIAMSON_ST for k in keys], "o" + ls, color="#56cfe1", label=f"St, move {name}")
        a3.plot(x, [dom[k]["cd"].mean() / np.mean(CD_REF) for k in keys], "s" + ls, color="#f3722c",
                label=f"C_D, move {name}")
    a3.axhline(1, color="#adb5bd", lw=0.8)
    a3.set(xscale="log", xlabel="distance from cylinder to that boundary (D)",
           ylabel="LBM / reference (St 0.164, C_D 1.34)", title="Domain-size study: inlet matters most")
    a3.set_xticks([5, 6, 10, 12, 20], ["5", "6", "10", "12", "20"])
    a3.minorticks_off()
    a3.legend(fontsize=7, ncol=2)
    for a in (a1, a2, a3):
        a.grid(alpha=0.2)
    fig.savefig(OUT / "validate_cylinder.png", dpi=120, facecolor=BG)
    plt.close(fig)
    return dict(runs=runs, study=study, oc=oc, big=big, dom=dom)


def polar_figures(sw):
    rows, c, U = sw["rows"], sw["c"], sw["U"]
    a, cl, cls, cd, cds, st = rows[:, :6].T
    fig, ax = plt.subplots(1, 3, figsize=(16, 4.8))
    aa = np.linspace(0, 12, 20)
    ax[0].plot(aa, 2 * np.pi * np.radians(aa), "--", color="#adb5bd", label="thin-airfoil theory 2 pi alpha")
    ax[0].errorbar(a, cl, yerr=cls, fmt="o-", color="#56cfe1", ms=7, capsize=3, lw=2,
                   label=f"LBM, Re = {sw['Re']} (bars: unsteady s.d.)")
    ck = sw["check"]
    ax[0].plot(ck[:, 0], ck[:, 1], "D", ms=9, mfc="none", mec="w", label="same, 1.5x larger tunnel")
    ax[1].plot(ck[:, 2], ck[:, 1], "D", ms=9, mfc="none", mec="w")
    ax[0].set(xlabel="angle of attack alpha (deg)", ylabel="lift coefficient C_L", title="Lift curve")
    ax[1].errorbar(cd, cl, xerr=cds, yerr=cls, fmt="o-", color="#f3722c", ms=7, capsize=3, lw=2)
    for ai, x, y in zip(a, cd, cl):
        ax[1].annotate(f"{ai:g} deg", (x, y), textcoords="offset points", xytext=(7, -4), fontsize=8, color="#adb5bd")
    ax[1].set(xlabel="drag coefficient C_D", ylabel="lift coefficient C_L", title="Drag polar")
    ax[2].plot(a, cl / cd, "o-", color="#ffe8a3", ms=7, lw=2)
    ax[2].set(xlabel="angle of attack alpha (deg)", ylabel="L / D", title="Glide ratio")
    for x in ax:
        x.grid(alpha=0.2)
    ax[0].legend(fontsize=8)
    fig.suptitle(f"NACA 0012 at Re = {sw['Re']}, chord {c} cells, D2Q9 BGK lattice Boltzmann", color="w")
    fig.tight_layout()
    fig.savefig(OUT / "airfoil_polar.png", dpi=130)
    plt.close(fig)

    n = len(sw["snaps"])
    cols = 3
    fig, axs = plt.subplots((n + cols - 1) // cols, cols, figsize=(15, 2.35 * ((n + cols - 1) // cols)),
                            facecolor=BG, squeeze=False)
    for axx, (al, f, solid, poly) in zip(axs.flat, sw["snaps"]):
        w = vort_field(f, solid, U, c)
        y0, y1, x0, x1 = int(1.5 * c), int(3.5 * c), int(3 * c), int(7.8 * c)
        draw_field(axx, w[y0:y1, x0:x1], body_patch(poly - [x0, y0]))
        axx.text(0.02, 0.86, f"alpha = {al:g} deg", transform=axx.transAxes, color="w", fontsize=10)
    for axx in axs.flat[n:]:
        axx.set_axis_off()
    fig.tight_layout(pad=0.4)
    fig.savefig(OUT / "airfoil_gallery.png", dpi=110, facecolor=BG)
    plt.close(fig)


def showcase_movie(alpha=20.0, c=128, Re=1000, U=0.08, frames=480, every=20):
    """High-resolution run of the stalled airfoil: vorticity MP4 (1280 x 720, 30 fps) + the hero still.
    The 1000-wide tunnel (not 1024: see the stream_collide docstring) is 7.8 chords long, 3.5 high."""
    if QUICK:
        c, frames = 64, 60
    ny, nx = int(3.5 * c), int(7.8125 * c)
    warm = int(10 * c / U)  # spin-up 10 c/U, then film 6 c/U: the same windows as the sweep
    if not QUICK:  # grid check: the same tunnel (in chords) at half the resolution
        h = airfoil_tunnel(alpha, 64, Re, U, int(3.5 * 64), int(7.8125 * 64), int(16 * 64 / U))
        m = int(6 * 64 / U)
        grid64 = (h["cl"][-m:].mean(), h["cd"][-m:].mean())
    r = airfoil_tunnel(alpha, c, Re, U, ny, nx, warm)  # spin-up, then film
    solid, poly = r["solid"], r["poly"]
    fig = plt.figure(figsize=(12.8, 7.2), dpi=100, facecolor=BG)
    ax = fig.add_axes((0, 0.2, 1, 0.8))
    im = draw_field(ax, np.zeros((ny, nx)), body_patch(poly))
    y0, y1 = 0.13 * ny, 0.87 * ny  # the test section: absorbing bands (10 %) and their edge shear layers left out
    x0 = 0.85 * c
    ax.set(xlim=(x0, x0 + (y1 - y0) * 1280 / 576), ylim=(y0, y1))
    ax.text(0.012, 0.93, f"NACA 0012  |  alpha = {alpha:g} deg  |  Re = {Re}", transform=ax.transAxes,
            color="w", fontsize=16)
    ax.text(0.012, 0.87, "vorticity, D2Q9 lattice Boltzmann (BGK) in Numba", transform=ax.transAxes,
            color="#8d99ae", fontsize=11)
    clock = ax.text(0.988, 0.93, "", transform=ax.transAxes, color="w", fontsize=14, ha="right", family="monospace")
    tr = fig.add_axes((0.06, 0.07, 0.9, 0.11), facecolor=BG)
    tt = np.arange(frames * every) / (c / U)
    l1, = tr.plot([], [], color="#56cfe1", lw=1.5, label="C_L")
    l2, = tr.plot([], [], color="#f3722c", lw=1.5, label="C_D")
    tr.set(xlim=(0, tt[-1]), ylim=(0, 2.4))
    tr.set_xlabel("t U / c (since filming started)", color="#8d99ae", fontsize=9, labelpad=1)
    tr.tick_params(colors="#8d99ae", labelsize=8)
    for sp in tr.spines.values():
        sp.set_color("#39414f")
    tr.legend(loc="upper right", ncol=2, fontsize=9, frameon=False, labelcolor="w")
    writer = FFMpegWriter(fps=30, bitrate=6000, extra_args=["-pix_fmt", "yuv420p"])
    t_warm = warm / (c / U)

    def on_frame(t, f, F):
        im.set_data(vort_field(f, solid, U, c))
        coef = 2 * F[: t + 1] / (U**2 * c)
        l1.set_data(tt[: t + 1], coef[:, 1])
        l2.set_data(tt[: t + 1], coef[:, 0])
        clock.set_text(f"t U/c = {t_warm + tt[t]:5.2f}")
        writer.grab_frame(facecolor=BG)

    with writer.saving(fig, str(OUT / "wind_tunnel.mp4"), dpi=100):
        r2 = airfoil_tunnel(alpha, c, Re, U, ny, nx, frames * every, every=every, on_frame=on_frame, f0=r["f"])
    plt.close(fig)
    cl, cd = r2["cl"], r2["cd"]
    print(f"[movie] alpha = {alpha:g}, chord {c}: C_L = {cl.mean():.3f} +- {cl.std():.3f}, C_D = {cd.mean():.3f} "
          f"({r2['mlups']:.0f} MLUPS incl. rendering)")
    if not QUICK:
        print(f"[movie] grid check, 7.8c x 3.5c tunnel, alpha = {alpha:g}: chord 64 -> C_L = {grid64[0]:.3f}, "
              f"C_D = {grid64[1]:.3f};  chord {c} -> C_L = {cl.mean():.3f} ({100 * (cl.mean() / grid64[0] - 1):+.1f} %), "
              f"C_D = {cd.mean():.3f} ({100 * (cd.mean() / grid64[1] - 1):+.1f} %)")

    w = vort_field(r2["f"], solid, U, c)  # hero still: same state, rendered big
    x1 = 0.9 * nx  # outlet sponge starts here
    Wp, Hp = 2400, int(2400 * (y1 - y0) / (x1 - x0))
    fig = plt.figure(figsize=(Wp / 150, Hp / 150), dpi=150, facecolor=BG)
    ax = fig.add_axes((0, 0, 1, 1))
    draw_field(ax, w, body_patch(poly))
    ax.set(xlim=(x0, x1), ylim=(y0, y1))
    ax.text(0.015, 0.93, "Virtual wind tunnel", transform=ax.transAxes, color="w", fontsize=24)
    ax.text(0.015, 0.875, f"NACA 0012 at {alpha:g} deg, Re = {Re}: the stalled wing sheds a vortex street. "
            f"Lattice Boltzmann, {ny * nx / 1e3:.0f}k cells.", transform=ax.transAxes, color="#8d99ae", fontsize=12)
    fig.savefig(OUT / "hero.png", dpi=150, facecolor=BG)
    plt.close(fig)
    return dict(cl=cl.mean(), cl_sd=cl.std(), cd=cd.mean(), c=c)


def verify(ver, val):
    """Assertions: the verification must pass; the validation must land in the right physical ballpark."""
    assert np.all(np.abs(ver["p"] - 2) < 0.1), f"Poiseuille order {ver['p']} is not 2"
    assert ver["errs"][1] < 0.01, "Poiseuille error at N = 16 above 1 %"
    assert abs(ver["tau_zero"] - ver["tau_magic"]) < 0.01, "bounce-back wall error does not vanish at the magic tau"
    st_fine = val["runs"][-1]["St"]
    assert 0.295 <= st_fine <= 0.305, f"Schafer-Turek St {st_fine} outside the benchmark range"
    big = val["big"]
    assert 0.16 <= big["St"] <= 0.20, f"open-cylinder St {big['St']} outside 0.16-0.20"
    assert abs(big["St"] / WILLIAMSON_ST - 1) < 0.03, f"open-cylinder St {big['St']} not within 3 % of Williamson"
    assert 1.2 <= big["cd"].mean() <= 1.6, "open-cylinder C_D out of the plausible range"
    print("[verify] all assertions passed")


def main():
    OUT.mkdir(exist_ok=True)
    plt.style.use("dark_background")
    plt.rcParams.update({"figure.facecolor": BG, "axes.facecolor": BG, "savefig.facecolor": BG})
    T0 = time.perf_counter()
    stamp = lambda what: print(f"[time] {what} done, {time.perf_counter() - T0:.0f} s elapsed")
    ver = verify_poiseuille()
    stamp("code verification (Poiseuille)")
    val = validate_cylinders()
    stamp("cylinder verification + validation")
    verify(ver, val)
    lim = limits()
    error_budget(val, lim)
    stamp("limits (stability, Mach)")
    sw = airfoil_sweep()
    polar_figures(sw)
    stamp("airfoil sweep")
    showcase_movie()
    print(f"[movie] MP4 {(OUT / 'wind_tunnel.mp4').stat().st_size / 1e6:.1f} MB")
    stamp("movie + hero (total runtime)")


if __name__ == "__main__":
    main()