Input_NetCDF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8efdbeef",
   "metadata": {},
   "source": [
    "# Create Input File:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6ab36f96",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "from ipywidgets import HBox, Label\n",
    "from netCDF4 import Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "df6562bf",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# stop execution class in case of invalid input parameters\n",
    "class StopExecution(Exception):\n",
    "    def _render_tracebeck_(self):\n",
    "        pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e69edd9",
   "metadata": {},
   "source": [
    "### Create Input File"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97eda80f",
   "metadata": {},
   "source": [
    "The following cells allow you to set the input parameters required for the half-SPM simulation. You can set the input parameters within the suggested range and the current default values are taken from _Chen 2020 (https:/doi.org/10.1149/1945-7111/ab9050)_ for a NCM (Nickel-Cobalt-Manganese) positive electrode. Once run, a NetCDF file called \"__SPM_input.nc__\" will be created which can be read by the program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "aa6e7362",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ea5a1771792a41bcb40650f69cfbc025",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Temperature [°C]'), FloatSlider(value=21.0, max=50.0, min=-20.0, step=1.0)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2e416e8af7c2482787d4d4788e391f60",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Mean particle radius [$\\\\mu m$]'), FloatSlider(value=5.22, max=7.22, min=3.22, ste…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7affb87cc182486b8d43c70ceb6cf7ab",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Electrode thickness [$\\\\mu m$]'), FloatSlider(value=75.6, min=50.0)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f7cba4acda0347178f73ef9ad3ba74c5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Reaction rate coefficent[$Am^{-2}(m^3mol^{-1})^{1.5}$]'), FloatSlider(value=3.42, …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7e550d06561e42ecbfceca9704d86ba0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Diffusion coefficient [$10^{-15} m^2 s^{-1}$]'), FloatSlider(value=1.48, step=0.01…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6b3df624ba7443d5820d49d7bce0ec93",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Maximum lithium concentration [$mol m^{-3}$]'), FloatSlider(value=51765.0, max=100…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3ba34be64aff4be1aef9397cb4ae8bbc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='State of charge [%]'), FloatSlider(value=0.0, step=1.0)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "91a1c89eda024ce4b0089de6214110da",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Applied current [A]'), FloatSlider(value=5.0, max=20.0, min=-20.0)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c85ea588d387415689c824188162a0b4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Active material volume fraction [%]'), FloatSlider(value=66.5)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "76df7aaad1a4437fa9d632a59b9032f2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Electrode plate area [$m^2$]'), FloatSlider(value=0.1027, max=1.0, min=0.001, step…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1f36dc6eb0b04eafb6a2ae315b936fbf",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='No. of simulation steps'), BoundedIntText(value=1000, max=1000000000, min=1)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bfdc7dda7ed241b6b841f71e6d633276",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Time step (s)'), FloatSlider(value=2.0, min=0.01, step=0.01)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "68adc38ce4aa4c41b789aef63694093b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Output written every [n] steps'), BoundedIntText(value=5, max=1000000000, min=1)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dc183d42a1fd4e54bc206697a89590c8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='No. of space steps'), IntSlider(value=20, max=1000, min=2)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "da543c3ee58941fc87f15bfc53d09c1b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Output voltage data'), Select(options=('Yes', 'No'), value='Yes')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# set input parameter within given ranges\n",
    "Temp_ = HBox([Label('Temperature [°C]'), widgets.FloatSlider(min=-20.0, max=50.0, value=21.0, step=1)])\n",
    "Rad_ = HBox([Label('Mean particle radius [$\\mu m$]'), widgets.FloatSlider(min=3.22, max=7.22, value=5.22, step=0.01)])\n",
    "Thick_ = HBox([Label('Electrode thickness [$\\mu m$]'), widgets.FloatSlider(min=50.0, max=100.0, value=75.6, step=0.1)])\n",
    "Rr_coef_ = HBox([Label('Reaction rate coefficent[$Am^{-2}(m^3mol^{-1})^{1.5}$]'), widgets.FloatSlider(min=1.5, max=6.5, value=3.42, step=0.01)])\n",
    "Dif_coef_ = HBox([Label('Diffusion coefficient [$10^{-15} m^2 s^{-1}$]'), widgets.FloatSlider(min=0, max=100, value=1.48, step=0.01)])\n",
    "Max_c_ = HBox([Label('Maximum lithium concentration [$mol m^{-3}$]'), widgets.FloatSlider(min=0.0, max=100000.0, value=51765.0, step=1)])\n",
    "SOC_ = HBox([Label('State of charge [%]'), widgets.FloatSlider(min=0, max=100, value=0, step=1)])\n",
    "Current_ = HBox([Label('Applied current [A]'), widgets.FloatSlider(min=-20.0, max=20.0, value=5.0, step=0.1)])\n",
    "Vol_per_ = HBox([Label('Active material volume fraction [%]'), widgets.FloatSlider(min = 0.0, max=100.0, value=66.5, step=0.1)])\n",
    "Area_ = HBox([Label('Electrode plate area [$m^2$]'), widgets.FloatSlider(min=0.001, max=1.0, value=0.1027, step=0.001)])\n",
    "Sim_steps_ = HBox([Label('No. of simulation steps'), widgets.BoundedIntText(value=1000, min=1, max=10**9, disabled=False)])\n",
    "Dt_ = HBox([Label('Time step (s)'), widgets.FloatSlider(min=0.01, max=100.0, value=2.0, step=0.01)]) \n",
    "Out_steps_ = HBox([Label('Output written every [n] steps'), widgets.BoundedIntText(value=5, min=1, max=10**9, disabled=False)])\n",
    "Space_steps_ = HBox([Label('No. of space steps'), widgets.IntSlider(min=2, max=1000, value=20, step=1)])\n",
    "Volt_do_ =  HBox([Label('Output voltage data'), widgets.Select(options=['Yes', 'No'], value='Yes', disabled=False)])\n",
    "\n",
    "display(Temp_, Rad_, Thick_, Rr_coef_, Dif_coef_, Max_c_, SOC_, Current_, Vol_per_, Area_, Sim_steps_, Dt_, Out_steps_, Space_steps_, Volt_do_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "821d4560",
   "metadata": {},
   "source": [
    "__Set parameters to whatever you want:__\n",
    "\n",
    "If you feel confident in using the program, you can manually override the sliders using the boxes below. This is recommended for experienced users and only physical valid values should be set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "29766fd7",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1ef61fbdf1a34b5ea7e279e887d71bf3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Temperature [°C]'), FloatText(value=21.0)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ba072507adb7411e8aceb14dfef7de22",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Mean particle radius [$\\\\mu m$]'), FloatText(value=5.22)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1d2310e692db433fb98bc79813f04743",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Electrode thickness [$\\\\mu m$]'), FloatText(value=75.6)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "63193b7225fb41ed9a44df2855015861",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Reaction rate coeffic[ent($Am^{-2}(m^3mol^{-1})^{1.5}$]'), FloatText(value=3.42)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2b5b6b48676749bb955251a4023e9ed8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Diffusion coefficient [$10^{-15} m^2 s^{-1}$]'), FloatText(value=1.48)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "81c6858216244302a981d2a3d3af0b18",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Maximum lithium concentration [$mol m^{-3}$]'), FloatText(value=51765.0)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a071da39ea3f4048ab4050cc943d4dff",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Applied current [A]'), FloatText(value=5.0)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1b28b50c46d746ed8315bf077e393220",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Electrode area [$m^2$]'), FloatText(value=0.1027)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "94a334514bd64798b664331180e0a1ac",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Time step (s)'), FloatText(value=2.0)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f87b688e5fa8401ba7be9110208b6472",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='No. of space steps'), IntText(value=20)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# unbounded boxes to set parameters outside of suggested ranges if wanted\n",
    "Temp_ = HBox([Label('Temperature [°C]'), widgets.FloatText(value=Temp_.children[1].value, disabled=False)])\n",
    "Rad_ = HBox([Label('Mean particle radius [$\\mu m$]'), widgets.FloatText(value=Rad_.children[1].value, disabled=False)])\n",
    "Thick_ = HBox([Label('Electrode thickness [$\\mu m$]'), widgets.FloatText(value=Thick_.children[1].value, disabled=False)])\n",
    "Rr_coef_ = HBox([Label('Reaction rate coeffic[ent($Am^{-2}(m^3mol^{-1})^{1.5}$]'), widgets.FloatText(value=Rr_coef_.children[1].value, disabled=False)])\n",
    "Dif_coef_ = HBox([Label('Diffusion coefficient [$10^{-15} m^2 s^{-1}$]'), widgets.FloatText(value=Dif_coef_.children[1].value, disabled=False)])\n",
    "Max_c_ = HBox([Label('Maximum lithium concentration [$mol m^{-3}$]'), widgets.FloatText(value=Max_c_.children[1].value, disabled=False)])\n",
    "Current_ = HBox([Label('Applied current [A]'), widgets.FloatText(value=Current_.children[1].value, disabled=False)])\n",
    "Area_ = HBox([Label('Electrode area [$m^2$]'), widgets.FloatText(value=Area_.children[1].value, disabled=False)])\n",
    "Dt_ = HBox([Label('Time step (s)'), widgets.FloatText(value=Dt_.children[1].value, disabled=False)])\n",
    "Space_steps_ = HBox([Label('No. of space steps'), widgets.IntText(value=Space_steps_.children[1].value, disabled=False)])\n",
    "\n",
    "display(Temp_, Rad_, Thick_, Rr_coef_, Dif_coef_, Max_c_, Current_, Area_, Dt_, Space_steps_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd8040c8",
   "metadata": {},
   "source": [
    "The remaining cells extract the set input parameters, convert them to SI units, and create a NetCDF input file that can be read by the programme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a5fd3743",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# saving the values set by the user and adjusting relevant parameters to be in SI units\n",
    "error = False\n",
    "Temp = Temp_.children[1].value + 273.15   #K\n",
    "if Temp <= 0:\n",
    "    error = True\n",
    "    print('Error, Temperature cannot be lower than 0K or -273.15 C, please reset!')\n",
    "Rad = Rad_.children[1].value * 10**(-6)   #m\n",
    "if Rad <= 0:\n",
    "    error = True\n",
    "    print('Error, Mean particle radius cannot be 0 or negative, please reset!')\n",
    "Thick = Thick_.children[1].value * 10**(-6)   #m\n",
    "if Thick <= 0:\n",
    "    error = True\n",
    "    print('Error, Electrode thickness cannot be 0 or negative, please reset!')\n",
    "Rr_coef =  Rr_coef_.children[1].value\n",
    "if Rr_coef <= 0:\n",
    "    error = True\n",
    "    print('Error, Reaction rate coefficient cannot be 0 or negative, please reset!')\n",
    "Dif_coef = Dif_coef_.children[1].value * 10**(-15)\n",
    "if Dif_coef <= 0:\n",
    "    error = True\n",
    "    print('Error, Doffusion coefficient cannot be 0 or negative, please reset!')\n",
    "Max_c = Max_c_.children[1].value\n",
    "if Max_c <= 0:\n",
    "    error = True\n",
    "    print('Error, Maximum lithium concentration cannot be 0 or negative, please reset!')\n",
    "# Init_c derived from values in Chen2020:\n",
    "# stoichiometry at 0% SOC = 0.2661, and at 100% SOC = 0.9084\n",
    "# assuming linear behaviour we get x = mcx +c with c = 0.9084\n",
    "# and m = 0.2661-0.9084 = -0.6423\n",
    "Init_c = Max_c * ((SOC_.children[1].value/100)*(-0.6423) + 0.9084)\n",
    "Iapp = Current_.children[1].value / Area_.children[1].value  # A/m^2\n",
    "Vol_per = Vol_per_.children[1].value\n",
    "Sim_steps = Sim_steps_.children[1].value\n",
    "if Sim_steps <= 0:\n",
    "    error = True\n",
    "    print('Error, the number of simulation steps cannot be 0 or negative, please reset!')\n",
    "Dt = Dt_.children[1].value\n",
    "if Dt <= 0:\n",
    "    error = True\n",
    "    print('Error, the time step cannot be 0 or negative, please reset!')\n",
    "Out_steps = Out_steps_.children[1].value\n",
    "if Out_steps <= 0:\n",
    "    error = True\n",
    "    print('Error, the number of output steps cannot be 0 or negative, please reset!')\n",
    "Space_steps = Space_steps_.children[1].value\n",
    "if Space_steps <= 0:\n",
    "    error = True\n",
    "    print('Error, the number of space steps cannot be 0 negative, please reset!')\n",
    "if Volt_do_.children[1].value == 'Yes':\n",
    "    Volt_do = 1\n",
    "else:\n",
    "    Volt_do = 0\n",
    "if error == True:\n",
    "    print('After the parameter has been changed, please rerun this cell.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "591b283a",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "#create new NetCDF file in 'writing mode' and 'NETCDF4 format'\n",
    "from netCDF4 import Dataset\n",
    "# check all input parameters are valid\n",
    "if error == True:\n",
    "    'Please adjust the input parameter indicated by the error message from the previous cell before creating the input file.'\n",
    "    raise StopExecution\n",
    "rootgrp = Dataset('SPM_input.nc', 'w', format='NETCDF4')\n",
    "#creating dimensions for vector holding all the input variables\n",
    "temp_dim = rootgrp.createDimension('temp_dim', 1)\n",
    "rad_dim = rootgrp.createDimension('rad_dim', 1)\n",
    "thick_dim = rootgrp.createDimension('thick_dim', 1)\n",
    "rr_coef_dim = rootgrp.createDimension('rr_coef_dim', 1)\n",
    "dif_coef_dim = rootgrp.createDimension('dif_coef_dim', 1)\n",
    "max_c_dim = rootgrp.createDimension('max_c_dim', 1)\n",
    "init_c_dim = rootgrp.createDimension('init_c_dim', 1)\n",
    "iapp_dim = rootgrp.createDimension('iapp_dim', 1)\n",
    "vol_per_dim = rootgrp.createDimension('vol_per_dim', 1)\n",
    "sim_steps_dim = rootgrp.createDimension('sim_steps_dim', 1)\n",
    "dt_dim = rootgrp.createDimension('dt_dim', 1)\n",
    "out_steps_dim = rootgrp.createDimension('out_steps_dim', 1)\n",
    "space_steps_dim = rootgrp.createDimension('space_steps_dim', 1)\n",
    "volt_do_dim = rootgrp.createDimension('volt_do_dim', 1)\n",
    "checkpoint_dim = rootgrp.createDimension('checkpoint_dim', 1)\n",
    "#creating variable \n",
    "temp = rootgrp.createVariable('temp', 'f8', ('temp_dim',))\n",
    "rad = rootgrp.createVariable('rad', 'f8', ('rad_dim',))\n",
    "thick = rootgrp.createVariable('thick', 'f8', ('thick_dim',))\n",
    "rr_coef = rootgrp.createVariable('rr_coef', 'f8', ('rr_coef_dim',))\n",
    "dif_coef = rootgrp.createVariable('dif_coef', 'f8', ('dif_coef_dim',))\n",
    "max_c = rootgrp.createVariable('max_c', 'f8', ('max_c_dim',))\n",
    "init_c = rootgrp.createVariable('init_c', 'f8', ('init_c_dim',))\n",
    "iapp = rootgrp.createVariable('iapp', 'f8', ('iapp_dim',))\n",
    "vol_per = rootgrp.createVariable('vol_per', 'f8', ('vol_per_dim',))\n",
    "sim_steps = rootgrp.createVariable('sim_steps', 'i4', ('sim_steps_dim',))\n",
    "dt = rootgrp.createVariable('dt', 'f8', ('dt_dim',))\n",
    "out_steps = rootgrp.createVariable('out_steps', 'i4', ('out_steps_dim',))\n",
    "space_steps = rootgrp.createVariable('space_steps', 'i4', ('space_steps_dim',))\n",
    "volt_do = rootgrp.createVariable('volt_do', 'i4', ('volt_do_dim',))\n",
    "checkpoint = rootgrp.createVariable('checkpoint', 'i4', ('checkpoint_dim',))\n",
    "# attributes\n",
    "rootgrp.description = 'Input parameters for SMP model'\n",
    "temp.description = 'Temperature'\n",
    "temp.units = 'K'\n",
    "rad.description = 'Mean particle radius'\n",
    "rad.units = 'm'\n",
    "thick.description = 'Electrode thickness'\n",
    "thick.units = 'm'\n",
    "rr_coef.description = 'Reaction rate coefficient'\n",
    "rr_coef.units = '$Am^{-2}(m^3mol^{-1})^{1.5}$'\n",
    "dif_coef.description ='Diffusion coefficient'\n",
    "dif_coef.units = '$m^2 s^{-1}$'\n",
    "max_c.description = 'Maximum lithium concentration'\n",
    "max_c.units = '$mol m^{-3}$'\n",
    "init_c.description = 'Initial lithium concentration'\n",
    "init_c.units = '$mol m^{-3}$'\n",
    "iapp.description = 'Applied current density'\n",
    "iapp.units = '$A/m^2$'\n",
    "vol_per.description = 'Active material volume fraction'\n",
    "vol_per.units = '%'\n",
    "sim_steps.description = 'Total number of simultion steps'\n",
    "sim_steps.units = 'unitless'\n",
    "dt.description = 'Time step'\n",
    "dt.units = 's'\n",
    "out_steps.description = 'Output written every [n] number of steps'\n",
    "out_steps.units = 'untiless'\n",
    "space_steps.description = 'No. of space steps'\n",
    "space_steps.units = 'unitless'\n",
    "volt_do.description = 'Write voltage data'\n",
    "volt_do.units = 'unitless'\n",
    "checkpoint.description = 'Starting from checkpont file'\n",
    "checkpoint.units = 'unitless'\n",
    "# writing data to input_parameters variable\n",
    "temp[0] = Temp\n",
    "rad[0] = Rad\n",
    "thick[0] = Thick\n",
    "rr_coef[0] = Rr_coef\n",
    "dif_coef[0] = Dif_coef\n",
    "max_c[0] = Max_c\n",
    "init_c[0] = Init_c\n",
    "iapp[0] = Iapp\n",
    "vol_per[0] = Vol_per\n",
    "sim_steps[0] = Sim_steps\n",
    "dt[0] = Dt\n",
    "out_steps[0] = Out_steps\n",
    "space_steps[0] = Space_steps\n",
    "volt_do[0] = Volt_do\n",
    "checkpoint[0] = 0\n",
    "#closing the NetCDF file\n",
    "rootgrp.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "591981a5",
   "metadata": {},
   "source": [
    "### Start From Checkpoint"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e0e70b8",
   "metadata": {},
   "source": [
    "To run from a checkpoint you need to set \"__Starting from checkpoint file__\" to \"__Yes__\". You should then enter in the number of extra simulation steps you want to perform along with how often you want to write to the output file. This cell produces a modified input file that you can now put back into the program.\n",
    "\n",
    "It is important that the input and output file from your earlier run still exists, so that the new data has a place to be written.\n",
    "\n",
    "A new simulation is started by setting \"__Starting from checkpoint file__\" to \"__No__\", which is the default."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7b62a15d",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "74f9e1934c9343dda69870650c9deb81",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Starting from checkpoint file'), Select(index=1, options=('Yes', 'No'), value='No'…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6427dfb5f29f43eba0d9906f877ea3ee",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='No. of simulation steps'), BoundedIntText(value=1000, max=1000000000, min=1)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d7f738cd887f44618d974db44f5142d1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Output written every [n] steps'), BoundedIntText(value=5, max=1000000000, min=1)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c83863c86fa744ec9e787fe9f0042c9c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Applied current [A]'), BoundedFloatText(value=5.0, max=20.0, min=-20.0)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# parameters for running from checkpoint file\n",
    "Checkpoint_ =  HBox([Label('Starting from checkpoint file'), widgets.Select(options=['Yes', 'No'], value='No', disabled=False)])\n",
    "Sim_steps_c = HBox([Label('No. of simulation steps'), widgets.BoundedIntText(value=1000, min=1, max=10**9, disabled=False)])\n",
    "Out_steps_c = HBox([Label('Output written every [n] steps'), widgets.BoundedIntText(value=5, min=1, max=10**9, disabled=False)])\n",
    "Iapp_c = HBox([Label('Applied current [A]'), widgets.BoundedFloatText(value=5, min=-20, max=20, disabled=False)])\n",
    "display(Checkpoint_, Sim_steps_c, Out_steps_c, Iapp_c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1e839b98",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# write chaged parameters to the input file\n",
    "if Checkpoint_.children[1].value == 'Yes':\n",
    "        Checkpoint = 1\n",
    "        rootgrp = Dataset('SPM_input.nc', 'r+', format='NETCDF4')\n",
    "        rootgrp['checkpoint'][:] = Checkpoint\n",
    "        rootgrp['sim_steps'][:] = Sim_steps_c.children[1].value \n",
    "        rootgrp['out_steps'][:] = Out_steps_c.children[1].value \n",
    "        rootgrp['iapp'][:] = Iapp_c.children[1].value\n",
    "        rootgrp.close()\n",
    "else:\n",
    "        Checkpoint = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2bbb371",
   "metadata": {},
   "source": [
    "### Uncertainty Propagation and Sensitivity Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41a6dccc",
   "metadata": {},
   "source": [
    "To perform uncertainty propagation you need to select \"__Yes__\" for \"__conduct uncertainty propagation__\" and then enter in the number of samples you want to take to approximate the uncertainty. You then need to enter in the standard deviations associated with your input parameters.\n",
    "\n",
    "If you don't want to consider the uncertainty created by a variable, simply set its standard deviation to 0 for it to be ignored.\n",
    "\n",
    "Just with the checkpoint system, your input file \"__SMP_input.nc__\" must already be created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e36606d2",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "052bb512460046fdbbfa852e41a6091d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Conduct uncertainty propagation'), Select(index=1, options=('Yes', 'No'), value='N…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cef214c9376144df890da926a6b8fb16",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='No. of samples [unitless]'), IntText(value=10)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3a2ae90d4f984e1da646cd929ac0da38",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Standard devidation (std) of temperature [°C]'), FloatText(value=14.9075)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ed8b7b73a6e0458a90d43a5a960e9cb3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Std of mean particle radius [$\\\\mu m$]'), FloatText(value=0.261)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a3272d0944ea4e94a331e7842fff8758",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Std of electrode thickness [$\\\\mu m$]'), FloatText(value=3.78)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4f8ac58b65bd47738ffe1e8be3a6e5cd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Std of reaction reate coefficient [$Am^{-2}(m^3mol^{-1})^{1.5}$]'), FloatText(valu…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d1fa71499a8544b7a6ef9979bd3fa594",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Std of diffusion coefficient [$10^{-15} m^2 s^{-1}$]'), FloatText(value=0.074)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3c44016d7cf34308a678631c3923f39e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Std of initial concentration [$mol m^{-3}$]'), FloatText(value=0.0047023326)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "beb2e2e4708f412cbb82185844d70d06",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Std of maximum concentration [$mol m^{-3}$]'), FloatText(value=0.0051765)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "30941756c4af4bc9bebf5d36011bba5b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Std of active material volume fraction [%]'), FloatText(value=3.325)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ebd3ed471cb34a9780ac2215595d19cc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Label(value='Std of applied current [A]'), FloatText(value=0.25)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# option to conduct uncertainty propagation\n",
    "UQ_ = HBox([Label('Conduct uncertainty propagation'), widgets.Select(options=['Yes', 'No'], value='No', disabled=False)])\n",
    "No_samples_ = HBox([Label('No. of samples [unitless]'), widgets.IntText(value=10, disabled=False)])\n",
    "Temp_std_ = HBox([Label('Standard devidation (std) of temperature [°C]'), widgets.FloatText(value=14.9075, disabled=False)])\n",
    "Rad_std_ = HBox([Label('Std of mean particle radius [$\\mu m$]'), widgets.FloatText(value=5.22*0.05, disabled=False)])\n",
    "Thick_std_ = HBox([Label('Std of electrode thickness [$\\mu m$]'), widgets.FloatText(value=75.6*0.05, disabled=False)])\n",
    "Rr_coef_std_ = HBox([Label('Std of reaction reate coefficient [$Am^{-2}(m^3mol^{-1})^{1.5}$]'), widgets.FloatText(value=6.5*0.05, disabled=False)])\n",
    "Dif_coef_std_ = HBox([Label('Std of diffusion coefficient [$10^{-15} m^2 s^{-1}$]'), widgets.FloatText(value=1.48*0.05, disabled=False)])\n",
    "Init_c_std_ = HBox([Label('Std of initial concentration [$mol m^{-3}$]'), widgets.FloatText(value=47023.326*(10**-7), disabled=False)])\n",
    "Max_c_std_ = HBox([Label('Std of maximum concentration [$mol m^{-3}$]'), widgets.FloatText(value=51765*(10**-7), disabled=False)])\n",
    "Vol_per_std_ = HBox([Label('Std of active material volume fraction [%]'), widgets.FloatText(value=66.5*0.05, disabled=False)])\n",
    "Iapp_std_ = HBox([Label('Std of applied current [A]'), widgets.FloatText(value=5.0*0.05, disabled=False)])\n",
    "display(UQ_, No_samples_, Temp_std_, Rad_std_, Thick_std_, Rr_coef_std_, Dif_coef_std_, Init_c_std_, Max_c_std_, Vol_per_std_, Iapp_std_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ff5a55d1",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# saving the set uQ values\n",
    "No_samples = No_samples_.children[1].value\n",
    "if No_samples <= 0:\n",
    "    error = True\n",
    "    print('Error, Number of samples cannot be 0 or negative, please reset!')\n",
    "Temp_std = Temp_std_.children[1].value #K\n",
    "if Temp_std <= 0:\n",
    "    error = True\n",
    "    print('Error, Std of temperature cannot be 0 or negative, please reset!')\n",
    "Rad_std = Rad_std_.children[1].value * 10**(-6)   #m\n",
    "if Rad_std <= 0:\n",
    "    error = True\n",
    "    print('Error, Std of mean particle radius cannot be 0 or negative, please reset!')\n",
    "Thick_std = Thick_std_.children[1].value * 10**(-6)   #m\n",
    "if Thick_std <= 0:\n",
    "    error = True\n",
    "    print('Error, Std of electrode thickness cannot be 0 or negative, please reset!')\n",
    "Rr_coef_std = Rr_coef_std_.children[1].value\n",
    "if Rr_coef_std <= 0:\n",
    "    error = True\n",
    "    print('Error, Std of reaction rate coefficient cannot be 0 or negative, please reset!')\n",
    "Dif_coef_std = Dif_coef_std_.children[1].value * 10**(-15)\n",
    "if Dif_coef_std <= 0:\n",
    "    error = True\n",
    "    print('Error, Std of diffusion coefficient cannot be 0 or negative, please reset!')\n",
    "Init_c_std = Init_c_std_.children[1].value\n",
    "if Init_c_std <= 0:\n",
    "    error = True\n",
    "    print('Error, Std of the initial lithium concentration cannot be 0 or negative, please reset!')\n",
    "Max_c_std = Max_c_std_.children[1].value\n",
    "if Max_c_std <= 0:\n",
    "    error = True\n",
    "    print('Error, Std of maximum lithium cannot be 0 or negative, please reset!')\n",
    "Vol_per_std = Vol_per_std_.children[1].value\n",
    "if Vol_per_std <= 0:\n",
    "    error = True\n",
    "    print('Error, Std of active material volume fraction cannot be 0 or negative, please reset!')\n",
    "Iapp_std = Iapp_std_.children[1].value\n",
    "if Iapp_std <= 0:\n",
    "    error = True\n",
    "    print('Error, Std of the applied current cannot be 0 or negative, please reset!')\n",
    "if error == True:\n",
    "    print('After the parameter has been changed, please rerun this cell.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "55822dbf",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# write UQ parameters to input file\n",
    "if UQ_.children[1].value == 'Yes':\n",
    "    if error == True:\n",
    "        print('Please adjust the std indicated by the error message from the previous cell before running this cell.')\n",
    "        raise StopExecution\n",
    "    rootgrp = Dataset('SPM_input.nc', 'r+', format='NETCDF4')\n",
    "    # no_samples\n",
    "    no_samples_dim = rootgrp.createDimension('no_samples_dim', 1)\n",
    "    no_samples = rootgrp.createVariable('no_samples', 'i4', ('no_samples_dim',))\n",
    "    no_samples.description = 'No. of samples'\n",
    "    no_samples.units = 'unitless'\n",
    "    no_samples[0] = No_samples\n",
    "    # Temp_std\n",
    "    if Temp_std != 0:\n",
    "        temp_std_dim = rootgrp.createDimension('temp_std_dim', 1)\n",
    "        temp_std = rootgrp.createVariable('temp_std', 'f8', ('temp_std_dim',))\n",
    "        temp_std.description = 'Std of temperature'\n",
    "        temp_std.units = 'K'\n",
    "        temp_std[0] = Temp_std\n",
    "    # Rad_std\n",
    "    if Rad_std != 0:\n",
    "        rad_std_dim = rootgrp.createDimension('rad_std_dim', 1)\n",
    "        rad_std = rootgrp.createVariable('rad_std', 'f8', ('rad_std_dim',))\n",
    "        rad_std.description = 'Std of mean particle radius'\n",
    "        rad_std.units = 'm'\n",
    "        rad_std[0] = Rad_std\n",
    "    # Thick_std\n",
    "    if Thick_std != 0:\n",
    "        thick_std_dim = rootgrp.createDimension('thick_std_dim', 1)\n",
    "        thick_std = rootgrp.createVariable('thick_std', 'f8', ('thick_std_dim',))\n",
    "        thick_std.description = 'Std of electrode thickness'\n",
    "        thick_std.units = 'm'\n",
    "        thick_std[0] = Thick_std\n",
    "    # Rr_std\n",
    "    if Rr_coef_std != 0:\n",
    "        rr_coef_std_dim = rootgrp.createDimension('rr_coef_std_dim', 1)\n",
    "        rr_coef_std = rootgrp.createVariable('rr_coef_std', 'f8', ('rr_coef_std_dim',))\n",
    "        rr_coef_std.description = 'Std of reaction rate coefficient'\n",
    "        rr_coef_std.units = '$Am^{-2}(m^3mol^{-1})^{1.5}$'\n",
    "        rr_coef_std[0] = Rr_coef_std\n",
    "    # Dif_coef\n",
    "    if Dif_coef_std != 0:\n",
    "        dif_coef_std_dim = rootgrp.createDimension('dif_coef_std_dim', 1)\n",
    "        dif_coef_std = rootgrp.createVariable('dif_coef_std', 'f8', ('dif_coef_std_dim',))\n",
    "        dif_coef_std.description = 'Std of diffusion coefficient'\n",
    "        dif_coef_std.units = '$m^2 s^{-1}$'\n",
    "        dif_coef_std[0] = Dif_coef_std\n",
    "    # Init_c_std\n",
    "    if Init_c_std != 0:\n",
    "        init_c_std_dim = rootgrp.createDimension('init_c_std_dim', 1)\n",
    "        init_c_std = rootgrp.createVariable('init_c_std', 'f8', ('init_c_std_dim',))\n",
    "        init_c_std.description = 'Std of initial concentration'\n",
    "        init_c_std.units = '$mol m^{-3}$'\n",
    "        init_c_std[0] = Init_c_std\n",
    "    # Max_c_std\n",
    "    if Max_c_std != 0:\n",
    "        max_c_std_dim = rootgrp.createDimension('max_c_std_dim', 1)\n",
    "        max_c_std = rootgrp.createVariable('max_c_std', 'f8', ('max_c_std_dim',))\n",
    "        max_c_std.description = 'Std of maximum concentration'\n",
    "        max_c_std.units = '$mol m^{-3}$'\n",
    "        max_c_std[0] = Max_c_std\n",
    "    # Vol_per_std\n",
    "    if Vol_per_std != 0:\n",
    "        vol_per_std_dim = rootgrp.createDimension('vol_per_std_dim', 1)\n",
    "        vol_per_std = rootgrp.createVariable('vol_per_std', 'f8', ('vol_per_std_dim',))\n",
    "        vol_per_std.description = 'Std of active material volume fraction'\n",
    "        vol_per_std.units = '%'\n",
    "        vol_per_std[0] = Vol_per_std\n",
    "    # Iapp_std\n",
    "    if Iapp_std != 0:\n",
    "        iapp_std_dim = rootgrp.createDimension('iapp_std_dim', 1)\n",
    "        iapp_std = rootgrp.createVariable('iapp_std', 'f8', ('iapp_std_dim',))\n",
    "        iapp_std.description = 'Std of initial concentration'\n",
    "        iapp_std.units = 'A'\n",
    "        iapp_std[0] = Iapp_std\n",
    "    rootgrp.close()"
   ]
  }
 ],
 "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}