{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Case 06: Stratified isolated building (Zhou et al. 2021)\n",
    "\n",
    "> **STATUS: SCAFFOLD - unrun.** This notebook is analysis structure\n",
    "> only. It runs before the GPU result and the digitized reference\n",
    "> data exist (every loader is guarded), and it makes NO physics claim.\n",
    "> Reference-data digitization + GPU validation are pending (epic\n",
    "> #1034 P5, issue #1039).\n",
    "\n",
    "Validation of the stratified isolated-building case for\n",
    "`validation/wind_engineering/06_stratified_building/06_stratified_building.nassu.yaml`\n",
    "(`stratifiedBuildingZhouC2RiN01`, Case 2, Ri = -0.1). An isolated\n",
    "rectangular building in an unstable ABL, with a ground hole source\n",
    "behind the leeward wall and a buoyant temperature scalar driving the\n",
    "Monin-Obukhov ground wall model.\n",
    "\n",
    "**Reference:** Zhou, X. et al. (2021), *J. Wind Eng. Ind. Aerodyn.*\n",
    "211, 104526 (TPU Database 2006 wind tunnel for Cases 1 and 2)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scaling and normalization\n",
    "\n",
    "Building height H = 0.16 m, footprint 0.5H x 0.5H, u_H = 1.37 m/s,\n",
    "Re = 14811, Ri = -0.1 (Case 2). In lattice units (level 0) H_lat = 32,\n",
    "u_H_lat = 0.05 (Ma = 0.087). Reported quantities:\n",
    "\n",
    "- velocity <u> / u_H,\n",
    "- temperature phi = (T - T_r) / dT, dT = 34 K,\n",
    "- concentration <c> / c0 with c0 = 0.05 q / (u_H H^2).\n",
    "\n",
    "Only Case 1 (Ri=0) and Case 2 (Ri=-0.1) are wind-tunnel-backed; Cases\n",
    "3-6 are CFD-only (Zhou's own LES). Label every reference series with its\n",
    "case and never present a CFD-only case as measurement data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pathlib\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "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 from {here}\")\n",
    "\n",
    "\n",
    "project_root = _find_project_root()\n",
    "case = 'validation/wind_engineering/06_stratified_building'\n",
    "results_dir = project_root / 'output' / case\n",
    "reference_dir = project_root / case / 'reference'\n",
    "\n",
    "# Scaling constants (see the case YAML header).\n",
    "H_LAT = 32.0\n",
    "U_H_LAT = 0.05\n",
    "PLANE_HEIGHT = 5.01\n",
    "DT_K = 34.0\n",
    "print('results dir:', results_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load model outputs (guarded)\n",
    "\n",
    "The case records the inlet line, the behind-building lines a,b,c, the\n",
    "probe points A-E, the symmetry plane and the full-domain statistics.\n",
    "Wire these readers to the actual export format once the GPU run exists;\n",
    "until then they return `None` and the plots degrade gracefully."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_series(export_name: str, entity: str):\n",
    "    \"\"\"Return a model series or None if absent.\"\"\"\n",
    "    if not results_dir.exists():\n",
    "        return None\n",
    "    # TODO(P5): point this at the real line/point/plane export and\n",
    "    # time-average over the statistics window.\n",
    "    hits = list(results_dir.rglob(f'*{export_name}*{entity}*'))\n",
    "    if not hits:\n",
    "        return None\n",
    "    raise NotImplementedError('Wire the series reader to the run output.')\n",
    "\n",
    "\n",
    "model_lines = {ln: load_series('behind_building_lines', f'line_{ln}')\n",
    "               for ln in ('a', 'b', 'c')}\n",
    "model_inlet = load_series('inlet_profile', 'inlet')\n",
    "have_model = any(v is not None for v in model_lines.values())\n",
    "print('model data present:', have_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load digitized reference data (guarded)\n",
    "\n",
    "Reference CSVs live under `reference/`. They are EMPTY PLACEHOLDERS\n",
    "until digitized from the Zhou et al. figures / TPU Database 2006 at P5\n",
    "(see `reference/README.md`). Each carries a `source` column naming its\n",
    "figure and case; that source is used verbatim in every plot legend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_reference(name: str):\n",
    "    path = reference_dir / name\n",
    "    if not path.exists():\n",
    "        return None\n",
    "    df = pd.read_csv(path, comment='#')\n",
    "    return df if not df.empty else None\n",
    "\n",
    "\n",
    "def source_label(df, fallback='Zhou et al. (2021)'):\n",
    "    if df is not None and 'source' in df and len(df):\n",
    "        return str(df['source'].iloc[0])\n",
    "    return fallback\n",
    "\n",
    "\n",
    "ref_inlet = load_reference('inlet_profiles.csv')\n",
    "ref_lines = load_reference('behind_building_lines.csv')\n",
    "ref_recirc = load_reference('recirc_length_vs_Ri.csv')\n",
    "have_ref = ref_lines is not None\n",
    "print('reference data present:', have_ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Behind-building lines a, b, c (Fig. 7)\n",
    "\n",
    "Streamwise velocity and concentration at x/H = 0.375, 0.625, 1.0 behind\n",
    "the leeward wall vs height, model against the digitized Zhou et al.\n",
    "(2021) Case 2 data with its reported +/-15% band. Plots only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)\n",
    "for ax, ln in zip(axes, ('a', 'b', 'c')):\n",
    "    if ref_lines is not None and 'line' in ref_lines:\n",
    "        sub = ref_lines[ref_lines['line'] == ln]\n",
    "        if len(sub):\n",
    "            err = sub['err_frac'] * sub['c_over_c0'] if 'err_frac' in sub else None\n",
    "            ax.errorbar(sub['c_over_c0'], sub['z_over_H'], xerr=err, fmt='ks',\n",
    "                        label=source_label(ref_lines))\n",
    "    if model_lines.get(ln) is not None:\n",
    "        pass  # TODO(P5): overlay model <c>/c0 vs z/H for this line.\n",
    "    ax.set_title(f'line {ln}')\n",
    "    ax.set_xlabel('<c> / c0')\n",
    "    if ref_lines is None and not have_model:\n",
    "        ax.text(0.5, 0.5, 'pending', ha='center', va='center',\n",
    "                transform=ax.transAxes)\n",
    "    if ax.get_legend_handles_labels()[0]:\n",
    "        ax.legend()\n",
    "axes[0].set_ylabel('z / H')\n",
    "fig.suptitle('Behind-building concentration (Zhou et al. 2021, Fig. 7) - SCAFFOLD')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inlet profiles (Fig. 5)\n",
    "\n",
    "Inlet u, k and temperature profiles vs the TPU Database 2006. Tests the\n",
    "provisional SEM + mean-T inflow (P5 fits the inlet to the TPU profiles)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)\n",
    "for ax, col, title in zip(axes, ['U_over_uH', 'k_over_uH2', 'phi_T'],\n",
    "                          ['u / u_H', 'k / u_H^2', 'phi_T']):\n",
    "    if ref_inlet is not None and col in ref_inlet:\n",
    "        ax.plot(ref_inlet[col], ref_inlet['z_over_H'], 'k^',\n",
    "                label=source_label(ref_inlet, 'TPU Database 2006 (Zhou et al. 2021)'))\n",
    "    if model_inlet is not None:\n",
    "        pass  # TODO(P5): overlay model inlet profile.\n",
    "    ax.set_title(title)\n",
    "    ax.set_xlabel(title)\n",
    "    if ref_inlet is None and model_inlet is None:\n",
    "        ax.text(0.5, 0.5, 'pending', ha='center', va='center',\n",
    "                transform=ax.transAxes)\n",
    "    if ax.get_legend_handles_labels()[0]:\n",
    "        ax.legend()\n",
    "axes[0].set_ylabel('z / H')\n",
    "fig.suptitle('Inlet profiles vs TPU Database 2006 (Zhou et al. 2021, Fig. 5) - SCAFFOLD')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recirculation length vs Ri (Fig. 9)\n",
    "\n",
    "Normalized recirculation length behind the building vs Richardson\n",
    "number (Zhou reports it shrinking ~32% from Ri 0 to -1.5). Model points\n",
    "come from the Ri ladder (P5); the reference is the digitized Zhou\n",
    "series, with wind-tunnel-backed points (Ri 0, -0.1) marked distinctly\n",
    "from the CFD-only points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(7, 4))\n",
    "if ref_recirc is not None and 'Ri' in ref_recirc:\n",
    "    wt = ref_recirc[ref_recirc.get('backing') == 'wind_tunnel']\n",
    "    cfd = ref_recirc[ref_recirc.get('backing') == 'CFD_only']\n",
    "    if len(wt):\n",
    "        ax.plot(wt['Ri'], wt['Lr_over_H'], 'ks',\n",
    "                label=source_label(ref_recirc, 'Zhou et al. 2021 (wind tunnel)'))\n",
    "    if len(cfd):\n",
    "        ax.plot(cfd['Ri'], cfd['Lr_over_H'], 'kx',\n",
    "                label='Zhou et al. 2021 (LES, CFD-only)')\n",
    "# TODO(P5): overlay Nassu ladder points (Ri vs recirc length).\n",
    "ax.set_xlabel('Ri')\n",
    "ax.set_ylabel('L_r / H')\n",
    "ax.set_title('Recirculation length vs Ri (Zhou et al. 2021, Fig. 9) - SCAFFOLD')\n",
    "if ref_recirc is None:\n",
    "    ax.text(0.5, 0.5, 'pending GPU ladder + reference digitization',\n",
    "            ha='center', va='center', transform=ax.transAxes)\n",
    "if ax.get_legend_handles_labels()[0]:\n",
    "    ax.legend()\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fields + PSDs (Figs. 8, 10-13)\n",
    "\n",
    "Symmetry-plane streamlines/contours (Fig. 8), vertical velocity\n",
    "(Fig. 10), u'w'/u_H^2 (Fig. 11) and the point PSDs at A-E (Figs. 12-13).\n",
    "Wire the plane and point-series readers at P5 and render contourf /\n",
    "spectrum plots (not tables)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO(P5): load 'symmetry_plane/vertical_symmetry' -> contourf u/T/c;\n",
    "# load 'probe_points' A-E -> Welch PSD vs Strouhal, one line per point,\n",
    "# each legend-labelled with its measurement point and the Zhou source.\n",
    "print('Field + PSD plots: pending GPU result.')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}