{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7fb27b941602401d91542211134fc71a",
   "metadata": {},
   "source": [
    "# Case 16 - CEDVAL A1-5 isolated building pollutant dispersion\n",
    "\n",
    "Validation notebook for `validation/wind_engineering/03_cedval_building/03_cedval_building.nassu.yaml`\n",
    "(GitHub issue #652). An isolated rectangular building sits in a neutral\n",
    "atmospheric boundary layer (ABL); a near-ground tracer source on the leeward\n",
    "wall base emits a passive scalar `pollutant`. We compare the simulated\n",
    "dimensionless surface/near-wake concentration against the CEDVAL A1-5\n",
    "wind-tunnel dataset (Univ. Hamburg EWTL) using the COST 732 / VDI 3783-9 CWE\n",
    "validation metrics.\n",
    "\n",
    "**This notebook is authored to run AFTER (a) a GPU run has produced the\n",
    "time-averaged result and (b) the maintainer has populated the reference data\n",
    "`validation/wind_engineering/03_cedval_building/reference/c_star_points.csv` from the\n",
    "password-protected EWTL dataset. It is intentionally NOT executed in CI.**\n",
    "\n",
    "## What is validated\n",
    "\n",
    "- Scalar transport (advection-diffusion of the pollutant DDF).\n",
    "- The **voxel surface source** (on-body / near-ground emission), isolated\n",
    "  cleanly in this simple geometry.\n",
    "- Impermeable building faces (voxel zero-flux scalar BC `ScalarHWBB`).\n",
    "- LES + turbulent Schmidt number coupling (`Sc_t = 0.7`).\n",
    "\n",
    "## Scaling (model -> lattice)\n",
    "\n",
    "CEDVAL A1-5 building L x W x H = 100 x 150 x 125 mm (model scale 1:200).\n",
    "The height is resolved with H = 50 lattice cells, giving a cell size\n",
    "`dx = 125 mm / 50 = 2.5 mm/cell` and L = 40, W = 60 lattice cells. The\n",
    "domain is ~15 H x 10 H x 6 H, rounded to 752 x 504 x 304 cells.\n",
    "\n",
    "The dimensionless concentration follows the CEDVAL convention:\n",
    "\n",
    "$$c^{*} = \\frac{C \\, U_{ref} \\, H^{2}}{Q}$$\n",
    "\n",
    "## FLAGGED unverified parameters\n",
    "\n",
    "The following must be confirmed against the EWTL dataset before the\n",
    "absolute comparison is meaningful:\n",
    "\n",
    "- **`Q`** (emission rate): the case uses a unit Dirichlet wall value\n",
    "  `phi_w = 1.0`; the effective `Q` is recovered from the integrated wall\n",
    "  flux of the converged field (computed below). The physical `Q` from the\n",
    "  dataset is still needed to place `C` in absolute units.\n",
    "  `# VERIFY: CEDVAL A1-5 emission rate Q from the EWTL dataset`\n",
    "- **`U_ref` + reference height**: the lattice `U_ref = 0.06` is chosen for\n",
    "  Ma/Re only; the dataset's physical `U_ref` and reference height set the\n",
    "  `c*` normalisation. `# VERIFY: CEDVAL A1-5 U_ref value + reference height`\n",
    "\n",
    "## Reference\n",
    "\n",
    "CEDVAL A1-5, Univ. Hamburg EWTL (Leitl 1998/2000). Acceptance metrics per\n",
    "COST Action 732 / VDI 3783 Part 9 (Schatzmann et al. 2010; Chang & Hanna\n",
    "2004): FAC2 >= 0.5, NMSE < 4, FB in [-0.3, 0.3], hit-rate q >= 0.66\n",
    "(D = 0.25).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acae54e37e7d407bbb7b55eff062a284",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pathlib\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "try:\n",
    "    import pyvista as pv\n",
    "\n",
    "    pv.OFF_SCREEN = True\n",
    "except Exception:  # pragma: no cover - pyvista optional for metric-only runs\n",
    "    pv = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a63283cbaf04dbcab1f6479b197f3a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _find_project_root() -> pathlib.Path:\n",
    "    here = pathlib.Path.cwd().resolve()\n",
    "    for cand in [here, *here.parents]:\n",
    "        if (cand / \"pyproject.toml\").exists() and (cand / \"nassu\").is_dir():\n",
    "            return cand\n",
    "    raise RuntimeError(f\"Could not locate Nassu project root upward from {here}\")\n",
    "\n",
    "\n",
    "project_root = _find_project_root()\n",
    "case_path = (\n",
    "    project_root / \"validation/wind_engineering/03_cedval_building/03_cedval_building.nassu.yaml\"\n",
    ")\n",
    "results_dir = project_root / \"validation/wind_engineering/03_cedval_building/results\"\n",
    "comparison_dir = project_root / \"validation/wind_engineering/03_cedval_building/reference\"\n",
    "ref_csv = comparison_dir / \"c_star_points.csv\"\n",
    "print(\"case:\", case_path)\n",
    "print(\"results:\", results_dir)\n",
    "print(\"reference:\", ref_csv)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dd0d8092fe74a7c96281538738b07e2",
   "metadata": {},
   "source": [
    "## Scaling constants\n",
    "\n",
    "H, dx and the lattice reference velocity follow the YAML header. The\n",
    "physical `U_ref` and `Q` are FLAGGED placeholders to be confirmed from the\n",
    "dataset; for the dimensionless comparison we use the lattice values and the\n",
    "emission rate recovered from the wall flux.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72eea5119410473aa328ad9291626812",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lattice geometry (see the YAML header).\n",
    "H_LAT = 50  # building height in lattice cells\n",
    "L_LAT = 40  # along-wind length\n",
    "W_LAT = 60  # cross-wind width\n",
    "DX_MM = 2.5  # lattice cell size, mm/cell (model frame)\n",
    "BODY_X = 250  # building centre x (lattice)\n",
    "DOMAIN_W = 504  # domain width (lattice)\n",
    "PLANE_H = 5.01  # ground plane height (lattice)\n",
    "\n",
    "# Lattice reference velocity (top-of-profile). # VERIFY physical U_ref.\n",
    "U_REF_LAT = 0.06\n",
    "\n",
    "# Emission rate Q. The Dirichlet patch uses phi_w = 1.0; the effective\n",
    "# lattice Q is the integrated scalar wall flux of the converged field,\n",
    "# computed below. # VERIFY: absolute physical Q from the EWTL dataset.\n",
    "Q_LAT = None  # set from the integrated wall flux (see the flux cell)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8edb47106e1a46a883d545849b8ab81b",
   "metadata": {},
   "source": [
    "## Load the time-averaged result\n",
    "\n",
    "Loads the exported statistics / instantaneous XDMF for the converged\n",
    "`pollutant_phi` field. Adjust the glob to the actual exported artefact name\n",
    "once the GPU run exists.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10185d26023b46108eb7d9f57d49d2b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_timeavg_field():\n",
    "    \"\"\"Load the time-averaged pollutant field from the exported results.\n",
    "\n",
    "    Returns a pyvista MultiBlock (or dataset). The exact artefact name\n",
    "    depends on the export config; prefer the statistics mean export, else\n",
    "    fall back to the last instantaneous snapshot.\n",
    "    \"\"\"\n",
    "    if pv is None:\n",
    "        raise RuntimeError(\"pyvista is required to load the field\")\n",
    "    candidates = sorted(results_dir.rglob(\"*.xdmf\")) + sorted(results_dir.rglob(\"*.pvd\"))\n",
    "    if not candidates:\n",
    "        raise FileNotFoundError(f\"No XDMF/PVD result under {results_dir}. Run the GPU case first.\")\n",
    "    # Prefer a statistics/mean artefact when present.\n",
    "    stats = [c for c in candidates if \"stat\" in c.name.lower() or \"mean\" in c.name.lower()]\n",
    "    target = stats[-1] if stats else candidates[-1]\n",
    "    print(\"loading:\", target)\n",
    "    return pv.read(str(target))\n",
    "\n",
    "\n",
    "# field = load_timeavg_field()  # uncomment once results exist"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8763a12b2bbd4a93a75aff182afb95dc",
   "metadata": {},
   "source": [
    "## Recover the effective emission rate Q from the wall flux\n",
    "\n",
    "Scalar mass-balance check (acceptance criterion): the scalar emitted at the\n",
    "source equals what is advected out plus what accumulates. For a converged\n",
    "(statistically steady) field the accumulation term vanishes, so the emitted\n",
    "rate equals the net advective+diffusive outflux. We use the source wall flux\n",
    "as the effective lattice `Q` for the `c*` normalisation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7623eae2785240b9bd12b16a66d81610",
   "metadata": {},
   "outputs": [],
   "source": [
    "def effective_Q_from_field(field):\n",
    "    \"\"\"Estimate the effective lattice emission rate Q.\n",
    "\n",
    "    Integrates the scalar advective flux through a control surface enclosing\n",
    "    the source (or, equivalently, the net outflux through the outlet plane).\n",
    "    Implementation depends on the available exported macrs; placeholder\n",
    "    until the GPU result schema is fixed.\n",
    "    \"\"\"\n",
    "    raise NotImplementedError(\n",
    "        \"Populate from the converged field: integrate the scalar flux \"\n",
    "        \"through a control surface around the source / outlet plane.\"\n",
    "    )\n",
    "\n",
    "\n",
    "# Q_LAT = effective_Q_from_field(field)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cdc8c89c7104fffa095e18ddfef8986",
   "metadata": {},
   "source": [
    "## Sample the predicted c* at the dataset sampling points\n",
    "\n",
    "Maps the CEDVAL model-frame sampling coordinates (mm, building base origin)\n",
    "to lattice positions via `dx` and the building placement, samples the\n",
    "simulated `pollutant_phi`, and normalises to `c* = C * U_ref * H^2 / Q`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b118ea5561624da68c537baed56e602f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_mm_to_lattice(x_mm, y_mm, z_mm):\n",
    "    \"\"\"Map CEDVAL model-frame coords (mm, building base origin) to lattice.\n",
    "\n",
    "    The building base is centred at (BODY_X, DOMAIN_W/2) on the ground\n",
    "    plane (z = PLANE_H). Model +x is along-wind, +y cross-wind, +z up.\n",
    "    \"\"\"\n",
    "    x = BODY_X + x_mm / DX_MM\n",
    "    y = DOMAIN_W / 2.0 + y_mm / DX_MM\n",
    "    z = PLANE_H + z_mm / DX_MM\n",
    "    return np.column_stack([x, y, z])\n",
    "\n",
    "\n",
    "def predicted_c_star(field, ref_df, q_lat, u_ref=U_REF_LAT, h_lat=H_LAT):\n",
    "    pts = model_mm_to_lattice(\n",
    "        ref_df[\"x_model_mm\"].to_numpy(),\n",
    "        ref_df[\"y_model_mm\"].to_numpy(),\n",
    "        ref_df[\"z_model_mm\"].to_numpy(),\n",
    "    )\n",
    "    poly = pv.PolyData(pts)\n",
    "    sampled = poly.sample(field)\n",
    "    C = np.asarray(sampled[\"pollutant_phi\"])\n",
    "    return C * u_ref * h_lat**2 / q_lat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "938c804e27f84196a10c8828c723f798",
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_df = pd.read_csv(ref_csv, comment=\"#\")\n",
    "print(f\"reference sampling points: {len(ref_df)}\")\n",
    "if len(ref_df) == 0:\n",
    "    print(\n",
    "        \"Reference CSV is empty (placeholder). Populate \"\n",
    "        \"validation/wind_engineering/03_cedval_building/reference/c_star_points.csv from the \"\n",
    "        \"EWTL dataset before running the comparison.\"\n",
    "    )\n",
    "ref_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "504fb2a444614c0babb325280ed9130a",
   "metadata": {},
   "source": [
    "## COST 732 / VDI 3783-9 validation metrics\n",
    "\n",
    "`FAC2`, `NMSE`, `FB`, `R` and the hit-rate `q` (D = 0.25) over the paired\n",
    "(observed, predicted) `c*` samples. Definitions follow Chang & Hanna (2004)\n",
    "and VDI 3783 Part 9 (2005).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59bbdb311c014d738909a11f9e486628",
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost732_metrics(obs, pred, hitrate_D=0.25):\n",
    "    obs = np.asarray(obs, dtype=float)\n",
    "    pred = np.asarray(pred, dtype=float)\n",
    "    eps = np.finfo(float).tiny\n",
    "    mo, mp = obs.mean(), pred.mean()\n",
    "    fb = (mo - mp) / (0.5 * (mo + mp))\n",
    "    nmse = np.mean((obs - pred) ** 2) / (mo * mp)\n",
    "    ratio = pred / np.where(np.abs(obs) < eps, eps, obs)\n",
    "    fac2 = np.mean((ratio >= 0.5) & (ratio <= 2.0))\n",
    "    r = np.corrcoef(obs, pred)[0, 1] if len(obs) > 1 else np.nan\n",
    "    # VDI 3783-9 hit-rate: relative tolerance D (absolute fallback W omitted).\n",
    "    q = np.mean(np.abs(pred - obs) <= hitrate_D * np.abs(obs))\n",
    "    return {\"FAC2\": fac2, \"NMSE\": nmse, \"FB\": fb, \"R\": r, \"q\": q}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b43b363d81ae4b689946ece5c682cd59",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run once results + reference exist:\n",
    "#\n",
    "# field = load_timeavg_field()\n",
    "# Q_LAT = effective_Q_from_field(field)\n",
    "# pred = predicted_c_star(field, ref_df, Q_LAT)\n",
    "# obs = ref_df[\"c_star_obs\"].to_numpy()\n",
    "# metrics = cost732_metrics(obs, pred)\n",
    "# print(metrics)\n",
    "#\n",
    "# assert metrics[\"FAC2\"] >= 0.5, f\"FAC2 {metrics['FAC2']:.3f} < 0.5\"\n",
    "# assert metrics[\"NMSE\"] < 4.0, f\"NMSE {metrics['NMSE']:.3f} >= 4\"\n",
    "# assert -0.3 <= metrics[\"FB\"] <= 0.3, f\"FB {metrics['FB']:.3f} outside [-0.3, 0.3]\"\n",
    "# assert metrics[\"q\"] >= 0.66, f\"hit-rate q {metrics['q']:.3f} < 0.66\"\n",
    "print(\"Metrics cell ready (commented out until GPU results + dataset exist).\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a65eabff63a45729fe45fb5ade58bdc",
   "metadata": {},
   "source": [
    "## Plots\n",
    "\n",
    "(1) Surface / near-ground concentration field on the vertical symmetry\n",
    "plane and a pedestrian-level horizontal plane; (2) a scatter of predicted\n",
    "vs observed `c*` with the FAC2 band.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3933fab20d04ec698c2621248eb3be0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_scatter(obs, pred):\n",
    "    fig, ax = plt.subplots(figsize=(5, 5))\n",
    "    lo = min(obs.min(), pred.min())\n",
    "    hi = max(obs.max(), pred.max())\n",
    "    ax.scatter(obs, pred, s=20, alpha=0.7, label=\"sampling points\")\n",
    "    ax.plot([lo, hi], [lo, hi], \"k-\", lw=1, label=\"1:1\")\n",
    "    ax.plot([lo, hi], [0.5 * lo, 0.5 * hi], \"k--\", lw=0.8, label=\"FAC2 band\")\n",
    "    ax.plot([lo, hi], [2.0 * lo, 2.0 * hi], \"k--\", lw=0.8)\n",
    "    ax.set_xlabel(\"observed c*  (CEDVAL A1-5)\")\n",
    "    ax.set_ylabel(\"predicted c*  (Nassu)\")\n",
    "    ax.set_title(\"Predicted vs observed dimensionless concentration\")\n",
    "    ax.legend()\n",
    "    ax.set_aspect(\"equal\", \"box\")\n",
    "    return fig\n",
    "\n",
    "\n",
    "# fig = plot_scatter(obs, pred)  # uncomment once results exist\n",
    "print(\"Plot helpers ready.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dd4641cc4064e0191573fe9c69df29b",
   "metadata": {},
   "source": [
    "## Version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8309879909854d7188b41380fd92a7c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from nassu.cfg.model import ConfigScheme\n",
    "\n",
    "sim_cfg = next(iter(ConfigScheme.sim_cfgs_from_file_dct(str(case_path)).values()))\n",
    "sim_info = sim_cfg.output.read_info()\n",
    "\n",
    "nassu_commit = sim_info[\"commit\"]\n",
    "nassu_version = sim_info[\"version\"]\n",
    "print(\"Version:\", nassu_version)\n",
    "print(\"Commit hash:\", nassu_commit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ed186c9a28b402fb0bc4494df01f08d",
   "metadata": {},
   "source": [
    "## Configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb1e1581032b452c9409d6c6813c49d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Code\n",
    "\n",
    "Code(filename=str(case_path))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "nassu",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}