{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1644473e-ab1e-45dc-b601-228e06168a3b",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Primary Cooling System Coastdown\n",
    "\n",
    "In this example we wish to take a look at a research nuclear reactor's Primary Cooling System (PCS) which begins to shut down (say, in a Loss Of Flow event). Luckily, a flywheel is installed, which dampens the quick loss of flow and saves the day (supposedly).\n",
    "\n",
    "This example does not deal with any specific real reactor parameters, and includes only 3 effective components:\n",
    "\n",
    "\n",
    "1. Pump\n",
    "1. Flywheel\n",
    "1. Resistor, effectively representing PCS + core\n",
    "\n",
    "We will solve the problem both numerically and analytically and compare the two solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81aaaab7-ce60-4b12-919b-ad003752c47a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from matplotlib import rc\n",
    "\n",
    "rc(\"font\", family=\"serif\", size=15)\n",
    "rc(\"savefig\", dpi=600)\n",
    "rc(\"figure\", figsize=(8, 6))\n",
    "rc(\"mathtext\", fontset=\"dejavuserif\")\n",
    "rc(\"lines\", linewidth=3)\n",
    "\n",
    "from stream.calculations import Inertia, Pump\n",
    "from stream.composition import ResistorFromKnownPoint"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4ff2223-5f99-401c-9695-e22d99ec1c82",
   "metadata": {},
   "source": [
    "## Components\n",
    "#### Pump"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d18ce477-f3a3-456a-90b6-d6605380d781",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dp0 = 1.6e5\n",
    "pump = Pump(pressure=dp0, name=\"Pump\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2c1ef6b-162e-4260-ae9f-621c794be857",
   "metadata": {},
   "source": [
    "#### Flywheel (Electrically, an Inductor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44ea4524-5db7-402c-98f4-01d4886319ae",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "inertia = 8e3\n",
    "flywheel = Inertia(inertia=inertia, name=\"Flywheel\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22865d1b-c909-4273-aaad-cb94a755f41d",
   "metadata": {},
   "source": [
    "#### Resistor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7fa7abb-d79b-4ebc-a94c-354bcd824209",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from stream.substances import light_water\n",
    "from stream.units import hour\n",
    "\n",
    "T = 20.0\n",
    "Q0 = 2000 / hour\n",
    "rho0 = light_water.density(T)\n",
    "mdot0 = Q0 * rho0\n",
    "\n",
    "resistor = ResistorFromKnownPoint(dp=-dp0, mdot=mdot0, behavior=\"parabolic\", Tin=T, fluid=light_water, name=\"PCS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edc0874f-fb90-411b-b3e6-a6f95bf0dd0e",
   "metadata": {},
   "source": [
    "## Assembling a Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eebe54a0-0ff0-4e06-8e3e-e3231281d2a1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "from stream.calculations import KirchhoffWDerivatives\n",
    "from stream.calculations.kirchhoff import to_str\n",
    "from stream.composition import flow_edge, flow_graph, flow_graph_to_agr_and_k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7476eb9-1186-4357-93c6-71e994cee827",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fg = flow_graph(flow_edge((\"A\", \"B\"), pump, flywheel), flow_edge((\"B\", \"A\"), resistor))\n",
    "\n",
    "agr, K = flow_graph_to_agr_and_k(\n",
    "    fg,\n",
    "    inertial_comps=[flywheel],\n",
    "    k_constructor=KirchhoffWDerivatives,\n",
    "    funcs={resistor: dict(Tin=T)},\n",
    ")\n",
    "from stream.analysis.report import report\n",
    "\n",
    "report(agr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3da027b9-6482-49e9-b04b-44a48d7bb1fe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "steady = agr.save(\n",
    "    agr.solve_steady(\n",
    "        {\n",
    "            K.name: {\n",
    "                \"(A -> B, 0)\": mdot0,\n",
    "                \"(A -> B, mdot2 0)\": 0.0,\n",
    "                \"(B -> A, 0)\": mdot0,\n",
    "                \"(B -> A, mdot2 0)\": 0.0,\n",
    "            },\n",
    "            pump.name: dict(Tin=T, pressure=dp0),\n",
    "            resistor.name: dict(Tin=T, pressure=-dp0),\n",
    "            flywheel.name: dict(Tin=T, pressure=0.0),\n",
    "        }\n",
    "    )\n",
    ")\n",
    "\n",
    "pd.DataFrame(steady.filter_values(lambda v: not np.isclose(v, 0)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bae8f04b-50d6-451b-9025-25fea9b804e3",
   "metadata": {},
   "source": [
    "## Analytic Solution of the Transient\n",
    "\n",
    "Assuming the pump has shut down immediately (that is, imposes zero pressure drop), the Kirchhoff Voltage Law (KVL) of the system looks like:\n",
    "\n",
    "$$ L\\frac{d\\dot{m}}{dt} = - k \\dot{m} ^ 2 $$\n",
    "\n",
    "Using $\\alpha \\equiv k/L$:\n",
    "\n",
    "$$ \\frac{d\\dot{m}}{dt} = -\\alpha \\dot{m} ^ 2 $$\n",
    "\n",
    "It is easy to verify that the equation above is solved by:\n",
    "\n",
    "$$ \\dot{m}(t) = \\frac{\\dot{m}_0}{1 + \\dot{m}_0\\alpha t} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8959c7d1-5be7-4c0f-96e3-d40aff0d30eb",
   "metadata": {},
   "source": [
    "## Transient Calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb23c0c4-5770-4ed1-b8ef-ada6eb2a53ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "pump.p = 0.0\n",
    "t = np.linspace(0, 300, 100)\n",
    "transient = agr.solve(steady, time=t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff365bdf-af9f-45f5-8ced-ec35c8e6a63a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(t, -agr.at_times(transient, resistor, \"pressure\") / 1e5, label=\"resistor\")\n",
    "plt.plot(t, agr.at_times(transient, flywheel, \"pressure\") / 1e5, label=\"flywheel\")\n",
    "plt.xlabel(\"Time [s]\")\n",
    "plt.ylabel(r\"$\\Delta p$ [Bar]\")\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bb6384b-ef54-4120-9555-7760946ad8ed",
   "metadata": {},
   "source": [
    "Comparing analytic and calculation of $\\dot{m}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61edd27e-11ce-4614-9d58-8783e0417b0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mdotc = agr.at_times(transient, K, to_str((\"A\", \"B\", 0)))\n",
    "mdot0 = steady[K.name][to_str((\"A\", \"B\", 0))]\n",
    "α = np.abs(resistor.dp_out(mdot=1.0, Tin=T) / inertia)\n",
    "mdota = mdot0 / (1 + α * mdot0 * t)\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(t, mdotc / mdot0, label=\"Calculation\")\n",
    "plt.plot(t, mdota / mdot0, label=\"Analytic\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"Time [s]\")\n",
    "plt.ylabel(r\"$\\dot{m} / \\dot{m}_0$ [Kg/s]\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59ec8be8-982f-4995-bd57-64fb589232ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(t, 100 * (mdotc - mdota) / mdota)\n",
    "plt.xlabel(\"Time [s]\")\n",
    "plt.ylabel(\"Relative Error [%]\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86ae71ae-ead7-4b50-8019-bc0d4c6e6667",
   "metadata": {},
   "source": [
    "As one can see, the analytic solution is well reproduced by the calculation."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}