diff --git a/notebooks/disc_cooling_sim.ipynb b/notebooks/disc_cooling_sim.ipynb
new file mode 100644
index 0000000000..e9c0537077
--- /dev/null
+++ b/notebooks/disc_cooling_sim.ipynb
@@ -0,0 +1,712 @@
+{
+ "nbformat": 4,
+ "nbformat_minor": 5,
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.10.0"
+  },
+  "colab": {
+   "provenance": [],
+   "gpuType": "T4"
+  },
+  "accelerator": "GPU"
+ },
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "badge-cell",
+   "metadata": {},
+   "source": [
+    "# Disc Injection-Moulding Cooling Simulation\n",
+    "\n",
+    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/Tuesdaythe13th/warp/blob/claude/disc-config-struct-ndXr3/notebooks/disc_cooling_sim.ipynb)\n",
+    "\n",
+    "A GPU-accelerated physics simulation of disc cooling after injection moulding, using [NVIDIA Warp](https://github.com/NVIDIA/warp).\n",
+    "\n",
+    "**What this notebook covers:**\n",
+    "- 2-D axisymmetric heat diffusion in a polymer disc (cylindrical coordinates)\n",
+    "- Avrami-based crystallinity evolution (PET + CSR + PTFE blend)\n",
+    "- Warp-risk scoring: thermal gradient × crystallinity asymmetry → predicted deflection\n",
+    "- Radial plots and 2-D field visualisations\n",
+    "\n",
+    "**Physics upgrade path (beyond this notebook):**\n",
+    "1. Replace the placeholder Avrami kinetics with a Nakamura model fitted to DSC data for your actual formulation.\n",
+    "2. Add a thermoelastic plate / eigenstrain solver so `warp_risk` becomes out-of-plane deflection in millimetres."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "install-header",
+   "metadata": {},
+   "source": [
+    "## 1 — Install dependencies"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "install-cell",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "!pip install -q warp-lang matplotlib"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "imports-header",
+   "metadata": {},
+   "source": [
+    "## 2 — Imports and Warp initialisation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "imports-cell",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import math\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib.colors as mcolors\n",
+    "import warp as wp\n",
+    "\n",
+    "wp.config.quiet = True\n",
+    "wp.init()\n",
+    "\n",
+    "DEVICE = \"cuda\" if wp.get_cuda_device_count() > 0 else \"cpu\"\n",
+    "print(f\"Running on: {DEVICE}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "structs-header",
+   "metadata": {},
+   "source": [
+    "## 3 — Warp structs\n",
+    "\n",
+    "Two structs are used to keep kernel signatures short:\n",
+    "\n",
+    "| Struct | Purpose |\n",
+    "|---|---|\n",
+    "| `DiscConfig` | Geometry and material properties (fixed for a given disc) |\n",
+    "| `CoolingParams` | Process parameters (vary between simulation runs) |"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "structs-cell",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@wp.struct\n",
+    "class DiscConfig:\n",
+    "    \"\"\"Geometry and material properties of the disc.\n",
+    "\n",
+    "    Attributes:\n",
+    "        radius:    Outer radius (m).\n",
+    "        thickness: Disc thickness (m).\n",
+    "        dr:        Radial grid spacing (m).\n",
+    "        dz:        Axial grid spacing (m).\n",
+    "        k:         Thermal conductivity (W/m·K).\n",
+    "        rho:       Density (kg/m³).\n",
+    "        cp:        Specific heat capacity (J/kg·K).\n",
+    "        T_g:       Glass-transition temperature (K).\n",
+    "        T_m:       Melting temperature (K).\n",
+    "        nx:        Number of radial grid cells.\n",
+    "        nz:        Number of axial grid cells.\n",
+    "    \"\"\"\n",
+    "    radius:    float\n",
+    "    thickness: float\n",
+    "    dr:        float\n",
+    "    dz:        float\n",
+    "    k:         float\n",
+    "    rho:       float\n",
+    "    cp:        float\n",
+    "    T_g:       float\n",
+    "    T_m:       float\n",
+    "    nx:        int\n",
+    "    nz:        int\n",
+    "\n",
+    "\n",
+    "@wp.struct\n",
+    "class CoolingParams:\n",
+    "    \"\"\"Process parameters for one simulation run.\n",
+    "\n",
+    "    Attributes:\n",
+    "        T_init:         Initial melt temperature (K).\n",
+    "        T_mold:         Mould wall temperature (K).\n",
+    "        dt:             Time step (s).\n",
+    "        total_time:     Total simulation time (s).\n",
+    "        avrami_k0:      Avrami rate pre-factor (1/s).\n",
+    "        avrami_n:       Avrami exponent (dimensionless).\n",
+    "        chi_max:        Maximum achievable crystallinity (0–1).\n",
+    "        warp_temp_coeff: Weight of thermal gradient in warp-risk score.\n",
+    "        warp_chi_coeff:  Weight of crystallinity asymmetry in warp-risk score.\n",
+    "    \"\"\"\n",
+    "    T_init:          float\n",
+    "    T_mold:          float\n",
+    "    dt:              float\n",
+    "    total_time:      float\n",
+    "    avrami_k0:       float\n",
+    "    avrami_n:        float\n",
+    "    chi_max:         float\n",
+    "    warp_temp_coeff: float\n",
+    "    warp_chi_coeff:  float"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "kernels-header",
+   "metadata": {},
+   "source": [
+    "## 4 — Warp kernels"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "kernels-cell",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@wp.func\n",
+    "def idx(i: int, j: int, nz: int) -> int:\n",
+    "    \"\"\"Row-major flat index for the (nx × nz) grid.\"\"\"\n",
+    "    return i * nz + j\n",
+    "\n",
+    "\n",
+    "@wp.func\n",
+    "def clamp_i(x: int, lo: int, hi: int) -> int:\n",
+    "    return wp.min(wp.max(x, lo), hi)\n",
+    "\n",
+    "\n",
+    "@wp.kernel\n",
+    "def init_temperature(\n",
+    "    T:      wp.array(dtype=wp.float32),\n",
+    "    params: CoolingParams,\n",
+    "):\n",
+    "    tid = wp.tid()\n",
+    "    T[tid] = wp.float32(params.T_init)\n",
+    "\n",
+    "\n",
+    "@wp.kernel\n",
+    "def init_scalar(\n",
+    "    a:     wp.array(dtype=wp.float32),\n",
+    "    value: float,\n",
+    "):\n",
+    "    tid = wp.tid()\n",
+    "    a[tid] = wp.float32(value)\n",
+    "\n",
+    "\n",
+    "@wp.kernel\n",
+    "def step_temperature(\n",
+    "    T_in:   wp.array(dtype=wp.float32),\n",
+    "    T_out:  wp.array(dtype=wp.float32),\n",
+    "    config: DiscConfig,\n",
+    "    params: CoolingParams,\n",
+    "):\n",
+    "    \"\"\"Explicit finite-difference heat diffusion in cylindrical coordinates.\n",
+    "\n",
+    "    Solves  ∂T/∂t = α (∂²T/∂r² + (1/r)∂T/∂r + ∂²T/∂z²)\n",
+    "    with Dirichlet mould-wall BCs on the top and bottom (z) faces and\n",
+    "    Neumann (zero-flux) BCs on the axis (r=0) and outer radius.\n",
+    "    \"\"\"\n",
+    "    tid = wp.tid()\n",
+    "    i = tid // config.nz\n",
+    "    j = tid - i * config.nz\n",
+    "\n",
+    "    alpha = config.k / (config.rho * config.cp)\n",
+    "    dr    = config.dr\n",
+    "    dz    = config.dz\n",
+    "    r     = (float(i) + 0.5) * dr\n",
+    "\n",
+    "    # Dirichlet BC on top/bottom mould walls.\n",
+    "    if j == 0 or j == config.nz - 1:\n",
+    "        T_out[tid] = wp.float32(params.T_mold)\n",
+    "        return\n",
+    "\n",
+    "    # Neumann BC: mirror stencil at axis and outer edge.\n",
+    "    im = clamp_i(i - 1, 0, config.nx - 1)\n",
+    "    ip = clamp_i(i + 1, 0, config.nx - 1)\n",
+    "    if i == 0:\n",
+    "        im = 1\n",
+    "    if i == config.nx - 1:\n",
+    "        ip = config.nx - 2\n",
+    "\n",
+    "    jm = j - 1\n",
+    "    jp = j + 1\n",
+    "\n",
+    "    Tc  = T_in[idx(i,  j,  config.nz)]\n",
+    "    Trm = T_in[idx(im, j,  config.nz)]\n",
+    "    Trp = T_in[idx(ip, j,  config.nz)]\n",
+    "    Tzm = T_in[idx(i,  jm, config.nz)]\n",
+    "    Tzp = T_in[idx(i,  jp, config.nz)]\n",
+    "\n",
+    "    d2Tdr2       = (Trp - 2.0 * Tc + Trm) / (dr * dr)\n",
+    "    dTdr_over_r  = 0.0\n",
+    "    if i > 0:\n",
+    "        dTdr_over_r = (Trp - Trm) / (2.0 * dr * r)\n",
+    "    d2Tdz2 = (Tzp - 2.0 * Tc + Tzm) / (dz * dz)\n",
+    "\n",
+    "    lap   = d2Tdr2 + dTdr_over_r + d2Tdz2\n",
+    "    Tnew  = Tc + params.dt * alpha * lap\n",
+    "    T_out[tid] = wp.float32(Tnew)\n",
+    "\n",
+    "\n",
+    "@wp.kernel\n",
+    "def update_crystallinity(\n",
+    "    T:       wp.array(dtype=wp.float32),\n",
+    "    chi_in:  wp.array(dtype=wp.float32),\n",
+    "    chi_out: wp.array(dtype=wp.float32),\n",
+    "    config:  DiscConfig,\n",
+    "    params:  CoolingParams,\n",
+    "):\n",
+    "    \"\"\"Simple Avrami-style crystallisation kinetics.\n",
+    "\n",
+    "    Crystal growth is fastest mid-way between T_g and T_m and saturates\n",
+    "    at chi_max.  Replace with a Nakamura model for production use.\n",
+    "    \"\"\"\n",
+    "    tid  = wp.tid()\n",
+    "    temp = T[tid]\n",
+    "    chi  = chi_in[tid]\n",
+    "\n",
+    "    if temp > config.T_g and temp < config.T_m:\n",
+    "        x    = (temp - config.T_g) / (config.T_m - config.T_g)\n",
+    "        x    = wp.max(0.0, wp.min(1.0, x))\n",
+    "        rate = params.avrami_k0 * x * (1.0 - chi / params.chi_max)\n",
+    "        chi  = chi + params.dt * params.avrami_n * rate\n",
+    "        chi  = wp.max(0.0, wp.min(params.chi_max, chi))\n",
+    "\n",
+    "    chi_out[tid] = wp.float32(chi)\n",
+    "\n",
+    "\n",
+    "@wp.kernel\n",
+    "def compute_warp_risk(\n",
+    "    T:         wp.array(dtype=wp.float32),\n",
+    "    chi:       wp.array(dtype=wp.float32),\n",
+    "    warp_risk: wp.array(dtype=wp.float32),\n",
+    "    config:    DiscConfig,\n",
+    "    params:    CoolingParams,\n",
+    "):\n",
+    "    \"\"\"Score each radial position by thermal gradient and crystallinity asymmetry.\n",
+    "\n",
+    "    Only threads at the mid-plane (j == nz/2) write to warp_risk[i].\n",
+    "    \"\"\"\n",
+    "    tid = wp.tid()\n",
+    "    i   = tid // config.nz\n",
+    "    j   = tid - i * config.nz\n",
+    "\n",
+    "    if j != config.nz // 2:\n",
+    "        return\n",
+    "\n",
+    "    top = T[idx(i, 0,              config.nz)]\n",
+    "    bot = T[idx(i, config.nz - 1,  config.nz)]\n",
+    "    mid = T[idx(i, j,              config.nz)]\n",
+    "\n",
+    "    chi_top = chi[idx(i, 1,              config.nz)]\n",
+    "    chi_bot = chi[idx(i, config.nz - 2,  config.nz)]\n",
+    "\n",
+    "    dT_thickness = wp.abs(top - bot)\n",
+    "    dT_mid       = wp.abs(mid - 0.5 * (top + bot))\n",
+    "    dchi         = wp.abs(chi_top - chi_bot)\n",
+    "\n",
+    "    risk = (\n",
+    "        params.warp_temp_coeff * (dT_thickness + dT_mid)\n",
+    "        + params.warp_chi_coeff * dchi\n",
+    "    )\n",
+    "    warp_risk[i] = wp.float32(risk)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "sim-class-header",
+   "metadata": {},
+   "source": [
+    "## 5 — Simulation class"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "sim-class-cell",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class ArtifexCoolingSim:\n",
+    "    \"\"\"Manages GPU arrays and orchestrates the time-stepping loop.\"\"\"\n",
+    "\n",
+    "    def __init__(self, config: DiscConfig, device: str = \"cpu\"):\n",
+    "        self.config = config\n",
+    "        self.device = device\n",
+    "        self.n      = config.nx * config.nz\n",
+    "\n",
+    "        self.T_a       = wp.zeros(self.n,      dtype=wp.float32, device=device)\n",
+    "        self.T_b       = wp.zeros(self.n,      dtype=wp.float32, device=device)\n",
+    "        self.chi_a     = wp.zeros(self.n,      dtype=wp.float32, device=device)\n",
+    "        self.chi_b     = wp.zeros(self.n,      dtype=wp.float32, device=device)\n",
+    "        self.warp_risk = wp.zeros(config.nx,   dtype=wp.float32, device=device)\n",
+    "\n",
+    "    def simulate_cooling(self, params: CoolingParams) -> dict:\n",
+    "        \"\"\"Run the full cooling simulation and return result metrics.\n",
+    "\n",
+    "        Args:\n",
+    "            params: Process parameters for this run.\n",
+    "\n",
+    "        Returns:\n",
+    "            Dictionary with scalar metrics and 2-D numpy field arrays.\n",
+    "        \"\"\"\n",
+    "        wp.launch(init_temperature, dim=self.n, inputs=[self.T_a, params],    device=self.device)\n",
+    "        wp.launch(init_scalar,      dim=self.n, inputs=[self.chi_a, 0.0],    device=self.device)\n",
+    "\n",
+    "        steps = int(params.total_time / params.dt)\n",
+    "        for _ in range(steps):\n",
+    "            wp.launch(\n",
+    "                kernel=step_temperature,\n",
+    "                dim=self.n,\n",
+    "                inputs=[self.T_a, self.T_b, self.config, params],\n",
+    "                device=self.device,\n",
+    "            )\n",
+    "            self.T_a, self.T_b = self.T_b, self.T_a\n",
+    "\n",
+    "            wp.launch(\n",
+    "                kernel=update_crystallinity,\n",
+    "                dim=self.n,\n",
+    "                inputs=[self.T_a, self.chi_a, self.chi_b, self.config, params],\n",
+    "                device=self.device,\n",
+    "            )\n",
+    "            self.chi_a, self.chi_b = self.chi_b, self.chi_a\n",
+    "\n",
+    "        wp.launch(\n",
+    "            kernel=compute_warp_risk,\n",
+    "            dim=self.n,\n",
+    "            inputs=[self.T_a, self.chi_a, self.warp_risk, self.config, params],\n",
+    "            device=self.device,\n",
+    "        )\n",
+    "\n",
+    "        T_np    = self.T_a.numpy().reshape((self.config.nx, self.config.nz))\n",
+    "        chi_np  = self.chi_a.numpy().reshape((self.config.nx, self.config.nz))\n",
+    "        risk_np = self.warp_risk.numpy()\n",
+    "\n",
+    "        max_delta_t      = float(np.max(np.abs(T_np[:, 0] - T_np[:, -1])))\n",
+    "        avg_chi_groove   = float(np.mean(chi_np[:, 1]))\n",
+    "        max_warp_risk    = float(np.max(risk_np))\n",
+    "        is_ok = avg_chi_groove < 0.15 and max_warp_risk < 15.0\n",
+    "\n",
+    "        return {\n",
+    "            \"max_delta_t\":       max_delta_t,\n",
+    "            \"avg_chi_groove\":    avg_chi_groove,\n",
+    "            \"max_warp_risk\":     max_warp_risk,\n",
+    "            \"is_ok\":             is_ok,\n",
+    "            \"final_T_field\":     T_np,\n",
+    "            \"final_chi_field\":   chi_np,\n",
+    "            \"warp_risk_radial\":  risk_np,\n",
+    "        }"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "params-header",
+   "metadata": {},
+   "source": [
+    "## 6 — Set up geometry and process parameters\n",
+    "\n",
+    "Edit the cells below to match your disc geometry and material properties.\n",
+    "The default values approximate a 120 mm PET optical disc cooled in a 25 °C mould."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "params-cell",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ---------------------------------------------------------------------------\n",
+    "# Geometry\n",
+    "# ---------------------------------------------------------------------------\n",
+    "RADIUS_MM    = 60.0   # outer radius (mm)\n",
+    "THICKNESS_MM = 1.2    # disc thickness (mm)\n",
+    "NX           = 60     # radial grid cells\n",
+    "NZ           = 20     # axial grid cells (must be ≥ 4)\n",
+    "\n",
+    "radius    = RADIUS_MM    * 1e-3\n",
+    "thickness = THICKNESS_MM * 1e-3\n",
+    "dr        = radius    / NX\n",
+    "dz        = thickness / NZ\n",
+    "\n",
+    "# ---------------------------------------------------------------------------\n",
+    "# Material (PET + 10 % CSR + 5 % PTFE blend, approximate values)\n",
+    "# ---------------------------------------------------------------------------\n",
+    "config        = DiscConfig()\n",
+    "config.radius    = radius\n",
+    "config.thickness = thickness\n",
+    "config.dr        = dr\n",
+    "config.dz        = dz\n",
+    "config.k         = 0.29      # W/m·K\n",
+    "config.rho       = 1360.0    # kg/m³\n",
+    "config.cp        = 1250.0    # J/kg·K\n",
+    "config.T_g       = 348.0     # K  (~75 °C)\n",
+    "config.T_m       = 533.0     # K  (~260 °C)\n",
+    "config.nx        = NX\n",
+    "config.nz        = NZ\n",
+    "\n",
+    "# Fourier stability: dt < dr²/(2α) and dt < dz²/(2α)\n",
+    "alpha   = config.k / (config.rho * config.cp)\n",
+    "dt_max  = 0.4 * min(dr, dz) ** 2 / alpha\n",
+    "print(f\"Thermal diffusivity α = {alpha:.2e} m²/s\")\n",
+    "print(f\"Stability limit dt_max = {dt_max*1e3:.3f} ms  →  using dt = {dt_max*1e3*0.9:.3f} ms\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cooling-params-cell",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ---------------------------------------------------------------------------\n",
+    "# Process parameters — edit these to explore different moulding conditions\n",
+    "# ---------------------------------------------------------------------------\n",
+    "params                   = CoolingParams()\n",
+    "params.T_init            = 553.0          # K  (melt inlet, ~280 °C)\n",
+    "params.T_mold            = 298.0          # K  (mould wall, ~25 °C)\n",
+    "params.dt                = dt_max * 0.9   # s  (90 % of stability limit)\n",
+    "params.total_time        = 8.0            # s  (total cooling window)\n",
+    "params.avrami_k0         = 0.005          # 1/s  — placeholder; fit from DSC\n",
+    "params.avrami_n          = 2.5            # —     placeholder Avrami exponent\n",
+    "params.chi_max           = 0.35           # max crystallinity for this blend\n",
+    "params.warp_temp_coeff   = 0.2            # K⁻¹  — normalisation weight\n",
+    "params.warp_chi_coeff    = 50.0           # —     crystallinity weight\n",
+    "\n",
+    "steps = int(params.total_time / params.dt)\n",
+    "print(f\"Time steps: {steps:,}   (dt = {params.dt*1e3:.3f} ms, total = {params.total_time} s)\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "run-header",
+   "metadata": {},
+   "source": [
+    "## 7 — Run the simulation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "run-cell",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim     = ArtifexCoolingSim(config, device=DEVICE)\n",
+    "results = sim.simulate_cooling(params)\n",
+    "\n",
+    "print(\"\\n=== Results ===\")\n",
+    "print(f\"  Max top-bottom ΔT         : {results['max_delta_t']:.2f} K\")\n",
+    "print(f\"  Mean groove crystallinity : {results['avg_chi_groove']:.4f}\")\n",
+    "print(f\"  Max warp-risk score       : {results['max_warp_risk']:.3f}\")\n",
+    "print(f\"  Quality gate passed       : {results['is_ok']}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "viz-header",
+   "metadata": {},
+   "source": [
+    "## 8 — Visualise the results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "viz-2d-cell",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "T_field   = results[\"final_T_field\"]    # (nx, nz)\n",
+    "chi_field = results[\"final_chi_field\"]  # (nx, nz)\n",
+    "risk      = results[\"warp_risk_radial\"] # (nx,)\n",
+    "\n",
+    "r_mm = np.linspace(dr / 2, radius - dr / 2, NX) * 1e3   # mm\n",
+    "z_mm = np.linspace(dz / 2, thickness - dz / 2, NZ) * 1e3  # mm\n",
+    "R, Z = np.meshgrid(r_mm, z_mm, indexing=\"ij\")\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(16, 4))\n",
+    "\n",
+    "# ---- Temperature field ----\n",
+    "ax = axes[0]\n",
+    "im = ax.pcolormesh(R, Z, T_field, cmap=\"inferno\", shading=\"auto\")\n",
+    "fig.colorbar(im, ax=ax, label=\"Temperature (K)\")\n",
+    "ax.set_xlabel(\"Radius (mm)\")\n",
+    "ax.set_ylabel(\"Thickness (mm)\")\n",
+    "ax.set_title(f\"Final temperature field  (t = {params.total_time} s)\")\n",
+    "\n",
+    "# ---- Crystallinity field ----\n",
+    "ax = axes[1]\n",
+    "im = ax.pcolormesh(R, Z, chi_field, cmap=\"viridis\", shading=\"auto\",\n",
+    "                   vmin=0, vmax=params.chi_max)\n",
+    "fig.colorbar(im, ax=ax, label=\"Crystallinity χ\")\n",
+    "ax.set_xlabel(\"Radius (mm)\")\n",
+    "ax.set_ylabel(\"Thickness (mm)\")\n",
+    "ax.set_title(\"Final crystallinity field\")\n",
+    "\n",
+    "# ---- Warp-risk radial profile ----\n",
+    "ax = axes[2]\n",
+    "ax.plot(r_mm, risk, color=\"crimson\", linewidth=2)\n",
+    "ax.axhline(15.0, color=\"gray\", linestyle=\"--\", linewidth=1, label=\"Quality threshold\")\n",
+    "ax.fill_between(r_mm, risk, 15.0,\n",
+    "                where=(risk > 15.0), color=\"crimson\", alpha=0.2, label=\"At-risk zone\")\n",
+    "ax.set_xlabel(\"Radius (mm)\")\n",
+    "ax.set_ylabel(\"Warp-risk score\")\n",
+    "ax.set_title(\"Radial warp-risk profile\")\n",
+    "ax.legend(fontsize=8)\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "radial-header",
+   "metadata": {},
+   "source": [
+    "## 9 — Radial mid-plane profiles"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "radial-cell",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mid = NZ // 2\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
+    "\n",
+    "ax = axes[0]\n",
+    "ax.plot(r_mm, T_field[:, mid] - 273.15, label=\"Mid-plane\", color=\"tab:orange\")\n",
+    "ax.plot(r_mm, T_field[:, 0]  - 273.15, label=\"Top surface\", color=\"tab:blue\", linestyle=\"--\")\n",
+    "ax.plot(r_mm, T_field[:, -1] - 273.15, label=\"Bottom surface\", color=\"tab:green\", linestyle=\":\")\n",
+    "ax.set_xlabel(\"Radius (mm)\")\n",
+    "ax.set_ylabel(\"Temperature (°C)\")\n",
+    "ax.set_title(\"Radial temperature profiles\")\n",
+    "ax.legend()\n",
+    "ax.grid(alpha=0.3)\n",
+    "\n",
+    "ax = axes[1]\n",
+    "ax.plot(r_mm, chi_field[:, mid],  label=\"Mid-plane\",     color=\"tab:purple\")\n",
+    "ax.plot(r_mm, chi_field[:, 1],    label=\"Near top wall\", color=\"tab:blue\",  linestyle=\"--\")\n",
+    "ax.plot(r_mm, chi_field[:, -2],   label=\"Near bot wall\", color=\"tab:green\", linestyle=\":\")\n",
+    "ax.axhline(0.15, color=\"gray\", linestyle=\"--\", linewidth=1, label=\"χ threshold (0.15)\")\n",
+    "ax.set_xlabel(\"Radius (mm)\")\n",
+    "ax.set_ylabel(\"Crystallinity χ\")\n",
+    "ax.set_title(\"Radial crystallinity profiles\")\n",
+    "ax.legend(fontsize=8)\n",
+    "ax.grid(alpha=0.3)\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "sweep-header",
+   "metadata": {},
+   "source": [
+    "## 10 — Parameter sweep: mould temperature vs. warp risk\n",
+    "\n",
+    "Scan mould temperature from 15 °C to 80 °C and observe how the peak warp-risk score changes."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "sweep-cell",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "T_mold_range = np.linspace(288, 353, 12)  # 15 °C – 80 °C in K\n",
+    "sweep_results = []\n",
+    "\n",
+    "for T_mold_K in T_mold_range:\n",
+    "    p                = CoolingParams()\n",
+    "    p.T_init         = params.T_init\n",
+    "    p.T_mold         = float(T_mold_K)\n",
+    "    p.dt             = params.dt\n",
+    "    p.total_time     = params.total_time\n",
+    "    p.avrami_k0      = params.avrami_k0\n",
+    "    p.avrami_n       = params.avrami_n\n",
+    "    p.chi_max        = params.chi_max\n",
+    "    p.warp_temp_coeff = params.warp_temp_coeff\n",
+    "    p.warp_chi_coeff  = params.warp_chi_coeff\n",
+    "\n",
+    "    s   = ArtifexCoolingSim(config, device=DEVICE)\n",
+    "    res = s.simulate_cooling(p)\n",
+    "    sweep_results.append({\n",
+    "        \"T_mold_C\":    T_mold_K - 273.15,\n",
+    "        \"max_risk\":    res[\"max_warp_risk\"],\n",
+    "        \"avg_chi\":     res[\"avg_chi_groove\"],\n",
+    "        \"is_ok\":       res[\"is_ok\"],\n",
+    "    })\n",
+    "\n",
+    "T_mold_C = [r[\"T_mold_C\"]  for r in sweep_results]\n",
+    "max_risk  = [r[\"max_risk\"]  for r in sweep_results]\n",
+    "avg_chi   = [r[\"avg_chi\"]   for r in sweep_results]\n",
+    "colors    = [\"green\" if r[\"is_ok\"] else \"red\" for r in sweep_results]\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
+    "\n",
+    "ax = axes[0]\n",
+    "ax.bar(T_mold_C, max_risk, width=3.5, color=colors, edgecolor=\"black\", linewidth=0.5)\n",
+    "ax.axhline(15.0, color=\"gray\", linestyle=\"--\", label=\"Quality threshold\")\n",
+    "ax.set_xlabel(\"Mould temperature (°C)\")\n",
+    "ax.set_ylabel(\"Max warp-risk score\")\n",
+    "ax.set_title(\"Peak warp risk vs. mould temperature\")\n",
+    "ax.legend()\n",
+    "\n",
+    "ax = axes[1]\n",
+    "ax.bar(T_mold_C, avg_chi, width=3.5, color=colors, edgecolor=\"black\", linewidth=0.5)\n",
+    "ax.axhline(0.15, color=\"gray\", linestyle=\"--\", label=\"χ threshold\")\n",
+    "ax.set_xlabel(\"Mould temperature (°C)\")\n",
+    "ax.set_ylabel(\"Mean groove crystallinity χ\")\n",
+    "ax.set_title(\"Groove crystallinity vs. mould temperature\")\n",
+    "ax.legend()\n",
+    "\n",
+    "# Legend patch\n",
+    "from matplotlib.patches import Patch\n",
+    "legend_els = [Patch(facecolor=\"green\", label=\"Pass\"), Patch(facecolor=\"red\", label=\"Fail\")]\n",
+    "for ax in axes:\n",
+    "    ax.legend(handles=ax.get_legend().legend_handles + legend_els)\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "next-steps-header",
+   "metadata": {},
+   "source": [
+    "## Next steps\n",
+    "\n",
+    "| Upgrade | Why it matters |\n",
+    "|---|---|\n",
+    "| **Nakamura kinetics from DSC data** | Current Avrami constants are placeholders; real PET+CSR+PTFE kinetics differ substantially. |\n",
+    "| **Thermoelastic / eigenstrain plate solver** | Convert `warp_risk` from a dimensionless score into predicted out-of-plane deflection (mm). |\n",
+    "| **Separate top/bottom mould BCs** | Real injection tools often have asymmetric cooling channels. |\n",
+    "| **Isaac Sim USD export** | Map the final temperature and crystallinity fields onto a `wp.vec3` mesh for visualisation in Isaac Sim. |\n",
+    "| **Automatic time-step control** | Use an adaptive dt to skip the manual Fourier-stability calculation. |"
+   ]
+  }
+ ]
+}