Case 07: Stratified-ABL neutral-limit regression

STATUS: SCAFFOLD - unrun (P3). This notebook is analysis structure only. It runs before the GPU result exists (every loader is guarded and degrades to a ‘pending’ placeholder), and it makes NO physics claim. This case is a SELF-REGRESSION: the baseline run IS the reference. GPU validation is pending (epic #1034 P5, issue #1039).

Mandatory neutral-limit regression for the stratified-ABL feature (Monin-Obukhov IBM ground wall model + Boussinesq buoyancy). With the feature OFF (baseline) or at its INERT parameter (neutral_limit: surface_heat_flux = 0, beta = 0), the solver must recover its pre-feature behaviour to the solver noise floor. The same pair is run across a 2:1 refinement slab (baseline_mb, neutral_limit_mb) to prove the inert buoyancy / temperature field stays continuous across the interface.

Config: validation/wind_engineering/07_stratified_neutral_limit/07_stratified_neutral_limit.nassu.yaml.

Acceptance criterion (stated up front)

The GPU solver is NOT bit-reproducible run-to-run (atomicAdd ordering, ~1e-6 relative in float32), so we validate STATISTICS, not instantaneous fields, and prove parity against the solver’s own noise floor. Using the time-averaged mean fields over [10000, 20000] steps, x-y-averaged into z-profiles of u, rho and temperature:

• Delta_feature = || neutral_limit(z) - baseline(z) ||

• Delta_noise = || baseline_rerun(z) - baseline(z) || (solver noise)

The regression PASSES iff Delta_feature <~ Delta_noise (same order, ~1e-6 relative) for every quantity. For the multiblock variant, the Delta_feature bars must stay at the same solver-noise floor as the single-block case, showing the inert feature composes across the 2:1 refinement interface (F2C/C2F/border reconstruction) with no discontinuity.

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/07_stratified_neutral_limit"
results_dir = project_root / "output" / case

# The regression sims. The `_rerun` sims establish the solver noise floor.
SINGLE = ["baseline", "baseline_rerun", "neutral_limit"]
MULTIBLOCK = ["baseline_mb", "baseline_mb_rerun", "neutral_limit_mb"]
PROFILE_MACRS = ["u", "rho", "temperature_phi"]
print("results dir:", results_dir, "(exists:", results_dir.exists(), ")")

results dir: /home/waine/repos/nassu/output/validation/wind_engineering/07_stratified_neutral_limit (exists: True )

Load mean fields and reduce to z-profiles (guarded)

Each sim exports first-order stats of rho, u, temperature_phi (data.exports.mean_fields). x and y are periodic (homogeneous), so the z-profile is the x-y average of the time-mean field. Every loader returns None when the run is absent, so the plots below degrade to a ‘pending’ placeholder without error.

def load_profile(sim_name: str):
    """Return {macr: z_profile_1d} for a sim, or None if the run is absent.

    z-profile = x-y average of the time-mean field over the stats window.
    """
    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
    # The accumulated 1st-order mean field (level 0), one time group.
    # Axis order is Fortran (0=x, 1=y, 2=z); x, y are periodic homogeneous
    # directions, so the z-profile is the x-y average. 'u' <- ux (the
    # streamwise, forced component).
    src = {"u": "ux", "rho": "rho", "temperature_phi": "temperature_phi"}
    with h5py.File(files[-1], "r") as h:
        blk = h[sorted(h.keys())[-1]]["block0"]
        return {m: np.asarray(blk[src[m]][:]).mean(axis=(0, 1)) for m in PROFILE_MACRS}


def load_group(names):
    return {n: load_profile(n) for n in names}


def have(group):
    return group and all(v is not None for v in group.values())


single = load_group(SINGLE)
multiblock = load_group(MULTIBLOCK)
print("single-block results present:", have(single))
print("multiblock results present:", have(multiblock))

single-block results present: True
multiblock results present: True

Profile overlay: neutral_limit vs baseline (single block)

The feature-inert run (neutral_limit) overlaid on the feature-off reference (baseline) and the noise-floor rerun (baseline_rerun) for u(z), rho(z) and temperature(z). If the feature is truly inert the three curves are indistinguishable. Plots only.

fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
titles = {"u": "u(z)", "rho": "rho(z)", "temperature_phi": "temperature(z)"}
styles = {
    "baseline": ("k-", "baseline (feature OFF)"),
    "baseline_rerun": ("0.6", "baseline rerun (noise floor)"),
    "neutral_limit": ("C1--", "neutral_limit (feature inert)"),
}
for ax, macr in zip(axes, PROFILE_MACRS):
    if have(single):
        for sim, (fmt, lbl) in styles.items():
            prof = single[sim][macr]
            z = np.arange(len(prof))
            ax.plot(prof, z, fmt, label=lbl)
    else:
        ax.text(0.5, 0.5, "pending GPU run", ha="center", va="center", transform=ax.transAxes)
    ax.set_title(titles[macr])
    ax.set_xlabel(titles[macr])
    if ax.get_legend_handles_labels()[0]:
        ax.legend(fontsize=8)
axes[0].set_ylabel("z (lattice)")
fig.suptitle("Neutral-limit regression: profile overlay (single block) - SCAFFOLD")
fig.tight_layout()
plt.show()

Parity bar chart: Delta_feature vs Delta_noise

For each quantity, the profile L2 difference of the feature-inert run against baseline (Delta_feature) next to the baseline-vs-rerun solver noise floor (Delta_noise). The regression PASSES when the orange bar does not exceed the grey bar (same order). A grouped bar chart - a scalar decomposition shown as bars, never a bare table.

def profile_l2(a, b):
    a, b = np.asarray(a), np.asarray(b)
    denom = np.linalg.norm(b) or 1.0
    return float(np.linalg.norm(a - b) / denom)  # relative L2


def parity_bars(group, label):
    fig, ax = plt.subplots(figsize=(7, 4))
    x = np.arange(len(PROFILE_MACRS))
    w = 0.38
    if have(group):
        base = group[[k for k in group if k.startswith("baseline") and "rerun" not in k][0]]
        rerun = group[[k for k in group if "rerun" in k][0]]
        feat = group[[k for k in group if k.startswith("neutral")][0]]
        d_noise = [profile_l2(rerun[m], base[m]) for m in PROFILE_MACRS]
        d_feat = [profile_l2(feat[m], base[m]) for m in PROFILE_MACRS]
        ax.bar(x - w / 2, d_noise, w, label="Delta_noise (baseline rerun)", color="0.6")
        ax.bar(x + w / 2, d_feat, w, label="Delta_feature (neutral_limit)", color="C1")
        ax.set_yscale("log")
        ax.axhline(max(d_noise), ls=":", color="0.4", lw=1)
        ax.legend()
    else:
        ax.text(0.5, 0.5, "pending GPU run", ha="center", va="center", transform=ax.transAxes)
    ax.set_xticks(x)
    ax.set_xticklabels(["u", "rho", "temperature"])
    ax.set_ylabel("relative L2 profile difference")
    ax.set_title(f"Parity vs solver noise floor ({label}) - SCAFFOLD")
    fig.tight_layout()
    plt.show()


parity_bars(single, "single block")

Multiblock variant: parity across the 2:1 refinement

The same parity check with a 2:1 refinement slab present. If the Delta_feature bars sit at the same solver-noise floor as the single-block case, the inert feature composes across the refinement interface (F2C/C2F/border reconstruction) without introducing any discontinuity - a term applied in bulk collision but omitted from border reconstruction or the level-transfer rescale would push Delta_feature above the floor here. The multiblock parity is the robust interface-composition check (the level-0 stitched export layout makes a raw per-node interface profile unreliable, so parity-vs-noise is used rather than a hand-placed interface line).

parity_bars(multiblock, "multiblock 2:1")

Verdict (populated at P5)

Once the run exists, the regression passes iff Delta_feature <~ Delta_noise for every quantity in BOTH the single-block and multiblock parity charts, and the multiblock temperature profile shows no sawtooth at the refinement interfaces. Any failure means the (supposedly inert) feature code path has moved the pre-feature physics and must be fixed before the 05 / 06 benchmark rungs are trusted.