Case 16 - CEDVAL A1-5 isolated building pollutant dispersion

Validation notebook for validation/wind_engineering/03_cedval_building/03_cedval_building.nassu.yaml (GitHub issue #652). An isolated rectangular building sits in a neutral atmospheric boundary layer (ABL); a near-ground tracer source on the leeward wall base emits a passive scalar pollutant. We compare the simulated dimensionless surface/near-wake concentration against the CEDVAL A1-5 wind-tunnel dataset (Univ. Hamburg EWTL) using the COST 732 / VDI 3783-9 CWE validation metrics.

This notebook is authored to run AFTER (a) a GPU run has produced the time-averaged result and (b) the maintainer has populated the reference data ``validation/wind_engineering/03_cedval_building/reference/c_star_points.csv`` from the password-protected EWTL dataset. It is intentionally NOT executed in CI.

What is validated

  • Scalar transport (advection-diffusion of the pollutant DDF).

  • The voxel surface source (on-body / near-ground emission), isolated cleanly in this simple geometry.

  • Impermeable building faces (voxel zero-flux scalar BC ScalarHWBB).

  • LES + turbulent Schmidt number coupling (Sc_t = 0.7).

Scaling (model -> lattice)

CEDVAL A1-5 building L x W x H = 100 x 150 x 125 mm (model scale 1:200). The height is resolved with H = 50 lattice cells, giving a cell size dx = 125 mm / 50 = 2.5 mm/cell and L = 40, W = 60 lattice cells. The domain is ~15 H x 10 H x 6 H, rounded to 752 x 504 x 304 cells.

The dimensionless concentration follows the CEDVAL convention:

\[c^{*} = \frac{C \, U_{ref} \, H^{2}}{Q}\]

FLAGGED unverified parameters

The following must be confirmed against the EWTL dataset before the absolute comparison is meaningful:

  • ``Q`` (emission rate): the case uses a unit Dirichlet wall value phi_w = 1.0; the effective Q is recovered from the integrated wall flux of the converged field (computed below). The physical Q from the dataset is still needed to place C in absolute units. # VERIFY: CEDVAL A1-5 emission rate Q from the EWTL dataset

  • ``U_ref`` + reference height: the lattice U_ref = 0.06 is chosen for Ma/Re only; the dataset’s physical U_ref and reference height set the c* normalisation. # VERIFY: CEDVAL A1-5 U_ref value + reference height

Reference

CEDVAL A1-5, Univ. Hamburg EWTL (Leitl 1998/2000). Acceptance metrics per COST Action 732 / VDI 3783 Part 9 (Schatzmann et al. 2010; Chang & Hanna 2004): FAC2 >= 0.5, NMSE < 4, FB in [-0.3, 0.3], hit-rate q >= 0.66 (D = 0.25).

[ ]:

import pathlib

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

try:
    import pyvista as pv

    pv.OFF_SCREEN = True
except Exception:  # pragma: no cover - pyvista optional for metric-only runs
    pv = None

[ ]:

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 upward from {here}")


project_root = _find_project_root()
case_path = (
    project_root / "validation/wind_engineering/03_cedval_building/03_cedval_building.nassu.yaml"
)
results_dir = project_root / "validation/wind_engineering/03_cedval_building/results"
comparison_dir = project_root / "validation/wind_engineering/03_cedval_building/reference"
ref_csv = comparison_dir / "c_star_points.csv"
print("case:", case_path)
print("results:", results_dir)
print("reference:", ref_csv)

Scaling constants

H, dx and the lattice reference velocity follow the YAML header. The physical U_ref and Q are FLAGGED placeholders to be confirmed from the dataset; for the dimensionless comparison we use the lattice values and the emission rate recovered from the wall flux.

[ ]:

# Lattice geometry (see the YAML header).
H_LAT = 50  # building height in lattice cells
L_LAT = 40  # along-wind length
W_LAT = 60  # cross-wind width
DX_MM = 2.5  # lattice cell size, mm/cell (model frame)
BODY_X = 250  # building centre x (lattice)
DOMAIN_W = 504  # domain width (lattice)
PLANE_H = 5.01  # ground plane height (lattice)

# Lattice reference velocity (top-of-profile). # VERIFY physical U_ref.
U_REF_LAT = 0.06

# Emission rate Q. The Dirichlet patch uses phi_w = 1.0; the effective
# lattice Q is the integrated scalar wall flux of the converged field,
# computed below. # VERIFY: absolute physical Q from the EWTL dataset.
Q_LAT = None  # set from the integrated wall flux (see the flux cell)

Load the time-averaged result

Loads the exported statistics / instantaneous XDMF for the converged pollutant_phi field. Adjust the glob to the actual exported artefact name once the GPU run exists.

[ ]:

def load_timeavg_field():
    """Load the time-averaged pollutant field from the exported results.

    Returns a pyvista MultiBlock (or dataset). The exact artefact name
    depends on the export config; prefer the statistics mean export, else
    fall back to the last instantaneous snapshot.
    """
    if pv is None:
        raise RuntimeError("pyvista is required to load the field")
    candidates = sorted(results_dir.rglob("*.xdmf")) + sorted(results_dir.rglob("*.pvd"))
    if not candidates:
        raise FileNotFoundError(f"No XDMF/PVD result under {results_dir}. Run the GPU case first.")
    # Prefer a statistics/mean artefact when present.
    stats = [c for c in candidates if "stat" in c.name.lower() or "mean" in c.name.lower()]
    target = stats[-1] if stats else candidates[-1]
    print("loading:", target)
    return pv.read(str(target))


# field = load_timeavg_field()  # uncomment once results exist

Recover the effective emission rate Q from the wall flux

Scalar mass-balance check (acceptance criterion): the scalar emitted at the source equals what is advected out plus what accumulates. For a converged (statistically steady) field the accumulation term vanishes, so the emitted rate equals the net advective+diffusive outflux. We use the source wall flux as the effective lattice Q for the c* normalisation.

[ ]:

def effective_Q_from_field(field):
    """Estimate the effective lattice emission rate Q.

    Integrates the scalar advective flux through a control surface enclosing
    the source (or, equivalently, the net outflux through the outlet plane).
    Implementation depends on the available exported macrs; placeholder
    until the GPU result schema is fixed.
    """
    raise NotImplementedError(
        "Populate from the converged field: integrate the scalar flux "
        "through a control surface around the source / outlet plane."
    )


# Q_LAT = effective_Q_from_field(field)

Sample the predicted c* at the dataset sampling points

Maps the CEDVAL model-frame sampling coordinates (mm, building base origin) to lattice positions via dx and the building placement, samples the simulated pollutant_phi, and normalises to c* = C * U_ref * H^2 / Q.

[ ]:

def model_mm_to_lattice(x_mm, y_mm, z_mm):
    """Map CEDVAL model-frame coords (mm, building base origin) to lattice.

    The building base is centred at (BODY_X, DOMAIN_W/2) on the ground
    plane (z = PLANE_H). Model +x is along-wind, +y cross-wind, +z up.
    """
    x = BODY_X + x_mm / DX_MM
    y = DOMAIN_W / 2.0 + y_mm / DX_MM
    z = PLANE_H + z_mm / DX_MM
    return np.column_stack([x, y, z])


def predicted_c_star(field, ref_df, q_lat, u_ref=U_REF_LAT, h_lat=H_LAT):
    pts = model_mm_to_lattice(
        ref_df["x_model_mm"].to_numpy(),
        ref_df["y_model_mm"].to_numpy(),
        ref_df["z_model_mm"].to_numpy(),
    )
    poly = pv.PolyData(pts)
    sampled = poly.sample(field)
    C = np.asarray(sampled["pollutant_phi"])
    return C * u_ref * h_lat**2 / q_lat

[ ]:

ref_df = pd.read_csv(ref_csv, comment="#")
print(f"reference sampling points: {len(ref_df)}")
if len(ref_df) == 0:
    print(
        "Reference CSV is empty (placeholder). Populate "
        "validation/wind_engineering/03_cedval_building/reference/c_star_points.csv from the "
        "EWTL dataset before running the comparison."
    )
ref_df.head()

COST 732 / VDI 3783-9 validation metrics

FAC2, NMSE, FB, R and the hit-rate q (D = 0.25) over the paired (observed, predicted) c* samples. Definitions follow Chang & Hanna (2004) and VDI 3783 Part 9 (2005).

[ ]:

def cost732_metrics(obs, pred, hitrate_D=0.25):
    obs = np.asarray(obs, dtype=float)
    pred = np.asarray(pred, dtype=float)
    eps = np.finfo(float).tiny
    mo, mp = obs.mean(), pred.mean()
    fb = (mo - mp) / (0.5 * (mo + mp))
    nmse = np.mean((obs - pred) ** 2) / (mo * mp)
    ratio = pred / np.where(np.abs(obs) < eps, eps, obs)
    fac2 = np.mean((ratio >= 0.5) & (ratio <= 2.0))
    r = np.corrcoef(obs, pred)[0, 1] if len(obs) > 1 else np.nan
    # VDI 3783-9 hit-rate: relative tolerance D (absolute fallback W omitted).
    q = np.mean(np.abs(pred - obs) <= hitrate_D * np.abs(obs))
    return {"FAC2": fac2, "NMSE": nmse, "FB": fb, "R": r, "q": q}

[ ]:

# Run once results + reference exist:
#
# field = load_timeavg_field()
# Q_LAT = effective_Q_from_field(field)
# pred = predicted_c_star(field, ref_df, Q_LAT)
# obs = ref_df["c_star_obs"].to_numpy()
# metrics = cost732_metrics(obs, pred)
# print(metrics)
#
# assert metrics["FAC2"] >= 0.5, f"FAC2 {metrics['FAC2']:.3f} < 0.5"
# assert metrics["NMSE"] < 4.0, f"NMSE {metrics['NMSE']:.3f} >= 4"
# assert -0.3 <= metrics["FB"] <= 0.3, f"FB {metrics['FB']:.3f} outside [-0.3, 0.3]"
# assert metrics["q"] >= 0.66, f"hit-rate q {metrics['q']:.3f} < 0.66"
print("Metrics cell ready (commented out until GPU results + dataset exist).")

Plots

  1. Surface / near-ground concentration field on the vertical symmetry plane and a pedestrian-level horizontal plane; (2) a scatter of predicted vs observed c* with the FAC2 band.

[ ]:

def plot_scatter(obs, pred):
    fig, ax = plt.subplots(figsize=(5, 5))
    lo = min(obs.min(), pred.min())
    hi = max(obs.max(), pred.max())
    ax.scatter(obs, pred, s=20, alpha=0.7, label="sampling points")
    ax.plot([lo, hi], [lo, hi], "k-", lw=1, label="1:1")
    ax.plot([lo, hi], [0.5 * lo, 0.5 * hi], "k--", lw=0.8, label="FAC2 band")
    ax.plot([lo, hi], [2.0 * lo, 2.0 * hi], "k--", lw=0.8)
    ax.set_xlabel("observed c*  (CEDVAL A1-5)")
    ax.set_ylabel("predicted c*  (Nassu)")
    ax.set_title("Predicted vs observed dimensionless concentration")
    ax.legend()
    ax.set_aspect("equal", "box")
    return fig


# fig = plot_scatter(obs, pred)  # uncomment once results exist
print("Plot helpers ready.")

Version

[ ]:

from nassu.cfg.model import ConfigScheme

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