{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MSRE Pump Transient Benchmark"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook simulates pump transients for the MSRE at zero power (10W). A dynamical system representing the MSRE is built using the [msrDynamics](https://github.com/LukeLabrie/msrDynamics) tool. The system is run at steady-state until $t = 2500$, where the fuel and coolant pumps are spun down at the rate defined by the ORNL data. After allowing the system to settle to a new equilibrium, at $t = 7500$, the pumps are spun up again. The model includes a reactivity response, `rho_control`, which models the response of the MSRE control rod to maintain constant power. It is modeled as a pure integrator of $\\frac{dn}{dt}$ from the point-kinetics equations with a coefficient of -1.0, i.e. it tries to exactly cancel reactivity changes introduced by flow changes, with the caveat that changes are limited to $\\pm$29 pcm/s, which was the maximum rate achievable by the MSRE control rods. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from parameters_U235_pump_transient import *\n",
    "from jitcdde import t\n",
    "from msrDynamics import Node, System\n",
    "import pandas as pd\n",
    "import os\n",
    "from scipy.interpolate import interp1d\n",
    "import sympy as sp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting Style "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting style \n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"font.family\"] = \"monospace\"\n",
    "\n",
    "colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown']\n",
    "\n",
    "# function to update the style of each axis\n",
    "def update_axis_style(ax, title = None, x_label='', y_label='', x_ticks=True,fsl = 12, fsb = 12, fst = 14):\n",
    "    if title:\n",
    "        ax.set_title(title,fontsize=fst)\n",
    "    ax.set_xlabel(x_label,fontsize=fsb)\n",
    "    ax.set_ylabel(y_label,fontsize=fsl)\n",
    "    ax.grid(True, which='both', linestyle='--', linewidth=0.5)\n",
    "    ax.tick_params(axis='x', which='both', bottom=x_ticks, top=False, labelbottom=x_ticks)\n",
    "    ax.tick_params(axis='y', which='both', left=True, right=False, labelleft=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sytem Definition using `msrDynamics`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up system \n",
    "\n",
    "MSRE = System()\n",
    "\n",
    "# define flow events\n",
    "df_spindown = pd.read_csv(os.getcwd() + '/data/spindown.csv', names = ['time', 'pct'], header = None)\n",
    "df_spinup = pd.read_csv(os.getcwd() + '/data/spinup.csv', names = ['time', 'pct'], header = None)\n",
    "spindown_fit = interp1d(df_spindown['time'], df_spindown['pct'])\n",
    "spinup_fit = interp1d(df_spinup['time'], df_spinup['pct'])\n",
    "t_spindown = 2500\n",
    "t_spinup = 7500\n",
    "\n",
    "# defining flow events\n",
    "def flow_fac(t):\n",
    "    if (t <= t_spindown):\n",
    "        return 1.0\n",
    "    elif (t < (t_spindown + 20.0)):\n",
    "        return max(spindown_fit(t-t_spindown)/100.0, 0.02)\n",
    "    elif (t <= t_spinup):\n",
    "        return 0.02\n",
    "    elif (t < t_spinup + 8.20):\n",
    "        return max(spinup_fit(t-t_spinup)/100.0, 0.02)\n",
    "    else:\n",
    "        return 1.0\n",
    "\n",
    "flow_pct = MSRE.add_input(flow_fac, T, max_anchors = 5000, input_tol = 32)\n",
    "\n",
    "# define nodes\n",
    "\n",
    "# radiator\n",
    "T_out_rc = Node(m=mn_rp, scp=mcp_rpn/mn_rp, W=W_rp, y0=T0_rp, name='T_out_rc')\n",
    "T_out_air = Node(m=mn_rs, scp=mcp_rsn/mn_rs, W=W_rs, y0=T0_rs, name='T_out_air')\n",
    "\n",
    "# heat exchanger\n",
    "T_hf1 = Node(m=mn_p, scp=mcp_pn/mn_p, W=W_p*flow_pct, y0=T0_p1, name='T_hf1')\n",
    "T_hf2 = Node(m=mn_p, scp=mcp_pn/mn_p, W=W_p*flow_pct, y0=T0_p2, name='T_hf2')\n",
    "T_hf3 = Node(m=mn_p, scp=mcp_pn/mn_p, W=W_p*flow_pct, y0=T0_p3, name='T_hf3')\n",
    "T_hf4 = Node(m=mn_p, scp=mcp_pn/mn_p, W=W_p*flow_pct, y0=T0_p4, name='T_hf4')\n",
    "T_ht1 = Node(m=m_tn, scp=scp_t, y0=T0_t1, name='T_ht1')\n",
    "T_ht2 = Node(m=m_tn, scp=scp_t, y0=T0_t2, name='T_ht2')\n",
    "T_hc1 = Node(m=mn_s, scp=mcp_sn/mn_s, W=W_s*flow_pct, y0=T0_s1, name='T_hc1')\n",
    "T_hc2 = Node(m=mn_s, scp=mcp_sn/mn_s, W=W_s*flow_pct, y0=T0_s2, name='T_hc2')\n",
    "T_hc3 = Node(m=mn_s, scp=mcp_sn/mn_s, W=W_s*flow_pct, y0=T0_s3, name='T_hc3')\n",
    "T_hc4 = Node(m=mn_s, scp=mcp_sn/mn_s, W=W_s*flow_pct, y0=T0_s4, name='T_hc4')\n",
    "\n",
    "# core \n",
    "n = Node(y0=n_frac0, name='n')\n",
    "C1 = Node(y0=C0[0], name='C1')\n",
    "C2 = Node(y0=C0[1], name='C2')\n",
    "C3 = Node(y0=C0[2], name='C3')\n",
    "C4 = Node(y0=C0[3], name='C4')\n",
    "C5 = Node(y0=C0[4], name='C5')\n",
    "C6 = Node(y0=C0[5], name='C6')\n",
    "rho = Node(y0=rho_0, name='rho')\n",
    "\n",
    "rho_control = Node(y0=0.0, name='rho_control')\n",
    "\n",
    "T_cg = Node(m=mcp_g1/scp_g, scp=scp_g, y0=T0_g1, name='T_cg')\n",
    "T_cf1 = Node(m=mn_f, scp=scp_f, W=W_f*flow_pct, y0=T0_f1, name='T_cf1')\n",
    "T_cf2 = Node(m=mn_f, scp=scp_f, W=W_f*flow_pct, y0=T0_f1, name='T_cf2')\n",
    "\n",
    "MSRE.add_nodes([T_out_rc,T_out_air,T_hf1,T_hf2,T_hf3,T_hf4,T_ht1,T_ht2,T_hc1,\n",
    "               T_hc2,T_hc3,T_hc4,n,C1,C2,C3,C4,C5,C6,T_cg,T_cf1,T_cf2,rho, rho_control])\n",
    "\n",
    "# dynamics\n",
    "\n",
    "# radiator\n",
    "T_out_rc.set_dTdt_advective(source = T_hc4.y(t-tau_hx_r/flow_pct))\n",
    "T_out_rc.set_dTdt_convective(source = [T_out_air.y()], hA = [hA_rpn*flow_pct])\n",
    "\n",
    "T_out_air.set_dTdt_advective(source = Trs_in)\n",
    "T_out_air.set_dTdt_convective(source = [T_out_rc.y()], hA = [hA_rsn])\n",
    "\n",
    "# heat exchanger\n",
    "T_hf1.set_dTdt_advective(source = T_cf2.y(t-tau_c_hx/flow_pct))\n",
    "T_hf1.set_dTdt_convective(source = [T_ht1.y()], hA = [hA_pn*flow_pct])\n",
    "\n",
    "T_hf2.set_dTdt_advective(source = T_hf1.y())\n",
    "T_hf2.dTdt_convective = T_hf1.dTdt_convective\n",
    "\n",
    "T_hf3.set_dTdt_advective(source = T_hf2.y())\n",
    "T_hf3.set_dTdt_convective(source = [T_ht2.y()], hA = [hA_pn*flow_pct])\n",
    "\n",
    "T_hf4.set_dTdt_advective(source = T_hf3.y())\n",
    "T_hf4.dTdt_convective = T_hf3.dTdt_convective\n",
    "\n",
    "T_ht1.set_dTdt_convective(source = [T_hf1.y(),T_hf1.y(),T_hc3.y(),T_hc3.y()], hA = [hA_pn*flow_pct,hA_pn*flow_pct,hA_sn*flow_pct,hA_sn*flow_pct])\n",
    "T_ht2.set_dTdt_convective(source = [T_hf3.y(),T_hf3.y(),T_hc1.y(),T_hc1.y()], hA = [hA_pn*flow_pct,hA_pn*flow_pct,hA_sn*flow_pct,hA_sn*flow_pct])\n",
    "\n",
    "T_hc1.set_dTdt_advective(source = T_out_rc.y(t-tau_r_hx/flow_pct))\n",
    "T_hc1.set_dTdt_convective(source = [T_ht2.y()], hA = [hA_sn*flow_pct])\n",
    "\n",
    "T_hc2.set_dTdt_advective(source = T_hc1.y())\n",
    "T_hc2.dTdt_convective = T_hc1.dTdt_convective\n",
    "\n",
    "T_hc3.set_dTdt_advective(source = T_hc2.y())\n",
    "T_hc3.set_dTdt_convective(source = [T_ht1.y()], hA = [hA_sn*flow_pct])\n",
    "\n",
    "T_hc4.set_dTdt_advective(source = T_hc3.y())\n",
    "T_hc4.dTdt_convective = T_hc3.dTdt_convective\n",
    "\n",
    "# core\n",
    "n.set_dndt(r = rho.y() + rho_control.y(), beta_eff = beta_t, Lambda = Lam, lam = lam, C = [C1.y(),C2.y(),C3.y(),C4.y(),C5.y(),C6.y()])\n",
    "C1.set_dcdt(n.y(), beta = beta[0],Lambda = Lam, lam = lam[0], t_c = tau_c/flow_pct, t_l = tau_l/flow_pct, flow = True)\n",
    "C2.set_dcdt(n.y(), beta = beta[1],Lambda = Lam, lam = lam[1], t_c = tau_c/flow_pct, t_l = tau_l/flow_pct, flow = True)\n",
    "C3.set_dcdt(n.y(), beta = beta[2],Lambda = Lam, lam = lam[2], t_c = tau_c/flow_pct, t_l = tau_l/flow_pct, flow = True)\n",
    "C4.set_dcdt(n.y(), beta = beta[3],Lambda = Lam, lam = lam[3], t_c = tau_c/flow_pct, t_l = tau_l/flow_pct, flow = True)\n",
    "C5.set_dcdt(n.y(), beta = beta[4],Lambda = Lam, lam = lam[4], t_c = tau_c/flow_pct, t_l = tau_l/flow_pct, flow = True)\n",
    "C6.set_dcdt(n.y(), beta = beta[5],Lambda = Lam, lam = lam[5], t_c = tau_c/flow_pct, t_l = tau_l/flow_pct, flow = True)\n",
    "\n",
    "# control rod response\n",
    "rho_control.dydt = -sp.Min(29.0e-5, sp.Max(-29.0e-5, n.dydt))\n",
    "\n",
    "\n",
    "T_cg.set_dTdt_convective(source = [T_cf1.y()], hA = [hA_fg*flow_pct])\n",
    "T_cg.set_dTdt_internal(source = [n.y()], k = [k_g*P])\n",
    "\n",
    "T_cf1.set_dTdt_advective(source = T_hf4.y(t-tau_hx_c/flow_pct))\n",
    "T_cf1.set_dTdt_convective(source = [T_cg.y()], hA = [k_1*hA_fg*flow_pct])\n",
    "T_cf1.set_dTdt_internal(source = [n.y()], k = [k_f1*P])\n",
    "\n",
    "T_cf2.set_dTdt_advective(source = T_cf1.y())\n",
    "T_cf2.dTdt_convective = T_cf1.dTdt_convective\n",
    "T_cf2.set_dTdt_internal(source = [n.y()], k = [k_f2*P])\n",
    "\n",
    "rho.set_drdt(sources = [T_cf1.dydt, T_cg.dydt], coeffs = [a_f,a_g])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "finalizing integrator...\n",
      "integrating...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Integration progress:   0%|          | 0/1000000 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating, compiling, and loading C code.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Integration progress: 100%|██████████| 1000000/1000000 [05:19<00:00, 3128.32it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "populating nodes objects solution vectors...\n"
     ]
    }
   ],
   "source": [
    "# solve \n",
    "sol_jit = MSRE.solve(T, max_delay = 30.0, abs_tol = 1e-12, rel_tol = 1e-10, min_step = 1e-13, populate_nodes = True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x75e806ae3d40>"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 83432 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_ornl_spindown = pd.read_csv(os.getcwd() + '/data/ornl_spindown.csv', header = None, names = ['time', 'pcm'])\n",
    "df_ornl_spinup = pd.read_csv(os.getcwd() + '/data/ornl_spinup.csv', header = None, names = ['time', 'pcm'])\n",
    "\n",
    "fig, axs = plt.subplots(1,2, figsize = (12,5))\n",
    "\n",
    "T_compare_spindown = T[(T >= 2500.0) & (T <= 2572.0)] -2500\n",
    "rho_compare_spindown = 200+(rho_control.y_out[(T >= 2500.0) & (T <= 2572.0)])*1e5\n",
    "\n",
    "T_compare_spinup = T[(T >= 7500.0) & (T <= 7547.0)] -7500\n",
    "rho_compare_spinup = 200+(rho_control.y_out[(T >= 7500.0) & (T <= 7547.0)])*1e5\n",
    "\n",
    "update_axis_style(axs[0], title = 'MSRE Pump Coastdown Test')\n",
    "axs[0].plot(T_compare_spindown, rho_compare_spindown, label = 'simulated, msrDynamics')\n",
    "axs[0].scatter(df_ornl_spindown['time'], df_ornl_spindown['pcm'], marker=\"2\",color=colors[2],label=\"measured, ORNL\")\n",
    "axs[0].set_xlim([0,72])\n",
    "axs[0].legend()\n",
    "\n",
    "update_axis_style(axs[1], title = 'MSRE Pump Startup Test')\n",
    "axs[1].plot(T_compare_spinup, rho_compare_spinup-rho_compare_spinup[0], label = 'simulated, msrDynamics')\n",
    "axs[1].scatter(df_ornl_spinup['time'], df_ornl_spinup['pcm'], marker=\"2\",color=colors[2],label=\"measured, ORNL\")\n",
    "axs[1].set_xlim([0,47])\n",
    "axs[1].legend()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "onion_dynamics_py312",
   "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.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
