Case 17: CODASC isolated street canyon - pollutant line source

Validation of the flagship pollutant-dispersion case for validation/wind_engineering/04_codasc_canyon/04_codasc_canyon.nassu.yaml (codascCanyonEmpty). An isolated, tree-free street canyon (W/H = 1, L/H = 10) under a neutral ABL perpendicular to the canyon, with a continuous ground-level line source of a passive tracer on the canyon floor. We compare the dimensionless wall concentration K = c+ along the leeward (wall A) and windward (wall B) building facades against the CODASC wind-tunnel database (KIT Karlsruhe), and report the COST 732 / VDI 3783 Part 9 hit-rate metrics.

This notebook is authored to run after a GPU result and the CODASC reference data exist. It skips gracefully (and the asserts are guarded) when either is missing, so it can live in the repo before the heavy run.

References: Gromke (2011) Environ. Pollut. 159:2094-2099; Gromke & Ruck (2012) Boundary-Layer Meteorol. 144:41-64; CODASC database (www.codasc.de); COST Action 732 (Schatzmann et al. 2010); VDI 3783 Part 9 (2005).

Scaling and the K normalization

Model scale 1:150, building height H = 0.12 m, roof velocity U_H = 4.65 m/s, Re_H = 37000. In lattice units (level 0) H_lat = 24, U_H_lat = 0.05 (Ma = 0.087). The dimensionless wall concentration is

\[K = c^{+} = \frac{c\, U_H\, H\, l_s}{Q_{SF6}}.\]

K is invariant to a constant rescaling of the source strength provided the same rate is used in K. The case injects a localized floor source of peak strength S0 = 1e-2 chosen for signal-to-noise; with the consistent lattice rate Q_lat = S0 * l_s_lat this gives K = c * U_H_lat * H_lat / S0 = 120 * c.

[ ]:

import pathlib

import matplotlib.pyplot as plt
import numpy as np
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()
results_dir = project_root / "validation/wind_engineering/04_codasc_canyon/results"
comparison_dir = project_root / "validation/wind_engineering/04_codasc_canyon/reference"

# Scaling constants (see the case YAML header).
H_LAT = 24.0
U_H_LAT = 0.05
S0 = 1.0e-2
PLANE_HEIGHT = 5.01
K_FACTOR = U_H_LAT * H_LAT / S0  # K = K_FACTOR * c   (= 120)
print("K_FACTOR =", K_FACTOR)

Load the time-averaged facade concentration profiles

The case records two vertical probe lines, one node off each facade, along the canyon mid span: facade_leeward_A/wall_A_profile and facade_windward_B/wall_B_profile. We time-average pollutant_phi over the statistics window and convert to K.

The loader below is a thin placeholder: adapt the read to the actual probe export format once the GPU run exists (HDF5 / XDMF / CSV series).

[ ]:

def load_facade_K(series_name: str, line_name: str):
    """Load a facade probe line, time-average phi, return (z_over_H, K).

    Returns (None, None) if the result is not present yet so the
    notebook can run before the GPU result exists.
    """
    if not results_dir.exists():
        return None, None
    # TODO(maintainer): point this at the real probe export for the run.
    # The historic-series probe writes one record per sampled step; we
    # average over the statistics window and map node height -> z/H.
    candidates = list(results_dir.rglob(f"*{series_name}*{line_name}*"))
    if not candidates:
        return None, None
    raise NotImplementedError("Wire the probe-series reader to the run output format here.")


z_lee, K_lee = load_facade_K("facade_leeward_A", "wall_A_profile")
z_wind, K_wind = load_facade_K("facade_windward_B", "wall_B_profile")
have_model = K_lee is not None and K_wind is not None
print("model facade data present:", have_model)

Load the CODASC reference profiles

The reference K(z/H) profiles live under validation/wind_engineering/04_codasc_canyon/reference/. They are placeholders until the maintainer populates them from the confirmed www.codasc.de download for the empty (tree-free) W/H = 1, alpha = 90 deg case.

[ ]:

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


ref_lee = load_reference("facade_K_leeward.csv")
ref_wind = load_reference("facade_K_windward.csv")
have_ref = ref_lee is not None and ref_wind is not None
print("CODASC reference data present:", have_ref)

COST 732 / VDI 3783 Part 9 metrics

Computed over the paired (model, observed) facade K samples, after interpolating the model profile onto the reference heights.

  • FAC2: fraction with 0.5 <= K_mod/K_obs <= 2 (target >= 0.5).

  • NMSE: mean((Ko-Km)^2)/(mean(Ko) mean(Km)) (target < 4).

  • FB: (mean(Ko)-mean(Km))/(0.5(mean(Ko)+mean(Km))) (in [-0.3, 0.3]).

  • R: Pearson correlation.

  • hit-rate q: fraction within |Km-Ko| <= W*|Ko| or <= D, with W = 0.25, D = 0.25 (target >= 0.66).

[ ]:

def cost732_metrics(K_obs, K_mod, W=0.25, D=0.25):
    K_obs = np.asarray(K_obs, dtype=float)
    K_mod = np.asarray(K_mod, dtype=float)
    mask = np.isfinite(K_obs) & np.isfinite(K_mod)
    K_obs, K_mod = K_obs[mask], K_mod[mask]
    n = K_obs.size
    if n == 0:
        return {"n": 0}
    ratio = np.divide(K_mod, K_obs, out=np.full_like(K_mod, np.inf), where=K_obs != 0)
    fac2 = np.mean((ratio >= 0.5) & (ratio <= 2.0))
    mo, mm = K_obs.mean(), K_mod.mean()
    nmse = np.mean((K_obs - K_mod) ** 2) / (mo * mm) if mo * mm > 0 else np.inf
    fb = (mo - mm) / (0.5 * (mo + mm)) if (mo + mm) != 0 else np.inf
    r = (
        np.corrcoef(K_obs, K_mod)[0, 1]
        if n > 1 and K_obs.std() > 0 and K_mod.std() > 0
        else np.nan
    )
    hit = np.mean(np.abs(K_mod - K_obs) <= np.maximum(W * np.abs(K_obs), D))
    return {"n": n, "FAC2": fac2, "NMSE": nmse, "FB": fb, "R": r, "q": hit}


def paired(z_ref, K_ref, z_mod, K_mod):
    """Interpolate the model K onto the reference heights."""
    order = np.argsort(z_mod)
    return np.interp(z_ref, np.asarray(z_mod)[order], np.asarray(K_mod)[order])

[ ]:

metrics = {}
if have_model and have_ref:
    for tag, z_m, K_m, ref in [
        ("leeward", z_lee, K_lee, ref_lee),
        ("windward", z_wind, K_wind, ref_wind),
    ]:
        K_m_on_ref = paired(ref["z_over_H"], ref["K"], z_m, K_m)
        metrics[tag] = cost732_metrics(ref["K"], K_m_on_ref)
    pooled_obs = np.concatenate([ref_lee["K"], ref_wind["K"]])
    pooled_mod = np.concatenate(
        [
            paired(ref_lee["z_over_H"], ref_lee["K"], z_lee, K_lee),
            paired(ref_wind["z_over_H"], ref_wind["K"], z_wind, K_wind),
        ]
    )
    metrics["pooled"] = cost732_metrics(pooled_obs, pooled_mod)
    for tag, m in metrics.items():
        print(tag, {k: round(v, 4) for k, v in m.items()})
else:
    print("Skipping metrics: model and/or reference data not present yet.")

Acceptance asserts

Per the issue acceptance criteria: FAC2 >= 0.5, NMSE < 4, hit-rate q >= 0.66 on the pooled leeward+windward facade samples. Guarded so the notebook is a no-op until both data sources exist.

[ ]:

if metrics.get("pooled", {}).get("n", 0) > 0:
    m = metrics["pooled"]
    assert m["FAC2"] >= 0.5, f"FAC2 {m['FAC2']:.3f} < 0.5"
    assert m["NMSE"] < 4.0, f"NMSE {m['NMSE']:.3f} >= 4"
    assert m["q"] >= 0.66, f"hit-rate q {m['q']:.3f} < 0.66"
    print("CODASC acceptance metrics PASSED")
else:
    print("Asserts skipped (no paired data yet).")

Facade K profiles

Leeward (wall A) and windward (wall B) dimensionless wall concentration vs height, model against CODASC.

[ ]:

fig, axes = plt.subplots(1, 2, figsize=(10, 5), sharey=True)
for ax, tag, z_m, K_m, ref in [
    (axes[0], "leeward (wall A)", z_lee, K_lee, ref_lee),
    (axes[1], "windward (wall B)", z_wind, K_wind, ref_wind),
]:
    if z_m is not None and K_m is not None:
        z_over_H = (np.asarray(z_m) - PLANE_HEIGHT) / H_LAT
        ax.plot(np.asarray(K_m), z_over_H, "-", label="Nassu")
    if ref is not None:
        ax.plot(ref["K"], ref["z_over_H"], "ks", label="CODASC")
    ax.set_title(tag)
    ax.set_xlabel("K = c+")
    ax.legend()
axes[0].set_ylabel("z / H")
fig.tight_layout()
plt.show()

Canyon concentration field

Time-averaged pollutant field on the vertical streamwise plane through the canyon mid span (plane_series/canyon_mid_span), showing the in-canyon recirculation transporting the floor release toward the leeward wall. Wire the plane reader to the run output once it exists.

[ ]:

# TODO(maintainer): load the `canyon_mid_span` plane export and
# imshow / contourf K = K_FACTOR * pollutant_phi here once results exist.
print("Canyon field plot: pending GPU result.")

Version

[ ]:

from nassu.cfg.model import ConfigScheme

case_path = (
    project_root / "validation/wind_engineering/04_codasc_canyon/04_codasc_canyon.nassu.yaml"
)
sim_cfg = next(iter(ConfigScheme.sim_cfgs_from_file_dct(str(case_path)).values()))
sim_info = sim_cfg.output.read_info()

nassu_commit = sim_info["commit"]
nassu_version = sim_info["version"]
print("Version:", nassu_version)
print("Commit hash:", nassu_commit)

Configuration

[ ]:

from IPython.display import Code

Code(filename=str(case_path))