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 numpy as np
import matplotlib.pyplot as plt


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/Documents/Codigos/AeroSim/nassu-stratified-abl/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.