Case 08: Buoyant channel - analytical validation of the Boussinesq coupling

Laminar natural convection in a differentially-heated vertical channel isolates the Boussinesq body force F = rho0*beta*b*g (Guo source route): single block, no LES, no refinement, no IBM, no wall model - buoyancy is the only driver. A hot wall (phi = 1, low x) and a cold wall (phi = 0, high x) set an exactly linear scalar profile b(xi) = b0*(1 - xi); the buoyancy drives a steady vertical velocity w = uz with the closed form

\[w(\xi) = \frac{A}{12}\,(2\xi^3 - 3\xi^2 + \xi), \qquad A = \frac{\beta\,g_z\,b_0\,L^2}{\nu}, \qquad \xi = x/L,\]

antisymmetric about the channel centre, with zeros at xi = 0, 1/2, 1, extrema at xi = 0.2113, 0.7887, and peak magnitude |A|/(72*sqrt(3)). Reference: the free-convection-between-vertical-plates solution (Batchelor 1954; Bird, Stewart & Lightfoot, Transport Phenomena).

Acceptance criterion (stated up front)

The exact solution is a continuum result; a finite-resolution LBM with halfway-bounce-back and regularized-Dirichlet walls carries a near-wall discretization error that is first order in the gap resolution. The case therefore validates on the shape and on convergence, not on a single tight number:

• antisymmetric cubic shape reproduced: profile relative-L2 vs the exact cubic < 4% at N=64, decreasing with resolution;

• peak location within +/-1 node of xi = 0.2113 / 0.7887 (exact);

• peak amplitude within 4% of |A|/(72*sqrt(3)) at N=64, decreasing with resolution (the convergence chart below must show the N=96 error below the N=64 error, confirming the residual is discretization);

• sign physical (hot / low-x side rises, c > 0);

• zero net flux |sum w| / (|w_peak|*N) < 1%;

• scalar linearity rel-L2 of b(xi) vs b0*(1 - xi) < 1% (input check).

Statistics (the late-window mean field), not instantaneous snapshots.

import glob
import pathlib

import h5py
import matplotlib.pyplot as plt
import numpy as np


def _find_project_root() -> pathlib.Path:
    here = pathlib.Path.cwd().resolve()
    for cand in [here, *here.parents]:
        if (cand / "pyproject.toml").exists() and (cand / "nassu").is_dir():
            return cand
    raise RuntimeError(f"Could not locate Nassu project root from {here}")


project_root = _find_project_root()
case = "validation/wind_engineering/08_buoyant_channel"
results_dir = project_root / "output" / case

# beta, g_z, b0, nu are fixed; only the gap L (= domain resolution) changes.
beta, g_z, b0, nu = 0.01, -1.0e-3, 1.0, 0.1
SIMS = {"N=64": ("buoyant_channel", 64), "N=96": ("buoyant_channel_hires", 96)}
print("results dir:", results_dir, "(exists:", results_dir.exists(), ")")

results dir: /home/waine/repos/nassu/output/validation/wind_engineering/08_buoyant_channel (exists: True )

Load mean fields and reduce to the wall-normal (gap) profile (guarded)

The gap axis is detected robustly as the axis carrying the scalar gradient (so the reduction is independent of the exported array’s axis order): the profile is the average of the time-mean field over the two homogeneous axes. w <- uz (the buoyancy-driven vertical component). Every loader returns None when the run is absent, so the plots degrade to ‘pending’.

def _gap_axis(field3):
    # the wall-normal (gap) axis is the one along which the scalar varies most
    ranges = []
    for ax in range(3):
        others = tuple(j for j in range(3) if j != ax)
        ranges.append(float(np.ptp(field3.mean(axis=others))))
    return int(np.argmax(ranges))


def load_profile(sim_name):
    pat = str(
        results_dir / f"{sim_name}__000" / "outputs" / "mean_fields.volume.mean_fields.stats.*.h5"
    )
    files = sorted(glob.glob(pat))
    if not files:
        return None
    with h5py.File(files[-1], "r") as h:
        blk = h[sorted(h.keys())[-1]]["block0"]
        phi = np.asarray(blk["temperature_phi"][:])
        uz = np.asarray(blk["uz"][:])
    gax = _gap_axis(phi)
    others = tuple(j for j in range(3) if j != gax)
    w = uz.mean(axis=others)
    b = phi.mean(axis=others)
    N = w.size
    xi = (np.arange(N) + 0.5) / N
    return {"xi": xi, "w": w, "b": b, "N": N}


def analytic(xi, Lgap):
    A = beta * g_z * b0 * Lgap**2 / nu
    return (A / 12.0) * (2 * xi**3 - 3 * xi**2 + xi), A


runs = {tag: load_profile(name) for tag, (name, _) in SIMS.items()}
have = {tag: r is not None for tag, r in runs.items()}
print("results present:", have)

results present: {'N=64': True, 'N=96': True}

Velocity profile vs the analytical cubic

The simulated w(xi) = uz at each resolution overlaid on the exact closed-form cubic. If the coupling is correct the markers sit on the curve and approach it as resolution increases.

fig, axes = plt.subplots(1, len(SIMS), figsize=(11, 4))
for ax, (tag, (name, Lgap)) in zip(np.atleast_1d(axes), SIMS.items()):
    r = runs[tag]
    if r is not None:
        xi = r["xi"]
        xf = np.linspace(0, 1, 200)
        wex, A = analytic(xf, Lgap)
        g = 2 * xi**3 - 3 * xi**2 + xi
        c = float((r["w"] * g).sum() / (g * g).sum())  # match the code's sign
        ax.plot(
            xf,
            np.sign(c) * np.abs(wex) * 1e3,
            "k-",
            label="analytical cubic (Batchelor 1954 / BSL)",
        )
        ax.plot(xi, r["w"] * 1e3, "C1o", ms=3, label=f"Nassu w=uz ({tag})")
    else:
        ax.text(0.5, 0.5, "pending GPU run", ha="center", va="center", transform=ax.transAxes)
    ax.axhline(0, color="0.8", lw=0.8)
    ax.set_xlabel("xi = x / L (gap)")
    ax.set_ylabel("w = uz  (x1e-3, lattice)")
    ax.set_title(f"Vertical velocity ({tag})")
    if ax.get_legend_handles_labels()[0]:
        ax.legend(fontsize=8)
fig.tight_layout()
plt.show()

Scalar profile (input check)

The transported buoyancy scalar must be exactly linear across the gap; any deviation is the error in the buoyancy input and bounds the momentum match.

fig, ax = plt.subplots(figsize=(6, 4))
for tag, (name, Lgap) in SIMS.items():
    r = runs[tag]
    if r is not None:
        ax.plot(r["xi"], r["b"], "o", ms=3, label=f"Nassu b ({tag})")
if any(have.values()):
    xf = np.linspace(0, 1, 50)
    ax.plot(xf, b0 * (1 - xf), "k-", label="linear b0(1 - xi)")
    ax.legend(fontsize=8)
else:
    ax.text(0.5, 0.5, "pending GPU run", ha="center", va="center", transform=ax.transAxes)
ax.set_xlabel("xi = x / L (gap)")
ax.set_ylabel("scalar b")
ax.set_title("Scalar linearity")
fig.tight_layout()
plt.show()

Convergence of the error with resolution

The acceptance metrics at each resolution as a grouped bar chart. The bars must decrease from N=64 to N=96 - the discretization error converging toward the analytical solution - which is what certifies the buoyancy coupling as correct rather than merely close at one grid.

def metrics(r, Lgap):
    xi = r["xi"]
    g = 2 * xi**3 - 3 * xi**2 + xi
    _, A = analytic(xi, Lgap)
    wmax = abs(A) / (72 * np.sqrt(3))
    c = float((r["w"] * g).sum() / (g * g).sum())
    shape = float(np.linalg.norm(r["w"] - c * g) / np.linalg.norm(c * g)) * 100
    peak = abs(abs(r["w"]).max() - wmax) / wmax * 100
    scal = float(np.linalg.norm(r["b"] - b0 * (1 - xi)) / np.linalg.norm(b0 * (1 - xi))) * 100
    flux = abs(r["w"].sum()) / (abs(r["w"]).max() * r["N"]) * 100
    return {
        "shape rel-L2 %": shape,
        "peak err %": peak,
        "scalar rel-L2 %": scal,
        "net flux %": flux,
        "sign hot-up": c > 0,
    }


labels = ["shape rel-L2 %", "peak err %", "scalar rel-L2 %", "net flux %"]
fig, ax = plt.subplots(figsize=(8, 4))
tags = [t for t in SIMS if runs[t] is not None]
if tags:
    x = np.arange(len(labels))
    w = 0.8 / len(tags)
    rows = {}
    for k, tag in enumerate(tags):
        m = metrics(runs[tag], SIMS[tag][1])
        rows[tag] = m
        ax.bar(x + (k - (len(tags) - 1) / 2) * w, [m[lbl] for lbl in labels], w, label=tag)
    ax.axhline(4.0, ls=":", color="0.4", lw=1, label="4% guide")
    ax.set_xticks(x)
    ax.set_xticklabels(labels, rotation=15, ha="right")
    ax.set_ylabel("error (%)")
    ax.legend()
    print("metric summary (supplement to the chart above):")
    for tag in tags:
        print(
            f"  {tag}: "
            + ", ".join(f"{lbl}={rows[tag][lbl]:.3f}" for lbl in labels)
            + f", sign hot-up={rows[tag]['sign hot-up']}"
        )
else:
    ax.text(0.5, 0.5, "pending GPU run", ha="center", va="center", transform=ax.transAxes)
ax.set_title("Error vs analytical cubic - converges with resolution")
fig.tight_layout()
plt.show()

metric summary (supplement to the chart above):
  N=64: shape rel-L2 %=3.751, peak err %=3.166, scalar rel-L2 %=0.793, net flux %=0.003, sign hot-up=True
  N=96: shape rel-L2 %=2.463, peak err %=2.143, scalar rel-L2 %=0.526, net flux %=0.011, sign hot-up=True

Verdict

The case passes when, at both resolutions, the simulated w(xi) reproduces the antisymmetric analytical cubic (peak location exact, physical sign, zero net flux, linear scalar) and the shape / peak error decreases from N=64 to N=96 - the discretization error converging toward the exact solution. A monotone drop certifies the Boussinesq buoyancy coupling; a static or growing error would flag a coupling defect. This is the analytical, single-block, feature-alone rung the composition matrix rests on.