{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Closed-domain thermodynamic-pressure rise (sealed heated box)\n",
    "\n",
    "**Results pending GPU execution.** This notebook is a scaffold: the analysis is\n",
    "complete but it has not yet been run against simulation output. Run the case with\n",
    "`uv run nassu run validation/buoyant_flows/05_closed_domain_pressure/05_closed_domain_pressure.nassu.yaml`\n",
    "on a GPU, then execute this notebook.\n",
    "\n",
    "This case validates the **closed-domain** low-Mach thermodynamic-pressure closure\n",
    "(`models.low_mach.domain_closure: closed`, Taha et al. 2024). A sealed, rigid,\n",
    "adiabatic box of ideal gas is heated at a known uniform volumetric rate `Q`. With no\n",
    "inlet/outlet the total mass is conserved, so the closure evolves the spatially-uniform\n",
    "thermodynamic pressure from total-mass conservation `P(t) = M r / I(t)`, with\n",
    "`I(t) = integral (1/T) dV`.\n",
    "\n",
    "The analytic target is the constant-volume first-law result for a sealed rigid\n",
    "adiabatic ideal gas heated at total rate `Q_total`:\n",
    "\n",
    "$$ \\frac{dP}{dt} = (\\gamma - 1)\\,\\frac{Q_\\text{total}}{V}, $$\n",
    "\n",
    "exact for adiabatic walls (the conduction term integrates to the zero wall heat flux).\n",
    "With whole-box uniform heating the volumetric rate `Q` gives `Q_total = Q V`, so the\n",
    "target reduces to `dP/dt = (gamma - 1) Q`. The notebook recovers `P(t)` node-wise from\n",
    "the EOS, fits the slope, and asserts it matches the target to within 2%."
   ],
   "id": "cell-0"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pathlib\n",
    "import glob\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import h5py\n",
    "\n",
    "from nassu.cfg.model import ConfigScheme\n",
    "\n",
    "import nassu.viz as common\n",
    "\n",
    "common.use_style()"
   ],
   "id": "cell-1"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "# Run from the project root: the simulation `save_path` (and every output path read\n",
    "# off the parsed config below) is repo-root-relative.\n",
    "project_root = _find_project_root()\n",
    "os.chdir(project_root)\n",
    "case_dir = project_root / \"validation/buoyant_flows/05_closed_domain_pressure\"\n",
    "config_path = case_dir / \"05_closed_domain_pressure.nassu.yaml\"\n",
    "\n",
    "sim_cfg = ConfigScheme.sim_cfgs_from_file(str(config_path))[0]\n",
    "n_steps = sim_cfg.n_steps\n",
    "print(f\"loaded {sim_cfg.name}: n_steps={n_steps}, domain={sim_cfg.domain.domain_size}\")"
   ],
   "id": "cell-2"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Analytic quantities read straight off the parsed low-Mach config, so the target is\n",
    "# never hard-coded: gamma = cp / cv with cv = cp - r, Q is the whole-box volumetric\n",
    "# source rate, and the target slope is (gamma - 1) Q.\n",
    "lm = sim_cfg.models.low_mach\n",
    "r_gas = lm.r\n",
    "cp = lm.cp\n",
    "cv = cp - r_gas\n",
    "gamma = cp / cv\n",
    "P0 = lm.P_thermo\n",
    "Q = lm.energy.source_region_rate  # volumetric heat-release rate (energy/vol/step)\n",
    "\n",
    "ds = sim_cfg.domain.domain_size\n",
    "V = ds.x * ds.y * ds.z\n",
    "Q_total = Q * V  # whole-box uniform heating\n",
    "\n",
    "dPdt_analytic = (gamma - 1.0) * Q_total / V  # = (gamma - 1) Q\n",
    "print(f\"r={r_gas}, cp={cp}, cv={cv}, gamma={gamma}\")\n",
    "print(f\"V={V}, Q={Q:.3e}, Q_total={Q_total:.3e}\")\n",
    "print(f\"analytic dP/dt = (gamma-1) Q_total/V = {dPdt_analytic:.6e} per step\")"
   ],
   "id": "cell-3"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Locate the instantaneous volume export off the parsed config and read each snapshot.\n",
    "# The volume export stores the full domain as a single block per time group:\n",
    "# `t<step>/block0/<macr>`, each a 3-D array in (z, y, x) layout.\n",
    "export = sim_cfg.output.instantaneous[\"default\"]\n",
    "h5_path = sorted(glob.glob(str(export.base_path) + \"*.h5\"))[0]\n",
    "print(\"reading\", h5_path)\n",
    "\n",
    "\n",
    "def read_snapshots(path):\n",
    "    with h5py.File(path, \"r\") as h:\n",
    "        times = sorted({k.split(\"/\")[0] for k in h.keys()}, key=lambda s: float(s[1:]))\n",
    "        steps = np.array([float(t[1:]) for t in times])\n",
    "        data = {}\n",
    "        for t in times:\n",
    "            g = h[t][\"block0\"]\n",
    "            data[t] = {m: g[m][:] for m in (\"rho\", \"ux\", \"uy\", \"uz\", \"temp\")}\n",
    "    return steps, times, data\n",
    "\n",
    "\n",
    "steps, times, snaps = read_snapshots(h5_path)\n",
    "print(f\"{len(steps)} snapshots, steps {steps[0]:.0f} .. {steps[-1]:.0f}\")"
   ],
   "id": "cell-4"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recover the spatially-uniform thermodynamic pressure node-wise from the EOS\n",
    "# rho = P/(r T)  =>  P = r * rho * T, exact since P is uniform. Also collect the peak\n",
    "# speed (must stay quiescent) and the spatial spread of rho/T (must stay uniform).\n",
    "cs = 1.0 / np.sqrt(3.0)\n",
    "P_t, umax, rho_spread, T_spread = [], [], [], []\n",
    "any_nan = False\n",
    "for t in times:\n",
    "    s = snaps[t]\n",
    "    rho, T = s[\"rho\"], s[\"temp\"]\n",
    "    P_t.append(float(r_gas * np.nanmean(rho * T)))\n",
    "    speed = np.sqrt(s[\"ux\"] ** 2 + s[\"uy\"] ** 2 + s[\"uz\"] ** 2)\n",
    "    umax.append(float(np.nanmax(speed)))\n",
    "    rho_spread.append(float(np.nanstd(rho) / np.nanmean(rho)))\n",
    "    T_spread.append(float(np.nanstd(T) / np.nanmean(T)))\n",
    "    any_nan = any_nan or bool(np.isnan(rho).any() or np.isnan(T).any())\n",
    "\n",
    "P_t = np.array(P_t)\n",
    "umax = np.array(umax)\n",
    "rho_spread = np.array(rho_spread)\n",
    "T_spread = np.array(T_spread)\n",
    "\n",
    "# Least-squares slope of P(t) over the whole run (uniform heating -> linear from step 0).\n",
    "slope, intercept = np.polyfit(steps, P_t, 1)\n",
    "rel_err = abs(slope - dPdt_analytic) / abs(dPdt_analytic)\n",
    "print(f\"P(0) recovered    = {P_t[0]:.6f}  (config P0 = {P0})\")\n",
    "print(f\"fitted   dP/dt    = {slope:.6e} per step\")\n",
    "print(f\"analytic dP/dt    = {dPdt_analytic:.6e} per step\")\n",
    "print(f\"relative error    = {rel_err * 100:.3f} %\")\n",
    "print(f\"peak |u| over run = {umax.max():.3e}  (Ma_max = {umax.max() / cs:.3e})\")"
   ],
   "id": "cell-5"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results"
   ],
   "id": "cell-6"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 3, figsize=(16, 4))\n",
    "\n",
    "# (1) P(t) recovered vs the analytic linear rise.\n",
    "P_analytic = P0 + dPdt_analytic * steps\n",
    "ax[0].plot(steps, P_t, \"o\", ms=4, label=r\"$P(t) = r\\,\\langle\\rho T\\rangle$ (recovered)\")\n",
    "ax[0].plot(steps, P_analytic, \"--\", c=\"k\", lw=1.5, label=r\"analytic $P_0 + (\\gamma-1)\\,Q\\,t$\")\n",
    "ax[0].set_xlabel(\"step\")\n",
    "ax[0].set_ylabel(\"thermodynamic pressure $P$\")\n",
    "ax[0].set_title(f\"Closed-domain pressure rise (err {rel_err * 100:.2f}%)\")\n",
    "ax[0].legend()\n",
    "ax[0].grid(alpha=0.3)\n",
    "\n",
    "# (2) Residual of the recovered pressure about the analytic rise, with the +/-2%\n",
    "# acceptance band scaled to the analytic rise at each step (absolute units).\n",
    "residual = P_t - P_analytic\n",
    "tol_band = 0.02 * dPdt_analytic * steps\n",
    "ax[1].plot(steps, residual, \"-o\", ms=4, c=\"C3\", label=r\"$P(t) - P_\\mathrm{analytic}$\")\n",
    "ax[1].fill_between(steps, -tol_band, tol_band, color=\"grey\", alpha=0.2, label=\"2% band\")\n",
    "ax[1].axhline(0.0, ls=\"--\", c=\"k\", lw=1)\n",
    "ax[1].set_xlabel(\"step\")\n",
    "ax[1].set_ylabel(r\"$P - P_\\mathrm{analytic}$\")\n",
    "ax[1].set_title(\"Pressure residual vs analytic rise\")\n",
    "ax[1].legend()\n",
    "ax[1].grid(alpha=0.3)\n",
    "\n",
    "# (3) Quiescence and spatial-uniformity diagnostics.\n",
    "ax[2].plot(steps, umax, \"-o\", ms=3, label=r\"$\\max|u|$\")\n",
    "ax[2].plot(steps, rho_spread, \"-s\", ms=3, label=r\"$\\sigma(\\rho)/\\langle\\rho\\rangle$\")\n",
    "ax[2].plot(steps, T_spread, \"-^\", ms=3, label=r\"$\\sigma(T)/\\langle T\\rangle$\")\n",
    "ax[2].axhline(0.1 * cs, ls=\"--\", c=\"grey\", lw=1, label=r\"$\\mathrm{Ma}=0.1$\")\n",
    "ax[2].set_xlabel(\"step\")\n",
    "ax[2].set_ylabel(\"value\")\n",
    "ax[2].set_title(\"Box stays quiescent and spatially uniform\")\n",
    "ax[2].legend()\n",
    "ax[2].grid(alpha=0.3)\n",
    "fig.tight_layout()"
   ],
   "id": "cell-7"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "Acceptance criteria for the closed-domain pressure closure:\n",
    "\n",
    "- **Pressure-rise rate** - the fitted `dP/dt` matches the analytic `(gamma - 1) Q`\n",
    "  to within 2%.\n",
    "- **No divergence** - all snapshots are finite (no NaN in `rho`/`T`).\n",
    "- **Quiescent box** - `max|u|` stays far below `Ma = 0.1` (no spurious dilatation\n",
    "  flow or reduced-pressure checkerboard mode).\n",
    "- **Spatial uniformity** - `rho` and `T` stay spatially uniform (small relative\n",
    "  spread), confirming the closure keeps the uniformly-heated sealed box homogeneous."
   ],
   "id": "cell-8"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert not any_nan, \"divergence: NaN found in the fields\"\n",
    "assert rel_err < 0.02, f\"dP/dt off analytic by {rel_err * 100:.2f}% (> 2% tolerance)\"\n",
    "assert umax.max() < 0.1 * cs, f\"velocity not subsonic: max|u|={umax.max():.3e}\"\n",
    "assert rho_spread.max() < 1e-2, f\"rho not spatially uniform: max spread {rho_spread.max():.3e}\"\n",
    "assert T_spread.max() < 1e-2, f\"T not spatially uniform: max spread {T_spread.max():.3e}\"\n",
    "print(\"VALIDATION PASSED: closed-domain dP/dt matches (gamma-1) Q within 2%.\")"
   ],
   "id": "cell-9"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Version"
   ],
   "id": "cell-10"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_info = sim_cfg.output.read_info()\n",
    "print(sim_info)"
   ],
   "id": "cell-11"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Configuration"
   ],
   "id": "cell-12"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Code\n",
    "\n",
    "Code(filename=str(config_path), language=\"yaml\")"
   ],
   "id": "cell-13"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}