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 numpy as np
import pandas as pd
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/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.

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.')