Case 06: Stratified isolated building (Zhou et al. 2021)

STATUS: SCAFFOLD - unrun. This notebook is analysis structure only. It runs before the GPU result and the digitized reference data exist (every loader is guarded), and it makes NO physics claim. Reference-data digitization + GPU validation are pending (epic #1034 P5, issue #1039).

Validation of the stratified isolated-building case for validation/wind_engineering/06_stratified_building/06_stratified_building.nassu.yaml (stratifiedBuildingZhouC2RiN01, Case 2, Ri = -0.1). An isolated rectangular building in an unstable ABL, with a ground hole source behind the leeward wall and a buoyant temperature scalar driving the Monin-Obukhov ground wall model.

Reference: Zhou, X. et al. (2021), J. Wind Eng. Ind. Aerodyn. 211, 104526 (TPU Database 2006 wind tunnel for Cases 1 and 2).

Scaling and normalization

Building height H = 0.16 m, footprint 0.5H x 0.5H, u_H = 1.37 m/s, Re = 14811, Ri = -0.1 (Case 2). In lattice units (level 0) H_lat = 32, u_H_lat = 0.05 (Ma = 0.087). Reported quantities:

• velocity / u_H,

• temperature phi = (T - T_r) / dT, dT = 34 K,

• concentration / c0 with c0 = 0.05 q / (u_H H^2).

Only Case 1 (Ri=0) and Case 2 (Ri=-0.1) are wind-tunnel-backed; Cases 3-6 are CFD-only (Zhou’s own LES). Label every reference series with its case and never present a CFD-only case as measurement data.

import pathlib

import matplotlib.pyplot as plt
import pandas as pd


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/06_stratified_building"
results_dir = project_root / "output" / case
reference_dir = project_root / case / "reference"

# Scaling constants (see the case YAML header).
H_LAT = 32.0
U_H_LAT = 0.05
PLANE_HEIGHT = 5.01
DT_K = 34.0
print("results dir:", results_dir)

Load model outputs (guarded)

The case records the inlet line, the behind-building lines a,b,c, the probe points A-E, the symmetry plane and the full-domain statistics. Wire these readers to the actual export format once the GPU run exists; until then they return None and the plots degrade gracefully.

def load_series(export_name: str, entity: str):
    """Return a model series or None if absent."""
    if not results_dir.exists():
        return None
    # TODO(P5): point this at the real line/point/plane export and
    # time-average over the statistics window.
    hits = list(results_dir.rglob(f"*{export_name}*{entity}*"))
    if not hits:
        return None
    raise NotImplementedError("Wire the series reader to the run output.")


model_lines = {ln: load_series("behind_building_lines", f"line_{ln}") for ln in ("a", "b", "c")}
model_inlet = load_series("inlet_profile", "inlet")
have_model = any(v is not None for v in model_lines.values())
print("model data present:", have_model)

Load digitized reference data (guarded)

Reference CSVs live under reference/. They are EMPTY PLACEHOLDERS until digitized from the Zhou et al. figures / TPU Database 2006 at P5 (see reference/README.md). Each carries a source column naming its figure and case; that source is used verbatim in every plot legend.

def load_reference(name: str):
    path = reference_dir / name
    if not path.exists():
        return None
    df = pd.read_csv(path, comment="#")
    return df if not df.empty else None


def source_label(df, fallback="Zhou et al. (2021)"):
    if df is not None and "source" in df and len(df):
        return str(df["source"].iloc[0])
    return fallback


ref_inlet = load_reference("inlet_profiles.csv")
ref_lines = load_reference("behind_building_lines.csv")
ref_recirc = load_reference("recirc_length_vs_Ri.csv")
have_ref = ref_lines is not None
print("reference data present:", have_ref)

Behind-building lines a, b, c (Fig. 7)

Streamwise velocity and concentration at x/H = 0.375, 0.625, 1.0 behind the leeward wall vs height, model against the digitized Zhou et al.

2021. Case 2 data with its reported +/-15% band. Plots only.

fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
for ax, ln in zip(axes, ("a", "b", "c")):
    if ref_lines is not None and "line" in ref_lines:
        sub = ref_lines[ref_lines["line"] == ln]
        if len(sub):
            err = sub["err_frac"] * sub["c_over_c0"] if "err_frac" in sub else None
            ax.errorbar(
                sub["c_over_c0"],
                sub["z_over_H"],
                xerr=err,
                fmt="ks",
                label=source_label(ref_lines),
            )
    if model_lines.get(ln) is not None:
        pass  # TODO(P5): overlay model <c>/c0 vs z/H for this line.
    ax.set_title(f"line {ln}")
    ax.set_xlabel("<c> / c0")
    if ref_lines is None and not have_model:
        ax.text(0.5, 0.5, "pending", ha="center", va="center", transform=ax.transAxes)
    if ax.get_legend_handles_labels()[0]:
        ax.legend()
axes[0].set_ylabel("z / H")
fig.suptitle("Behind-building concentration (Zhou et al. 2021, Fig. 7) - SCAFFOLD")
fig.tight_layout()
plt.show()

Inlet profiles (Fig. 5)

Inlet u, k and temperature profiles vs the TPU Database 2006. Tests the provisional SEM + mean-T inflow (P5 fits the inlet to the TPU profiles).

fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
for ax, col, title in zip(
    axes, ["U_over_uH", "k_over_uH2", "phi_T"], ["u / u_H", "k / u_H^2", "phi_T"]
):
    if ref_inlet is not None and col in ref_inlet:
        ax.plot(
            ref_inlet[col],
            ref_inlet["z_over_H"],
            "k^",
            label=source_label(ref_inlet, "TPU Database 2006 (Zhou et al. 2021)"),
        )
    if model_inlet is not None:
        pass  # TODO(P5): overlay model inlet profile.
    ax.set_title(title)
    ax.set_xlabel(title)
    if ref_inlet is None and model_inlet is None:
        ax.text(0.5, 0.5, "pending", ha="center", va="center", transform=ax.transAxes)
    if ax.get_legend_handles_labels()[0]:
        ax.legend()
axes[0].set_ylabel("z / H")
fig.suptitle("Inlet profiles vs TPU Database 2006 (Zhou et al. 2021, Fig. 5) - SCAFFOLD")
fig.tight_layout()
plt.show()

Recirculation length vs Ri (Fig. 9)

Normalized recirculation length behind the building vs Richardson number (Zhou reports it shrinking ~32% from Ri 0 to -1.5). Model points come from the Ri ladder (P5); the reference is the digitized Zhou series, with wind-tunnel-backed points (Ri 0, -0.1) marked distinctly from the CFD-only points.

fig, ax = plt.subplots(figsize=(7, 4))
if ref_recirc is not None and "Ri" in ref_recirc:
    wt = ref_recirc[ref_recirc.get("backing") == "wind_tunnel"]
    cfd = ref_recirc[ref_recirc.get("backing") == "CFD_only"]
    if len(wt):
        ax.plot(
            wt["Ri"],
            wt["Lr_over_H"],
            "ks",
            label=source_label(ref_recirc, "Zhou et al. 2021 (wind tunnel)"),
        )
    if len(cfd):
        ax.plot(cfd["Ri"], cfd["Lr_over_H"], "kx", label="Zhou et al. 2021 (LES, CFD-only)")
# TODO(P5): overlay Nassu ladder points (Ri vs recirc length).
ax.set_xlabel("Ri")
ax.set_ylabel("L_r / H")
ax.set_title("Recirculation length vs Ri (Zhou et al. 2021, Fig. 9) - SCAFFOLD")
if ref_recirc is None:
    ax.text(
        0.5,
        0.5,
        "pending GPU ladder + reference digitization",
        ha="center",
        va="center",
        transform=ax.transAxes,
    )
if ax.get_legend_handles_labels()[0]:
    ax.legend()
fig.tight_layout()
plt.show()

Fields + PSDs (Figs. 8, 10-13)

Symmetry-plane streamlines/contours (Fig. 8), vertical velocity (Fig. 10), u’w’/u_H^2 (Fig. 11) and the point PSDs at A-E (Figs. 12-13). Wire the plane and point-series readers at P5 and render contourf / spectrum plots (not tables).

# TODO(P5): load 'symmetry_plane/vertical_symmetry' -> contourf u/T/c;
# load 'probe_points' A-E -> Welch PSD vs Strouhal, one line per point,
# each legend-labelled with its measurement point and the Zhou source.
print("Field + PSD plots: pending GPU result.")