{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from parameters_U233 import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from jitcdde import t\n",
    "from msrDynamics import Node, System\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# unpack ORNL data\n",
    "df_ORNL = pd.read_csv(f\"./data/ORNL_msre_{int(P)}MW_U233_insertion.csv\",names=['t','dP'])\n",
    "df_ORNL = df_ORNL.sort_values(df_ORNL.columns[0])\n",
    "df_simulink = pd.read_excel(f\"./data/simulink_msre_{P}MW_U233_insertion.xlsx\")\n",
    "i_trans = [i for i in range(len(df_simulink)) if df_simulink['time'][i] >= 2500]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = df_simulink['time']\n",
    "MSRE = System()\n",
    "\n",
    "# radiator\n",
    "T_out_rc = Node(m = mn_rp, scp = mcp_rpn/mn_rp, W = W_rp, y0 = T0_rp)\n",
    "T_out_air = Node(m = mn_rs, scp = mcp_rsn/mn_rs, W = W_rs, y0 = T0_rs)\n",
    "\n",
    "# heat exchanger\n",
    "T_hf1 = Node(m = mn_p, scp = mcp_pn/mn_p, W = W_p, y0 = T0_p1)\n",
    "T_hf2 = Node(m = mn_p, scp = mcp_pn/mn_p, W = W_p, y0 = T0_p2)\n",
    "T_hf3 = Node(m = mn_p, scp = mcp_pn/mn_p, W = W_p, y0 = T0_p3)\n",
    "T_hf4 = Node(m = mn_p, scp = mcp_pn/mn_p, W = W_p, y0 = T0_p4)\n",
    "T_ht1 = Node(m = m_tn, scp = scp_t, y0 = T0_t1)\n",
    "T_ht2 = Node(m = m_tn, scp = scp_t, y0 = T0_t2)\n",
    "T_hc1 = Node(m = mn_s, scp = mcp_sn/mn_s, W = W_s, y0 = T0_s1)\n",
    "T_hc2 = Node(m = mn_s, scp = mcp_sn/mn_s, W = W_s, y0 = T0_s2)\n",
    "T_hc3 = Node(m = mn_s, scp = mcp_sn/mn_s, W = W_s, y0 = T0_s3)\n",
    "T_hc4 = Node(m = mn_s, scp = mcp_sn/mn_s, W = W_s, y0 = T0_s4)\n",
    "\n",
    "# core \n",
    "n = Node(y0 = n_frac0)\n",
    "C1 = Node(y0 = C0[0])\n",
    "C2 = Node(y0 = C0[1])\n",
    "C3 = Node(y0 = C0[2])\n",
    "C4 = Node(y0 = C0[3])\n",
    "C5 = Node(y0 = C0[4])\n",
    "C6 = Node(y0 = C0[5])\n",
    "rho = Node(y0 = 0.0)\n",
    "\n",
    "# add reactivity input\n",
    "t_ins = 2500\n",
    "def rho_insert(t):\n",
    "    if (t<t_ins):\n",
    "        return 0.0\n",
    "    else:\n",
    "        return inserted\n",
    "\n",
    "rho_ext = MSRE.add_input(rho_insert, T)\n",
    "\n",
    "T_cg = Node(m = mcp_g1/scp_g, scp = scp_g, y0 = T0_g1)\n",
    "T_cf1 = Node(m = mn_f, scp = scp_f, W = W_f, y0 = T0_f1)\n",
    "T_cf2 = Node(m = mn_f, scp = scp_f, W = W_f, y0 = T0_f2)\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])\n",
    "\n",
    "# dynamics \n",
    "\n",
    "# radiator\n",
    "T_out_rc.set_dTdt_advective(source = T_hc4.y(t-tau_hx_r))\n",
    "T_out_rc.set_dTdt_convective(source = [T_out_air.y()], hA = [hA_rpn])\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))\n",
    "T_hf1.set_dTdt_convective(source = [T_ht1.y()], hA = [hA_pn])\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])\n",
    "\n",
    "T_hf4.set_dTdt_advective(source = T_hf3.y())\n",
    "# T_hf4.set_dTdt_convective(source = [T_ht2.y()], hA = [hA_pn])\n",
    "T_hf4.dTdt_convective = T_hf3.dTdt_convective\n",
    "\n",
    "# T_ht1.set_dTdt_convective(source = [T_hf1.y(),T_hf2.y(),T_hc3.y(),T_hc4.y()], hA = [hA_pn,hA_pn,hA_sn,hA_sn])\n",
    "# T_ht2.set_dTdt_convective(source = [T_hf3.y(),T_hf4.y(),T_hc1.y(),T_hc2.y()], hA = [hA_pn,hA_pn,hA_sn,hA_sn])\n",
    "T_ht1.set_dTdt_convective(source = [T_hf1.y(),T_hf1.y(),T_hc3.y(),T_hc3.y()], hA = [hA_pn,hA_pn,hA_sn,hA_sn])\n",
    "T_ht2.set_dTdt_convective(source = [T_hf3.y(),T_hf3.y(),T_hc1.y(),T_hc1.y()], hA = [hA_pn,hA_pn,hA_sn,hA_sn])\n",
    "\n",
    "T_hc1.set_dTdt_advective(source = T_out_rc.y(t-tau_r_hx))\n",
    "T_hc1.set_dTdt_convective(source = [T_ht2.y()], hA = [hA_sn])\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])\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_ext, 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[0],Lam,lam[0],True,tau_c,tau_l)\n",
    "C2.set_dcdt(n.y(),beta[1],Lam,lam[1],True,tau_c,tau_l)\n",
    "C3.set_dcdt(n.y(),beta[2],Lam,lam[2],True,tau_c,tau_l)\n",
    "C4.set_dcdt(n.y(),beta[3],Lam,lam[3],True,tau_c,tau_l)\n",
    "C5.set_dcdt(n.y(),beta[4],Lam,lam[4],True,tau_c,tau_l)\n",
    "C6.set_dcdt(n.y(),beta[5],Lam,lam[5],True,tau_c,tau_l)\n",
    "\n",
    "T_cg.set_dTdt_convective(source = [T_cf1.y()], hA = [hA_fg])\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))\n",
    "T_cf1.set_dTdt_convective(source = [T_cg.y()], hA = [k_1*hA_fg])\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_cf2.dydt, T_cg.dydt], coeffs = [a_f/2,a_f/2,a_g])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating, compiling, and loading C code.\n",
      "Using default integration parameters.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/luke/git/envs/onion_env/lib/python3.11/site-packages/jitcdde/_jitcdde.py:795: UserWarning: You did not explicitly handle initial discontinuities. Proceed only if you know what you are doing. This is only fine if you somehow chose your initial past such that the derivative of the last anchor complies with the DDE. In this case, you can set the attribute `initial_discontinuities_handled` to `True` to suppress this warning. See https://jitcdde.rtfd.io/#discontinuities for details.\n",
      "  warn(\"You did not explicitly handle initial discontinuities. Proceed only if you know what you are doing. This is only fine if you somehow chose your initial past such that the derivative of the last anchor complies with the DDE. In this case, you can set the attribute `initial_discontinuities_handled` to `True` to suppress this warning. See https://jitcdde.rtfd.io/#discontinuities for details.\")\n",
      "/home/luke/git/envs/onion_env/lib/python3.11/site-packages/jitcdde/_jitcdde.py:795: UserWarning: You did not explicitly handle initial discontinuities. Proceed only if you know what you are doing. This is only fine if you somehow chose your initial past such that the derivative of the last anchor complies with the DDE. In this case, you can set the attribute `initial_discontinuities_handled` to `True` to suppress this warning. See https://jitcdde.rtfd.io/#discontinuities for details.\n",
      "  warn(\"You did not explicitly handle initial discontinuities. Proceed only if you know what you are doing. This is only fine if you somehow chose your initial past such that the derivative of the last anchor complies with the DDE. In this case, you can set the attribute `initial_discontinuities_handled` to `True` to suppress this warning. See https://jitcdde.rtfd.io/#discontinuities for details.\")\n",
      "/home/luke/git/envs/onion_env/lib/python3.11/site-packages/jitcdde/_jitcdde.py:792: UserWarning: The target time is smaller than the current time. No integration step will happen. The returned state will be extrapolated from the interpolating Hermite polynomial for the last integration step. You may see this because you try to integrate backwards in time, in which case you did something wrong. You may see this just because your sampling step is small, in which case there is no need to worry.\n",
      "  warn(\"The target time is smaller than the current time. No integration step will happen. The returned state will be extrapolated from the interpolating Hermite polynomial for the last integration step. You may see this because you try to integrate backwards in time, in which case you did something wrong. You may see this just because your sampling step is small, in which case there is no need to worry.\")\n"
     ]
    }
   ],
   "source": [
    "sol_jit = MSRE.solve(T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_199186/112360540.py:1: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
      "  duration = df_ORNL.iloc[-1][0]-df_ORNL.iloc[0][0]\n",
      "/tmp/ipykernel_199186/112360540.py:8: FutureWarning: ChainedAssignmentError: behaviour will change in pandas 3.0!\n",
      "You are setting values through chained assignment. Currently this works in certain cases, but when using Copy-on-Write (which will become the default behaviour in pandas 3.0) this will never work to update the original DataFrame or Series, because the intermediate object on which we are setting values will behave as a copy.\n",
      "A typical example is when you are setting values in a column of a DataFrame, like:\n",
      "\n",
      "df[\"col\"][row_indexer] = value\n",
      "\n",
      "Use `df.loc[row_indexer, \"col\"] = values` instead, to perform the assignment in a single step and ensure this keeps updating the original `df`.\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "\n",
      "  df_ORNL['t'][i] += delta\n"
     ]
    }
   ],
   "source": [
    "duration = df_ORNL.iloc[-1][0]-df_ORNL.iloc[0][0]\n",
    "i_insert = [i for i in range(len(T)) if (T[i] > t_ins) and (T[i] < t_ins + duration)]\n",
    "ref_P = P*n.y_out[i_insert[0]-100]\n",
    "dP = [(k*P)-ref_P for k in n.y_out]\n",
    "\n",
    "delta = t_ins - df_ORNL['t'][0] - t_ins\n",
    "for i in range(len(df_ORNL)):\n",
    "    df_ORNL['t'][i] += delta\n",
    "\n",
    "ref_P_simulink = df_simulink['Mux(4)'][i_trans[0]-100]*P\n",
    "i_window = [i for i in i_trans if df_simulink['time'][i]-t_ins <= duration]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 94556 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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, x_label='', y_label='', x_ticks=True):\n",
    "    ax.set_title(title,fontsize=20)\n",
    "    ax.set_xlabel(x_label,fontsize=14)\n",
    "    ax.set_ylabel(y_label,fontsize=14)\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)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "update_axis_style(ax,f\"Step Reactivity Insertion Response: {int(inserted*(10**5))}pcm, {int(P)}MW, U233 Fuel\")\n",
    "ax.plot(T[i_insert[0]:i_insert[-1]+1]-t_ins,dP[i_insert[0]:i_insert[-1]+1],label=\"msrDynamics\",color=colors[0])\n",
    "ax.plot(df_simulink['time'][i_window]-t_ins,df_simulink['Mux(4)'][i_window]*P-ref_P_simulink,color=colors[1],linestyle=\"--\",label=\"Simulink, Singh et al.\")\n",
    "ax.plot(df_ORNL.iloc[:,0],df_ORNL.iloc[:,1],label=\"ORNL-TM-2997\",color=colors[2],linestyle=\"--\")\n",
    "ax.set_ylabel(r\"$\\Delta P$\")\n",
    "ax.set_xlabel(r\"$t$(s)\")\n",
    "ax.legend()\n",
    "\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "thesis_env",
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
