{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "steckler-status",
   "metadata": {},
   "source": [
    "> **Status: planned - not yet validated against reference data.** The case is\n",
    "> posed on the variable-density low-Mach closure (`models.low_mach`, Taha et al.\n",
    "> 2024): EOS density `rho = P/(r T)`, the exact `(rho - rho_inf) g` buoyancy and\n",
    "> a finite-difference energy field carrying the burner heat-release source. A\n",
    "> validated run also needs closure stability at the fire-scale density ratio, a\n",
    "> thermal wall BC on the voxelized room shell, an EOS-consistent open doorway,\n",
    "> and the 62.9 kW heat-release-to-lattice mapping; see the case `README.md`.\n",
    ">\n",
    "> The analysis below reads the low-Mach `temp` energy macroscopic and the\n",
    "> `models.low_mach` configuration, and reconstructs the dimensional profiles\n",
    "> through the variable-density route (the EOS, the emergent real temperature\n",
    "> ratio). The plotted numbers are finalised with the case's first validated GPU\n",
    "> run, when the burner `rate` and lattice gravity are calibrated to the matched\n",
    "> `Q*_H` and the steady-state averaging window is measured."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "700165a6",
   "metadata": {},
   "source": [
    "# Case 03 - Steckler room fire (Steckler, Quintiere & Rinkinen 1982)\n",
    "\n",
    "Validation of the open, internally-heated **compartment** flow on the\n",
    "variable-density **low-Mach** closure (`models.low_mach`, Taha et al. 2024): a\n",
    "single room with a doorway and a steady floor-centred heat release develops a\n",
    "stratified two-layer flow that spills hot fluid out of the top of the doorway\n",
    "and draws ambient air in through the bottom, with a **neutral plane** in between.\n",
    "\n",
    "**Scope.** This is posed as a **real fire**, not a reduced-ratio Boussinesq\n",
    "surrogate: the hot layer carries its real `T_hot / T_inf ~ 2`, the temperature\n",
    "closes the equation of state `rho = P / (r T)`, and the resulting density\n",
    "variation drives the flow through the exact `(rho - rho_inf) g` buoyancy (see\n",
    "the case `README.md`). The heat release is a volumetric source on the\n",
    "finite-difference energy field and the comparison is on **non-dimensional**\n",
    "quantities: the doorway velocity / temperature profile shapes and the\n",
    "neutral-plane height fraction `z_n / H_door`.\n",
    "\n",
    "The notebook time-averages the dense doorway `plane_series`, extracts the\n",
    "door-centreline vertical profiles and the neutral-plane height, and compares\n",
    "them against the digitised Steckler reference when it is available\n",
    "(`reference/doorway_profiles.csv`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "206572e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pathlib\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import nassu.viz as common\n",
    "from nassu.cfg.model import ConfigScheme\n",
    "\n",
    "common.use_style()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2891d2e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import importlib.util\n",
    "import os\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 upward from {here}\")\n",
    "\n",
    "\n",
    "project_root = _find_project_root()\n",
    "# The simulation `save_path` (and therefore every probe/data path read from the\n",
    "# config below) is repo-root-relative, so run from the project root regardless\n",
    "# of where the notebook was launched (case dir under Jupyter / nbconvert).\n",
    "os.chdir(project_root)\n",
    "case_dir = project_root / \"validation/buoyant_flows/03_steckler_room_fire\"\n",
    "case_path = case_dir / \"03_steckler_room_fire.nassu.yaml\"\n",
    "reference_csv = case_dir / \"reference/doorway_profiles.csv\"\n",
    "\n",
    "# Geometry is in real metres and the grid is driven by one knob, `dx` (level-0\n",
    "# lattice spacing). Read `dx` from the YAML and map physical (metre) coordinates\n",
    "# to level-0 lattice with the case's `lattice_mapping` helper, so the analysis\n",
    "# follows `dx` automatically (no hardcoded lattice numbers).\n",
    "_spec = importlib.util.spec_from_file_location(\n",
    "    \"lattice_mapping\", case_dir / \"geometry\" / \"lattice_mapping.py\"\n",
    ")\n",
    "lm = importlib.util.module_from_spec(_spec)\n",
    "_spec.loader.exec_module(lm)\n",
    "\n",
    "from nassu.cfg.loader import load_cfg_yaml\n",
    "\n",
    "dx = float(load_cfg_yaml(case_path)[\"variables\"][\"dx\"])\n",
    "\n",
    "# Doorway opening in lattice: x in [-0.37, 0.37] m, z in [0, 1.83] m, in the\n",
    "# front wall y = 0 (door centre at the room centre x = 0).\n",
    "DOOR_X = (lm.to_lat(-0.37, \"x\", dx), lm.to_lat(0.37, \"x\", dx))\n",
    "H_DOOR = 1.83 / dx  # doorway height (z, lattice units)\n",
    "DOOR_XC = lm.to_lat(0.0, \"x\", dx)  # door centreline x (= room centre)\n",
    "\n",
    "(name, sim_id), cfg = next(iter(ConfigScheme.sim_cfgs_from_file_dct(str(case_path)).items()))\n",
    "print(f\"loaded sim {name} (#{sim_id}); dx = {dx} m, door height H_DOOR = {H_DOOR:g} lattice\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75145ff5",
   "metadata": {},
   "source": [
    "## Load the steady doorway plane\n",
    "\n",
    "The doorway `plane_series` (plane normal $y$, sitting in the front-wall doorway)\n",
    "carries the through-door velocity component `uy` and the low-Mach energy-field\n",
    "temperature `temp` over the $x$-$z$ opening. We average all snapshots after the\n",
    "spin-up to get the statistically-steady field, exactly as the other thermal\n",
    "plane notebooks do."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e227a034",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_doorway_plane(cfg):\n",
    "    \"\"\"Return (xs, zs, T, uy) time-mean grids for the doorway plane.\n",
    "\n",
    "    The plane normal is y, so the in-plane coordinates are x and z and the\n",
    "    through-door velocity is the out-of-plane component `uy`. The doorway faces\n",
    "    -y (the plenum is at y < 0, the room at y > 0), so `uy > 0` is inflow (cold\n",
    "    ambient drawn into the room through the lower door) and `uy < 0` is outflow\n",
    "    (hot layer spilling out of the upper door). Averages all snapshots after 60%\n",
    "    of the run.\n",
    "    \"\"\"\n",
    "    plane = cfg.output.exports[\"plane_series\"].series.planes[\"doorway\"]\n",
    "    pts = pd.read_csv(plane.points_filename)\n",
    "\n",
    "    xs = np.sort(pts[\"x\"].unique())\n",
    "    zs = np.sort(pts[\"z\"].unique())\n",
    "    ix = {v: i for i, v in enumerate(xs)}\n",
    "    iz = {v: i for i, v in enumerate(zs)}\n",
    "    cols = [str(int(i)) for i in pts[\"idx\"]]\n",
    "    col_xi = np.array([ix[x] for x in pts[\"x\"]])\n",
    "    col_zi = np.array([iz[z] for z in pts[\"z\"]])\n",
    "\n",
    "    def _mean_grid(macr):\n",
    "        df = plane.read_full_data(macr)\n",
    "        start = int(df[\"time_step\"].max() * 0.6)  # drop the spin-up\n",
    "        df = df[df[\"time_step\"] >= start]\n",
    "        arr = df[cols].to_numpy(dtype=float).mean(axis=0)\n",
    "        grid = np.full((zs.size, xs.size), np.nan)\n",
    "        grid[col_zi, col_xi] = arr\n",
    "        return grid\n",
    "\n",
    "    return xs, zs, _mean_grid(\"temp\"), _mean_grid(\"uy\")\n",
    "\n",
    "\n",
    "xs, zs, T_grid, uy_grid = load_doorway_plane(cfg)\n",
    "print(f\"doorway plane grid: {xs.size} (x) x {zs.size} (z)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e1d4e1f",
   "metadata": {},
   "source": [
    "## Door-centreline profiles and the neutral plane\n",
    "\n",
    "At the door centreline (the room centre $x = 0$ m) we take the vertical profiles\n",
    "of the through-door velocity and the temperature over the opening height,\n",
    "normalise them, and locate the **neutral plane** as the height where the\n",
    "through-door velocity changes sign (outflow above, inflow below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5f33d37",
   "metadata": {},
   "outputs": [],
   "source": [
    "def door_centreline_profiles(xs, zs, T_grid, uy_grid):\n",
    "    \"\"\"Vertical profiles at the door centreline over the opening height.\"\"\"\n",
    "    xi = int(np.argmin(np.abs(xs - DOOR_XC)))\n",
    "    mask = (zs >= 0.0) & (zs <= H_DOOR)\n",
    "    z = zs[mask]\n",
    "    uy = uy_grid[mask, xi]\n",
    "    T = T_grid[mask, xi]\n",
    "    return z, uy, T\n",
    "\n",
    "\n",
    "def neutral_plane_height(z, uy):\n",
    "    \"\"\"Height where the through-door velocity crosses zero (linear interp).\"\"\"\n",
    "    s = np.sign(uy)\n",
    "    crossings = np.where(np.diff(s) != 0)[0]\n",
    "    if crossings.size == 0:\n",
    "        return np.nan\n",
    "    # Take the dominant crossing (largest velocity swing around it).\n",
    "    k = crossings[np.argmax(np.abs(uy[crossings + 1] - uy[crossings]))]\n",
    "    z0, z1 = z[k], z[k + 1]\n",
    "    u0, u1 = uy[k], uy[k + 1]\n",
    "    return z0 - u0 * (z1 - z0) / (u1 - u0)\n",
    "\n",
    "\n",
    "z, uy, T = door_centreline_profiles(xs, zs, T_grid, uy_grid)\n",
    "z_n = neutral_plane_height(z, uy)\n",
    "\n",
    "# Non-dimensionalise: height by door height; velocity by its peak magnitude;\n",
    "# temperature by the centreline span (ambient -> hot layer). This matches the\n",
    "# normalisation of the digitised reference (`reference/doorway_profiles.csv`).\n",
    "z_over_H = z / H_DOOR\n",
    "u_norm = uy / np.nanmax(np.abs(uy))\n",
    "T_span = np.nanmax(T) - np.nanmin(T)\n",
    "T_norm = (T - np.nanmin(T)) / T_span if T_span > 0 else T * 0.0\n",
    "\n",
    "print(f\"neutral-plane height: z_n = {z_n:.1f} lattice  ->  z_n / H_door = {z_n / H_DOOR:.3f}\")\n",
    "print(\"(Steckler test-16 experiment: z_n / H_door = 0.560)\")\n",
    "\n",
    "# Through-door mass balance (uniform spacing). uy > 0 is inflow, uy < 0 is\n",
    "# outflow; np.trapezoid (NumPy 2.x; np.trapz was removed).\n",
    "inflow = np.trapezoid(np.clip(uy, 0.0, None), z)\n",
    "outflow = -np.trapezoid(np.clip(uy, None, 0.0), z)\n",
    "print(\n",
    "    f\"door flux balance: inflow={inflow:.3g}  outflow={outflow:.3g}  \"\n",
    "    f\"imbalance={abs(outflow - inflow) / max(outflow, 1e-30):.1%}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0a32dd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reference overlay, read defensively (empty until the profiles are digitised).\n",
    "ref = None\n",
    "try:\n",
    "    _r = pd.read_csv(reference_csv, comment=\"#\")\n",
    "    if len(_r) > 0:\n",
    "        ref = _r\n",
    "except Exception as exc:  # noqa: BLE001\n",
    "    print(f\"reference not loaded ({exc}); plotting Nassu profiles alone\")\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10, 5), sharey=True)\n",
    "\n",
    "axes[0].plot(u_norm, z_over_H, \"-o\", ms=3, label=\"Nassu\")\n",
    "if ref is not None:\n",
    "    axes[0].plot(ref[\"u_norm\"], ref[\"z_over_H\"], \"ks\", ms=4, label=\"Steckler 1982\")\n",
    "axes[0].axvline(0.0, color=\"0.6\", lw=0.8)\n",
    "if np.isfinite(z_n):\n",
    "    axes[0].axhline(z_n / H_DOOR, color=\"C3\", ls=\"--\", lw=1, label=\"neutral plane\")\n",
    "axes[0].set_xlabel(\"through-door velocity (normalised)\")\n",
    "axes[0].set_ylabel(r\"$z / H_\\mathrm{door}$\")\n",
    "axes[0].legend()\n",
    "\n",
    "axes[1].plot(T_norm, z_over_H, \"-o\", ms=3, label=\"Nassu\")\n",
    "if ref is not None:\n",
    "    axes[1].plot(ref[\"T_norm\"], ref[\"z_over_H\"], \"ks\", ms=4, label=\"Steckler 1982\")\n",
    "if np.isfinite(z_n):\n",
    "    axes[1].axhline(z_n / H_DOOR, color=\"C3\", ls=\"--\", lw=1)\n",
    "axes[1].set_xlabel(\"temperature (normalised)\")\n",
    "axes[1].legend()\n",
    "\n",
    "fig.suptitle(\"Steckler doorway profiles (door centreline)\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7398d34",
   "metadata": {},
   "source": [
    "## Dimensional reconstruction by similarity\n",
    "\n",
    "The run is set up for **dynamic similarity** with Steckler test 16: the\n",
    "dimensionless heat-release rate $Q^*_H$, the Prandtl number and the geometry\n",
    "length-ratios all match the experiment (see the case `README.md`), so the\n",
    "*dimensionless* fields the solver produces are the experiment's. That means the\n",
    "**dimensional** doorway profiles can be reconstructed and overlaid on the raw\n",
    "NIST/FDS measurements in $^\\circ$C and m/s.\n",
    "\n",
    "Because the closure is variable-density, the lattice temperature is already the\n",
    "temperature **ratio** $\\theta = T/T_\\mathrm{inf}$ (the ambient is the lattice\n",
    "reference $T_\\mathrm{ref} = 1$), so the dimensional temperature is simply\n",
    "$T = \\theta\\,T_\\mathrm{inf}$ - there is no Boussinesq\n",
    "$\\theta_\\mathrm{eff} = \\beta\\,\\phi$ expansion-coefficient step. The hot-layer\n",
    "excess ratio\n",
    "$\\theta_\\mathrm{hot} - 1 = (T_\\mathrm{hot} - T_\\mathrm{inf})/T_\\mathrm{inf}$ is\n",
    "read **emergently** off the door-centreline temperature (not fitted), and the\n",
    "density follows the EOS $\\rho = P/(rT) = 1/\\theta$ - the warm layer thins to\n",
    "$\\rho_\\mathrm{hot} \\approx \\rho_\\mathrm{inf}/2$ at $\\theta_\\mathrm{hot} \\approx 2$.\n",
    "The reconstruction scales are\n",
    "\n",
    "$$\\Delta T = (\\theta_\\mathrm{hot}-1)\\,T_\\mathrm{inf}, \\qquad\n",
    "  U_b = \\sqrt{g\\,H_\\mathrm{door}\\,(\\theta_\\mathrm{hot}-1)},$$\n",
    "\n",
    "with $T_\\mathrm{inf} = 301$ K (28 $^\\circ$C ambient, the experiment) and the\n",
    "physical $g = 9.81$ m/s$^2$, $H_\\mathrm{door} = 1.83$ m. With a real temperature\n",
    "ratio this reconstruction carries no Boussinesq linearization error; the\n",
    "residual modelling limit is the reduced Reynolds number. The temperature scale\n",
    "grows like $(\\theta_\\mathrm{hot}-1)$ but the velocity scale only like its square\n",
    "root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff447bc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Variable-density similarity reconstruction: turn the matched-Q* dimensionless\n",
    "# profiles into dimensional ones (deg C, m/s) using the run's OWN emergent\n",
    "# hot-layer temperature RATIO (no Boussinesq beta), and overlay on the raw\n",
    "# NIST/FDS dimensional data.\n",
    "T0_C = 28.0  # ambient temperature (deg C), Steckler test 16\n",
    "T_INF_K = 301.15  # ambient absolute temperature (K); the lattice T_ref = 1 maps here\n",
    "G_PHYS = 9.81  # gravity (m/s^2)\n",
    "H_DOOR_M = 1.83  # door height (m)\n",
    "\n",
    "# In the variable-density closure the lattice temperature IS the ratio\n",
    "# theta = T / T_inf (ambient T_ref = 1), so the EOS density is rho = 1 / theta\n",
    "# and the emergent hot-layer excess ratio is read directly off the centreline.\n",
    "theta_inf = 1.0\n",
    "theta_hot = np.nanmax(T)\n",
    "excess = theta_hot - theta_inf  # (T_hot - T_inf) / T_inf, dimensionless\n",
    "\n",
    "# Predicted physical scales (not fitted to the experiment).\n",
    "dT_pred = excess * T_INF_K  # hot-layer temperature rise (K)\n",
    "Ub_pred = np.sqrt(G_PHYS * H_DOOR_M * excess)  # door buoyant velocity (m/s)\n",
    "\n",
    "print(f\"emergent theta_hot = T_hot/T_inf = {theta_hot:.3f}   (experiment ~ 2)\")\n",
    "print(f\"EOS hot-layer density rho_hot = 1/theta_hot = {1.0 / theta_hot:.3f}   (~ 0.5)\")\n",
    "print(f\"-> predicted dT = {dT_pred:.0f} K     (experiment ~ 90 K)\")\n",
    "print(f\"-> predicted U_b = {Ub_pred:.2f} m/s  (experiment ~ 2.3 m/s, measured peak 2.12)\")\n",
    "\n",
    "# Reconstructed Nassu dimensional profiles. Temperature is EOS-direct\n",
    "# (T = theta * T_inf); the through-door velocity uses the buoyant scale (nassu\n",
    "# convention: u > 0 = inflow).\n",
    "u_phys = u_norm * Ub_pred  # m/s\n",
    "T_phys = T * T_INF_K - 273.15  # deg C, straight from the temperature ratio\n",
    "\n",
    "# Raw experimental dimensional profiles (door-centreline column *_C).\n",
    "fds = pd.read_csv(case_dir / \"reference\" / \"steckler_test16_fds_exp.csv\")\n",
    "fds = fds[(fds[\"Height\"] <= H_DOOR_M) & fds[\"V_C\"].notna() & fds[\"T_C\"].notna()].sort_values(\n",
    "    \"Height\"\n",
    ")\n",
    "z_over_H_exp = fds[\"Height\"].to_numpy() / H_DOOR_M\n",
    "u_exp = -fds[\"V_C\"].to_numpy()  # flip FDS (outflow +) to nassu (inflow +)\n",
    "T_exp = fds[\"T_C\"].to_numpy()  # deg C\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10, 5), sharey=True)\n",
    "\n",
    "axes[0].plot(u_phys, z_over_H, \"-o\", ms=3, label=\"Nassu (reconstructed)\")\n",
    "axes[0].plot(u_exp, z_over_H_exp, \"ks\", ms=4, label=\"Steckler 1982 (NIST/FDS)\")\n",
    "axes[0].axvline(0.0, color=\"0.6\", lw=0.8)\n",
    "if np.isfinite(z_n):\n",
    "    axes[0].axhline(z_n / H_DOOR, color=\"C3\", ls=\"--\", lw=1, label=\"neutral plane\")\n",
    "axes[0].set_xlabel(\"through-door velocity (m/s)\")\n",
    "axes[0].set_ylabel(r\"$z / H_\\mathrm{door}$\")\n",
    "axes[0].legend()\n",
    "\n",
    "axes[1].plot(T_phys, z_over_H, \"-o\", ms=3, label=\"Nassu (reconstructed)\")\n",
    "axes[1].plot(T_exp, z_over_H_exp, \"ks\", ms=4, label=\"Steckler 1982 (NIST/FDS)\")\n",
    "if np.isfinite(z_n):\n",
    "    axes[1].axhline(z_n / H_DOOR, color=\"C3\", ls=\"--\", lw=1)\n",
    "axes[1].set_xlabel(r\"temperature ($^\\circ$C)\")\n",
    "axes[1].legend()\n",
    "\n",
    "fig.suptitle(\"Steckler doorway profiles - dimensional reconstruction by similarity\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a7f2317",
   "metadata": {},
   "source": [
    "## Centreline plane - plume and stratification\n",
    "\n",
    "The room mid-plane through the burner and the door centre (the room centre\n",
    "$x = 0$ m, plane normal $x$) shows the buoyant plume, the hot upper layer and\n",
    "the exchange through the doorway into the ambient plenum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f8285dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_centreline_plane(cfg):\n",
    "    plane = cfg.output.exports[\"plane_series\"].series.planes[\"centreline\"]\n",
    "    pts = pd.read_csv(plane.points_filename)\n",
    "    ys = np.sort(pts[\"y\"].unique())\n",
    "    zs = np.sort(pts[\"z\"].unique())\n",
    "    iy = {v: i for i, v in enumerate(ys)}\n",
    "    iz = {v: i for i, v in enumerate(zs)}\n",
    "    cols = [str(int(i)) for i in pts[\"idx\"]]\n",
    "    col_yi = np.array([iy[y] for y in pts[\"y\"]])\n",
    "    col_zi = np.array([iz[z] for z in pts[\"z\"]])\n",
    "\n",
    "    df = plane.read_full_data(\"temp\")\n",
    "    start = int(df[\"time_step\"].max() * 0.6)\n",
    "    df = df[df[\"time_step\"] >= start]\n",
    "    arr = df[cols].to_numpy(dtype=float).mean(axis=0)\n",
    "    grid = np.full((zs.size, ys.size), np.nan)\n",
    "    grid[col_zi, col_yi] = arr\n",
    "    return ys, zs, grid\n",
    "\n",
    "\n",
    "ys, zc, Tc = load_centreline_plane(cfg)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "pc = ax.pcolormesh(ys, zc, Tc, shading=\"auto\", cmap=\"inferno\")\n",
    "fig.colorbar(pc, ax=ax, label=r\"$T$ (temperature)\")\n",
    "ax.set_xlabel(\"y (lattice) - plenum -> doorway -> room interior -> back wall\")\n",
    "ax.set_ylabel(\"z (lattice)\")\n",
    "ax.set_title(\"Time-mean temperature, room centreline (x = 0 m)\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5474b0c8",
   "metadata": {},
   "source": [
    "## Centreline plane - instantaneous velocity\n",
    "\n",
    "The same room mid-plane ($x = 0$ m), but the **last instantaneous** velocity\n",
    "snapshot rather than the time mean: the buoyant plume rising off the burner, the\n",
    "ceiling jet, and the two-way exchange through the doorway into the plenum. The\n",
    "colour is the velocity magnitude; the white unit arrows show the in-plane\n",
    "direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2981f192",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_centreline_velocity_last(cfg):\n",
    "    \"\"\"Last-snapshot velocity components on the room centreline plane (x = 0).\n",
    "\n",
    "    Unlike the temperature map above (time-mean), this returns the single most\n",
    "    recent instantaneous frame so the plume and doorway jet show their live\n",
    "    structure rather than a smoothed average.\n",
    "    \"\"\"\n",
    "    plane = cfg.output.exports[\"plane_series\"].series.planes[\"centreline\"]\n",
    "    pts = pd.read_csv(plane.points_filename)\n",
    "    ys = np.sort(pts[\"y\"].unique())\n",
    "    zs = np.sort(pts[\"z\"].unique())\n",
    "    iy = {v: i for i, v in enumerate(ys)}\n",
    "    iz = {v: i for i, v in enumerate(zs)}\n",
    "    cols = [str(int(i)) for i in pts[\"idx\"]]\n",
    "    col_yi = np.array([iy[y] for y in pts[\"y\"]])\n",
    "    col_zi = np.array([iz[z] for z in pts[\"z\"]])\n",
    "\n",
    "    comp = {}\n",
    "    step = 0\n",
    "    for macr in (\"ux\", \"uy\", \"uz\"):\n",
    "        df = plane.read_full_data(macr)\n",
    "        step = int(df[\"time_step\"].max())\n",
    "        row = df.loc[df[\"time_step\"].idxmax(), cols].to_numpy(dtype=float)\n",
    "        grid = np.full((zs.size, ys.size), np.nan)\n",
    "        grid[col_zi, col_yi] = row\n",
    "        comp[macr] = grid\n",
    "    return ys, zs, comp[\"ux\"], comp[\"uy\"], comp[\"uz\"], step\n",
    "\n",
    "\n",
    "ys, zs, ux, uy, uz, step = load_centreline_velocity_last(cfg)\n",
    "speed = np.sqrt(ux**2 + uy**2 + uz**2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "pc = ax.pcolormesh(ys, zs, speed, shading=\"auto\", cmap=\"viridis\")\n",
    "fig.colorbar(pc, ax=ax, label=r\"$|u|$ (lattice)\")\n",
    "\n",
    "# In-plane velocity (uy horizontal, uz vertical), coarsely subsampled and scaled\n",
    "# by magnitude so the plume and doorway jet stand out and the quiescent ambient\n",
    "# stays clean. Masked solid/exterior nodes are zeroed (no arrow).\n",
    "sy = max(1, ys.size // 32)\n",
    "sz = max(1, zs.size // 24)\n",
    "Yq, Zq = np.meshgrid(ys[::sy], zs[::sz])\n",
    "ax.quiver(\n",
    "    Yq, Zq, np.nan_to_num(uy[::sz, ::sy]), np.nan_to_num(uz[::sz, ::sy]), color=\"w\", width=0.003\n",
    ")\n",
    "\n",
    "ax.set_xlabel(\"y (lattice) - plenum -> doorway -> room interior -> back wall\")\n",
    "ax.set_ylabel(\"z (lattice)\")\n",
    "ax.set_title(f\"Instantaneous velocity magnitude, room centreline (x = 0), step {step}\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11a85e93",
   "metadata": {},
   "source": [
    "## Version and configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0195df2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Version / configuration, read from the run's info.yaml (same as the other\n",
    "# validation notebooks).\n",
    "sim_info = cfg.output.read_info()\n",
    "print(\"nassu version:\", sim_info[\"version\"])\n",
    "print(\"commit hash:\", sim_info[\"commit\"])\n",
    "print(\"case:\", case_path.relative_to(project_root))\n",
    "print(\"domain:\", cfg.domain.domain_size)\n",
    "\n",
    "lm = cfg.models.low_mach\n",
    "nu = (cfg.models.LBM.tau - 0.5) / 3.0\n",
    "rho_inf = lm.rho_inf if lm.rho_inf is not None else lm.P_thermo / (lm.r * lm.T_ref)\n",
    "print(\"operator:\", cfg.models.LBM.coll_oper, \" bulk_viscosity:\", cfg.models.LBM.bulk_viscosity)\n",
    "print(\"tau (fluid):\", cfg.models.LBM.tau, \" nu:\", nu)\n",
    "print(\"EOS: r =\", lm.r, \" P_thermo =\", lm.P_thermo, \" T_ref =\", lm.T_ref, \" rho_inf =\", rho_inf)\n",
    "print(\"Pr:\", lm.Pr, \" thermal diffusivity alpha = nu/Pr:\", nu / lm.Pr)\n",
    "print(\"gravity:\", lm.gravity)\n",
    "print(\"energy source rate:\", lm.energy.source_regions[0].rate)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "nassu",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}