Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8efdbeef",
   "metadata": {},
   "source": [
    "# Tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "074ddec4",
   "metadata": {},
   "source": [
    "This notebook serves as a tutorial for how to use the Single Particle Model Code. First, begin by running the next cell in order to import all necessary packages for executing the code within this jupyter notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ab36f96",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import necessary packages\n",
    "%matplotlib tk\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "from ipywidgets import HBox, Label\n",
    "from netCDF4 import Dataset\n",
    "import netCDF4 as NC\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import sys\n",
    "import matplotlib as mpl\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib import animation\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from matplotlib.ticker import (StrMethodFormatter, AutoMinorLocator)\n",
    "from matplotlib.animation import FuncAnimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "984e9995",
   "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": "93f1d44f",
   "metadata": {},
   "source": [
    "## Compiling the Program"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95577249",
   "metadata": {},
   "source": [
    "__How to compile your code:__ \n",
    "\n",
    "The following cells are used to compile all relevant files by executing terminal commands within the directory this file is located. The first cell executes `make clean`. This cleans the directory and ensures your directory is fresh before compilation. The next cell runs the `make` command. This compiles all relevant files and produces the program that will be run.\n",
    "\n",
    "You should only ever run these cells once, unless you need to do a fresh install."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "653a2387",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# If errors are being displayed, run this cell to delete all files created by make. \n",
    "# Next, run the first and second cells of this section again\n",
    "!make clean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6587b7b",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Run this cell to compile the code and create the executables\n",
    "!make"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c91fc56",
   "metadata": {},
   "source": [
    "## Generating the Input File"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f03b54c5",
   "metadata": {},
   "source": [
    "__How to create input files to run the program:__\n",
    "\n",
    "Now we have created the program, you need an input file to put into the program. The following cells allow you to set the input parameters required for the half-SPM simulation. You can set the input parameters within our 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 is then read by the programme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa6e7362",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "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": "3e3a2d78",
   "metadata": {},
   "source": [
    "__How to set parameters to whatever you want:__\n",
    "\n",
    "While the sliders above are restrictive, 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": null,
   "id": "29766fd7",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "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": "42bafd66",
   "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": null,
   "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": null,
   "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": "d7ffdacf",
   "metadata": {},
   "source": [
    "## Running the Code"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22e274a0",
   "metadata": {},
   "source": [
    "__How to run the program:__\n",
    "\n",
    "The following cell executes `make exe` which runs the executable created earlier using our newly created input file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb6d9db1",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Run this cell to execute the code and generate data\n",
    "!make exe"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e38b14d",
   "metadata": {},
   "source": [
    "## Starting from a Checkpoint"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75786402",
   "metadata": {},
   "source": [
    "__How to continue a run for a few more steps:__\n",
    "\n",
    "While you have now created an output file, with data ready to be analysed, the simulation can also be continued from a checkpoint file that was created during the first run above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b498ae41",
   "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": null,
   "id": "7bad8246",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "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": null,
   "id": "a4f9e36f",
   "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 / Area_.children[1].value\n",
    "        rootgrp.close()\n",
    "else:\n",
    "        Checkpoint = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "495ce202",
   "metadata": {},
   "source": [
    "If you have selected the checkpoint restart, you can now run `make exe` to progress your simulation from the last simulation step."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9c45436",
   "metadata": {},
   "source": [
    "## Visualising your Results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e524fbf",
   "metadata": {},
   "source": [
    "__How to visualise the output data__\n",
    "\n",
    "Now that you have all your results in an output file, display an animation of the data that was generated. The cell below executes `make visual` to visualise the results. \n",
    "\n",
    "The cell displays the profile of Lithium concentration across the particle radius as a function of time, alongside a coloured contour plot of the concentration over time. It also displays a plot of the voltage over time (if you selected the option to calculate the voltage)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca9f70d0",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Creates a pop up window with the plotted data\n",
    "plt.close()\n",
    "# reading in the NetCDF data from SP_output.nc\n",
    "dat = NC.Dataset(\"SP_output.nc\", \"r\", format=\"NETCDF4\")\n",
    "do_volt = dat['volt_do'][:][0]\n",
    "dt = dat['dt'][:][0]\n",
    "\n",
    "#creating figure and axes dependent on whether voltage data is to be written or not\n",
    "if do_volt == 1:\n",
    "\tfig = plt.figure(figsize=(10,8))\n",
    "\tax1 = plt.subplot2grid((2,2),(0,0))\n",
    "\tax2 = plt.subplot2grid((2,2),(0,1))\n",
    "\tax3 = plt.subplot2grid((2,2),(1,0), colspan=2)\n",
    "else:\n",
    "\tfig = plt.figure(figsize=(10,4))\n",
    "\tax1 = plt.subplot2grid((1,2),(0,0))\n",
    "\tax2 = plt.subplot2grid((1,2),(0,1))\n",
    "\n",
    "\n",
    "# Subplot 1 - Animtion of lithium concentration \n",
    "# accessing the 2d concentration data and creating a time and space axis\n",
    "c = dat['conc'][:]\n",
    "\n",
    "t_steps = np.shape(c)[0]\n",
    "r_steps = dat['space_steps'][:][0]\n",
    "rad = dat['rad'][:][0]*10**(6)\n",
    "time_axis = np.linspace(0, t_steps-1, t_steps-1).astype(int)\n",
    "x_axis = np.linspace(0, rad, r_steps, endpoint=True)\n",
    "\n",
    "# animation of the concentration data \n",
    "# the main parts of the  animation code were taken from:\n",
    "# https://brushingupscience.com/2016/06/21/matplotlib-animations-the-easy-way/\n",
    "\n",
    "intervaltime = 0.5   #time between each animations step in ms\n",
    "\n",
    "# plotting the first graph which is a 2D line of lithium ion concentration\n",
    "# with respect to the particle radius\n",
    "line, = ax1.plot(x_axis, c[0])\n",
    "\n",
    "# function evolving the animation, which is called each frame \n",
    "# the x-axis is the particle radius, so only the x-axis data\n",
    "# i.e. the lithium concantration changes with simulation time\n",
    "def animate(t):\n",
    "    line.set_ydata(c[t,:])\n",
    "    return line,\n",
    "\n",
    "# creating the animation using matplotlib's FunAnimation\n",
    "# the arguments to be passed:\n",
    "# figure of the graph  - fig from fix,ax = plt.subplots())\n",
    "# function evolving the animation - animate(t)\n",
    "# interval time between each timesetep - intervaltime\n",
    "# frames, the array to iterate over - time_axis\n",
    "# blit - blitting set to true\n",
    "animation = FuncAnimation(fig, animate, interval=intervaltime, frames=time_axis, blit=True)\n",
    "\n",
    "# customise the graph with axes labels, major and minor ticks, labels etc.\n",
    "ax1.set_xlabel('Distance from particle centre [$\\mu m$]', size=8)\n",
    "ax1.xaxis.set_minor_locator(AutoMinorLocator())\n",
    "ax1.tick_params(axis='x', labelsize=7)\n",
    "ax1.set_ylabel('Lithium concentration $mol*m^{-3}$', size=8)\n",
    "ax1.set_ylim(np.min(c)-1*10**(-9), np.max(c)+1*10**(-9))\n",
    "ax1.yaxis.set_minor_locator(AutoMinorLocator())\n",
    "ax1.tick_params(axis='y', labelsize=7)\n",
    "ax1.ticklabel_format(axis='both', style=\"sci\", useMathText=True)\n",
    "ax1.xaxis.offsetText.set_fontsize(7)\n",
    "ax1.yaxis.offsetText.set_fontsize(7)\n",
    "ax1.set_title('Lithium Concentration Across Particle Radius', size=10, pad=15.0)\n",
    "\n",
    "# Subplot 2 - pcolorplot of concentration data\n",
    "\n",
    "# most of the required data is the same as in subplot 1\n",
    "# get the discretisation steps of the radius\n",
    "dr = rad/(r_steps-1)\n",
    " \n",
    "# reversed virdis colormap, values over the range are black, under goes white\n",
    "colour = cm.get_cmap('viridis_r').copy()\n",
    "colour.set_under(color='w')\n",
    "colour.set_over(color='k')\n",
    "max_colour = np.max(c)\n",
    "min_colour = np.min(c)\n",
    "colour_range = max_colour - min_colour\n",
    "ticklist = np.linspace(min_colour, max_colour, 6)\n",
    "\n",
    "#create grid for plotting\n",
    "#invert arrays to make sure smallest circle is on top\n",
    "rs = np.linspace(0,rad,r_steps,endpoint=True)\n",
    "rs = rs[::-1]\n",
    "cinv = c[:,::-1]\n",
    "\n",
    "#draw patches of circles\n",
    "#zorder to ensure smallest circle is on top\n",
    "patches = []\n",
    "for i, r in enumerate(rs):\n",
    "    circle = mpl.patches.Circle((0,0), r, zorder=i)\n",
    "    patches.append(circle)\n",
    "\n",
    "#set limits for colourmap and set colours of circles for initial plot\n",
    "#add circles to collection    \n",
    "p = mpl.collections.PatchCollection(patches, cmap=colour)\n",
    "p.set_clim([min_colour,max_colour])\n",
    "colours = np.array(cinv[0,:])\n",
    "p.set_array(colours)\n",
    "ax2.add_collection(p)\n",
    "\n",
    "#set figure labels\n",
    "ax2.set_xlim(-rad,rad)\n",
    "ax2.set_ylim(-rad,rad)\n",
    "ax2.set_xlabel('Distance from particle centre [$\\mu m$]', size=8)\n",
    "ax2.xaxis.set_minor_locator(AutoMinorLocator())\n",
    "ax2.tick_params(axis='x', labelsize=7)\n",
    "ax2.set_ylabel('Distance from particle centre [$\\mu m$]', size=8)\n",
    "ax2.yaxis.set_minor_locator(AutoMinorLocator())\n",
    "ax2.tick_params(axis='y', labelsize=7)\n",
    "ax2.ticklabel_format(axis='both', style=\"sci\", useMathText=True)\n",
    "ax2.xaxis.offsetText.set_fontsize(7)\n",
    "ax2.yaxis.offsetText.set_fontsize(7)\n",
    "ax2.set_title('Contour Plot of Concentration inside Particle', size=10, pad=15.0)\n",
    "cbar = ax2.figure.colorbar(p, ax=ax2, cmap=colour, ticks=ticklist)\n",
    "cbar.ax.tick_params(labelsize=7)\n",
    "cbar.ax.set_ylabel('Lithium concentration $mol*m^{-3}$', size=8)\n",
    "cbar.ax.yaxis.set_major_formatter(StrMethodFormatter(\"{x:.10f}\"))\n",
    "\n",
    "# animation function for the circle plot animation\n",
    "#only changes the colours of the circles without redrawing them\n",
    "# the concentration data for each timestep is feeded into a colour map\n",
    "# the inteval and frames are the same as in subplot 1 and blitting is set to true\n",
    "def animate_pcol(t):\n",
    "    colours = np.array(cinv[t,:])\n",
    "    p.set_array(colours)\n",
    "    \n",
    "    return p,\n",
    "\n",
    "animation_pcol = FuncAnimation(fig, animate_pcol, interval=intervaltime, frames=time_axis, blit=True)\n",
    "\n",
    "# Subplot 3 - plot of voltage output if do_volt set to true\n",
    "if do_volt == 1:\n",
    "\n",
    "\t# accessing the voltage data and setting the number of time steps as x-axis\n",
    "\tvolt = dat['volt'][:][:,0]\n",
    "\ttime = np.linspace(0, t_steps, t_steps)\n",
    "\n",
    "\t# plotting the 2D graph and customising it\n",
    "\tax3.plot(time*dt, volt)\n",
    "\tax3.set_xlabel('Time [s]', size=8)\n",
    "\tax3.xaxis.set_minor_locator(AutoMinorLocator())\n",
    "\tax3.tick_params(axis='x', labelsize=7)\n",
    "\tax3.set_ylabel('Voltage [V]', size=8)\n",
    "\tax3.yaxis.set_minor_locator(AutoMinorLocator())\n",
    "\tax3.tick_params(axis='y', labelsize=7)\n",
    "\tax3.ticklabel_format(axis='both', style=\"sci\", useMathText=True)\n",
    "\tax3.xaxis.offsetText.set_fontsize(7)\n",
    "\tax3.yaxis.offsetText.set_fontsize(7)\n",
    "\tax3.set_title('Voltage output', size=10, pad=15.0)\n",
    "\n",
    "\t# animation of an curser that moves with the time axis with the same interval and frames\n",
    "\t# as the other two animations\n",
    "\tvl = ax3.axvline(time[0]*dt, color='black', linestyle=':')\n",
    "\tax3.set_xlim()\n",
    "\n",
    "\tdef animate_t_bar(t):\n",
    "\t\tvl.set_xdata(time[t]*dt)\n",
    "\t\treturn vl,\n",
    "\n",
    "\tanimation_t_bar = FuncAnimation(fig, animate_t_bar, interval=intervaltime, frames=time_axis, blit=True)\n",
    "\n",
    "plt.draw()\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "dat.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8ca56d1",
   "metadata": {},
   "source": [
    "## Uncertainty Propagation and Sensitivity Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b450b2e",
   "metadata": {},
   "source": [
    "__How to setup sensitivity analysis and uncertainty propagation:__\n",
    "\n",
    "At this point you will have calculated your concentration data, and the corresponding voltages. You can now investigate the sensitivity of your voltage calculation with respect to the input parameters and calculate the uncertainty associated with the voltage calculation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13fbcbbf",
   "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 (be wary, the more samples you take, the longer it takes to calculate uncertainty, though the more accurate your error bar will be). You then need to enter in the standard deviations associated with your input parameters. We expect your standard deviations to have been taken from your experiment.\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 and in this directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08ec2394",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# option to conduct uncertainty propagation\n",
    "UQ_ = HBox([Label('Conduct uncertainty propagation'), widgets.Select(options=['Yes', 'No'], value='Yes', 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": null,
   "id": "8d97a584",
   "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": null,
   "id": "dc3c292c",
   "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()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b035fe16",
   "metadata": {},
   "source": [
    "Now you have added the required variables to the input file, you can now perform sensitivity analysis and uncertainty propagation which is explained below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93918d0f",
   "metadata": {},
   "source": [
    "## Sensitivity Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86b2ce8b",
   "metadata": {},
   "source": [
    "__How to find out what parameters affect my results the most:__\n",
    "\n",
    "___Problem:___ You want to quantify the uncertainty associated with your voltage calculation, however, due to your selected inputs, the code takes a long time to run. \n",
    "\n",
    "___Solution:___ By performing sensitivity analysis, you can see what parameters are contributing least to the uncertainty of your results. Thus you can set the standard deviation of these parameters to 0 and take fewer samples."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69e56ca5",
   "metadata": {},
   "source": [
    "By executing `make sensitive`, with your input file present, you can calculate and display the first order sensitivity analysis over time for your chosen inputs. To re-display the results, if you close the window, execute `make vis_sens`. The following cells perform these operations.\n",
    "\n",
    "If you wish to analyse the results in a different visualisation program you can find the results being plotted in \"__uq_code/data_store_sens/sens_data.csv__\" along with your original input file (\"__SPM_input_ori.nc__\") and the input parameters used in the sensitivity analysis (\"__inputs.csv__\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8021e334",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Perform sensitivity analysis, following cell visualises output\n",
    "!make sensitive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fef573ef",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Creates pop up window that visualises the sensitivity analysis\n",
    "#Set mean of variables here\n",
    "params = ['temp','rad','thick','rr_coef','dif_coef','init_c','max_c','vol_per','iapp']\n",
    "\n",
    "#Import the original parameters from the original input file and save to a vector\n",
    "dat_inp = NC.Dataset(\"uq_code/data_store_sens/SPM_input_ori.nc\", \"r\", format=\"NETCDF4\")\n",
    "\n",
    "mu = np.array([dat_inp['temp'][:][0],\n",
    "               dat_inp['rad'][:][0],\n",
    "               dat_inp['thick'][:][0],\n",
    "               dat_inp['rr_coef'][:][0],\n",
    "               dat_inp['dif_coef'][:][0],\n",
    "               dat_inp['init_c'][:][0],\n",
    "               dat_inp['max_c'][:][0],\n",
    "               dat_inp['vol_per'][:][0],\n",
    "               dat_inp['iapp'][:][0]])\n",
    "\n",
    "sim_steps = dat_inp['sim_steps'][:][0]\n",
    "dt = dat_inp['dt'][:][0]\n",
    "\n",
    "#Auxilliary 2D array of variable means, identical columns\n",
    "Mu_s = np.tile(mu.reshape((9,1)),(1,sim_steps))\n",
    "\n",
    "#Set value of perturbation (eps) here\n",
    "#Set h as a matrix for element-wise division, identical columns with elements = h_i\n",
    "eps = 1e-6\n",
    "h_i = (eps*mu).reshape((9,1))\n",
    "h = np.tile(h_i,(1, sim_steps))\n",
    "\n",
    "#Obtain mean voltage curve V_0 here, each row identical\n",
    "#9 rows representing variables\n",
    "#columns represent time steps\n",
    "V_0 = np.zeros((9,sim_steps))\n",
    "dat_mu = NC.Dataset(\"uq_code/data_store_sens/SP_output_0.nc\", \"r\", format=\"NETCDF4\")\n",
    "volt_0 = np.array(dat_mu['volt'][:][:,0])\n",
    "V_0 = np.tile(volt_0, (9, 1))\n",
    "\n",
    "#Obtain perturbed voltages here\n",
    "#9 rows representing variables\n",
    "#Columns representing timesteps\n",
    "Vs = np.zeros((9,sim_steps))\n",
    "\n",
    "for i in range(1, 10, 1):\n",
    "    dat = NC.Dataset(f\"uq_code/data_store_sens/SP_output_{i}.nc\", \"r\", format=\"NETCDF4\")\n",
    "    Vs[i-1, :] = np.array(dat['volt'][:][:,0])\n",
    "    dat.close()\n",
    "\n",
    "#Compute dV/dx as a function of time\n",
    "dV_dx = np.divide((Vs - V_0),h)\n",
    "\n",
    "#Compute absolute scaled sensitivity using element-wise multiplication\n",
    "#each column i of the matrix represents the sensitivities for each paramter at time t=i\n",
    "V_first_sensitivities = np.abs(Mu_s*dV_dx)\n",
    "\n",
    "#Put the data into a data frame, and save to .csv\n",
    "dat_fram = {'temp_sens': V_first_sensitivities[0, :],\n",
    "            'rad_sens': V_first_sensitivities[1, :],\n",
    "            'thick_sens': V_first_sensitivities[2, :],\n",
    "            'rr_coef_sens': V_first_sensitivities[3, :],\n",
    "            'dif_coef_sens': V_first_sensitivities[4, :],\n",
    "            'init_c_sens': V_first_sensitivities[5, :],\n",
    "            'max_c_sens': V_first_sensitivities[6, :],\n",
    "            'vol_per_sens': V_first_sensitivities[7, :],\n",
    "            'iapp_sens': V_first_sensitivities[8, :],\n",
    "            'temp_dvdx': dV_dx[0, :],\n",
    "            'rad_dvdx': dV_dx[1, :],\n",
    "            'thick_dvdx': dV_dx[2, :],\n",
    "            'rr_coef_dvdx': dV_dx[3, :],\n",
    "            'dif_coef_dvdx': dV_dx[4, :],\n",
    "            'init_c_dvdx': dV_dx[5, :],\n",
    "            'max_c_dvdx': dV_dx[6, :],\n",
    "            'vol_per_dvdx': dV_dx[7, :],\n",
    "            'iapp_dvdx': dV_dx[8, :]}\n",
    "\n",
    "cols = ['temp_sens', 'rad_sens', 'thick_sens', 'rr_coef_sens', 'dif_coef_sens', 'init_c_sens', 'max_c_sens', 'vol_per_sens', 'iapp_sens',\n",
    "        'temp_dvdx', 'rad_dvdx', 'thick_dvdx', 'rr_coef_dvdx', 'dif_coef_dvdx', 'init_c_dvdx', 'max_c_dvdx', 'vol_per_dvdx', 'iapp_dvdx']\n",
    "df = pd.DataFrame(dat_fram, columns=cols)\n",
    "df.to_csv('sens_data.csv', index=False)\n",
    "\n",
    "\n",
    "#plot as animation\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "\n",
    "x = np.arange(0, 9, 1)\n",
    "\n",
    "#time between frames\n",
    "intervaltime = 0.5\n",
    "\n",
    "#plot first frame\n",
    "line, = plt.plot(x, V_first_sensitivities[:, 0], 'o-')\n",
    "\n",
    "#definition of animation\n",
    "def animate(i):\n",
    "    line.set_ydata(V_first_sensitivities[:, i])\n",
    "    time.set_text(('T=')+str(i*dt)+(' s'))\n",
    "    return line, time,\n",
    "\n",
    "#Design plot\n",
    "ax.set_xticks(x)\n",
    "ax.set_xticklabels(params, rotation=45)\n",
    "ax.set_ylabel('Absolute scaled sensitivity of $V$')\n",
    "ax.set_xlabel('Parameter')\n",
    "ax.set_ylim(np.min(V_first_sensitivities)+10**(-9)-0.1, np.max(V_first_sensitivities)+1)\n",
    "ax.grid()\n",
    "ax.set_title('First Order Sensitivities Over Time')\n",
    "\n",
    "time = ax.text(0.1,0.85, \"\", bbox={'facecolor':'w', 'alpha':0.5, 'pad':5}, transform=ax.transAxes)\n",
    "\n",
    "#plot animation\n",
    "animate_sensitivities = FuncAnimation(fig, animate,interval=10, frames=range(1,sim_steps), blit=True)\n",
    "\n",
    "dat_inp.close()\n",
    "dat_mu.close()\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f40ee9c",
   "metadata": {},
   "source": [
    "## Uncertainty Propagation\n",
    "\n",
    "### Using Sensitivity Analysis to Approximate Uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b51fed68",
   "metadata": {},
   "source": [
    "__How to get a quick approximation of uncertainty:__\n",
    "\n",
    "Provided you have just performed the sensitivity analysis above, you can get an approximation of the uncertainty using the standard deviations and sensitivity analysis data by executing `make uncer_from_sens`. This just gives a quick approximation of the voltage uncertainty. If you close the window and want to re-display execute `make vis_uncer_from_sens`. The following cells run these commands.\n",
    "\n",
    "If you want the data for analysis in a different program, the standard deviations of the voltage at each time step can be found in \"__uq_code/data_store_sens/std_V_dat.csv__\", along with the mean voltage in \"__SP_output_0.nc__\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33ac5c93",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Perform approximate uncertainty quantification using standard deviations and sensitivity analysis\n",
    "!make uncer_from_sens"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "941c8989",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Create pop up widow that visualises approximate uncertainty quatification\n",
    "#Get the results from the mean parameters and save to an array\n",
    "dat_mu = NC.Dataset(f\"uq_code/data_store_sens/SP_output_0.nc\", \"r\", format=\"NETCDF4\")\n",
    "volt_dat_mu = np.array(dat_mu['volt'][:][:,0])\n",
    "\n",
    "#Get the number of time steps and the step length in order to create a time x axis in seconds\n",
    "t_steps = dat_mu['sim_steps'][:][0]\n",
    "dt = dat_mu['dt'][:][0]\n",
    "x = (np.linspace(0, t_steps, t_steps).astype(int))*dt\n",
    "\n",
    "#Plot the mean voltage\n",
    "plt.plot(x, volt_dat_mu, color='black', label='Mean Voltage')\n",
    "plt.title('Uncertainty in Voltage Calculation', pad=15.0)\n",
    "plt.xlabel('Time [s]')\n",
    "plt.ylabel('Voltage [V]')\n",
    "plt.grid()\n",
    "\n",
    "std_V_dat = pd.read_csv('uq_code/data_store_sens/std_V_dat.csv')\n",
    "\n",
    "std_V = np.array(std_V_dat['std_V'][:])\n",
    "low_b = volt_dat_mu-(2*std_V)\n",
    "up_b = volt_dat_mu+(2*std_V)\n",
    "plt.fill_between(x, low_b, up_b, alpha=0.2, label='95% Confidence (SDs)', color='C3')\n",
    "plt.plot(x, low_b, alpha=0.5, color='C3')\n",
    "plt.plot(x, up_b, alpha=0.5, color='C3')\n",
    "\n",
    "#Display the plot\n",
    "plt.legend()\n",
    "\n",
    "dat_mu.close()\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6644e7e",
   "metadata": {},
   "source": [
    "### Using Random Latin Hypercube Sampling to Approximate Uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "781cf02e",
   "metadata": {},
   "source": [
    "Now that you know what parameters are the most important from the sensitivity analysis, and you want to get accurate error bars for the uncertainty they create, you can perform random Latin hypercube sampling using for the number of samples you selected by executing `make uncertain`. If you close the window, and want to re-display, you can execute `make vis_uncer`. \n",
    "\n",
    "This will give you accurate error bars provided you entered enough samples.\n",
    "\n",
    "The following cells execute these commands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8f467a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "!make uncertain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f3e85f5",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Creates pop up window that visualises the uncertainty quantification\n",
    "#Read the input parameters from the original input file\n",
    "inp_dat = NC.Dataset(f\"uq_code/data_store_up/SPM_input_ori.nc\", \"r\", format=\"NETCDF4\")\n",
    "num_samps = inp_dat['no_samples'][:][0]\n",
    "\n",
    "#Get the results from the mean parameters and save to an array\n",
    "dat_mu = NC.Dataset(f\"uq_code/data_store_up/SP_output_0.nc\", \"r\", format=\"NETCDF4\")\n",
    "volt_dat_mu = np.array(dat_mu['volt'][:][:,0])\n",
    "\n",
    "#Get the number of time steps and the step length in order to create a time x axis in seconds\n",
    "t_steps = dat_mu['sim_steps'][:][0]\n",
    "dt = dat_mu['dt'][:][0]\n",
    "x = (np.linspace(0, t_steps, t_steps).astype(int))*dt\n",
    "\n",
    "#Create an array to store all the voltage data from each run\n",
    "volt_array = np.empty((t_steps, (num_samps-1)))\n",
    "\n",
    "#Open the files from the data base, extract the voltage and save to the array\n",
    "for i in range(1, num_samps, 1):\n",
    "    dat = NC.Dataset(f\"uq_code/data_store_up/SP_output_{i}.nc\", \"r\", format=\"NETCDF4\")\n",
    "    volt_array[:, (i-1)] = np.array(dat['volt'][:][:,0])\n",
    "    dat.close()\n",
    "\n",
    "#Calculate the 2.5th and 97.5th percentile of the voltage at each time step to get a 95% confidence\n",
    "ana = np.percentile(volt_array, [2.5, 97.5], axis=1)\n",
    "\n",
    "#Plot the range onto a graph\n",
    "plt.fill_between(x, ana[0, :], ana[1, :], alpha=0.2, label='95% Confidence (random samps)', color='C0')\n",
    "plt.plot(x, ana[0, :], color='C0', alpha=0.5)\n",
    "plt.plot(x, ana[1, :], color='C0', alpha=0.5)\n",
    "\n",
    "#Save the data to a data frame, and a .csv file for the user to be able to move to a different plotting software\n",
    "dat_fram = {'mean volt': volt_dat_mu,\n",
    "            '5th percentile': ana[0,:],\n",
    "            '95th percentile': ana[1,:]}\n",
    "df = pd.DataFrame(dat_fram, columns=['mean', '5th percentile', '95th percentile'])\n",
    "df.to_csv('voltage_confidence_up.csv', index=False)\n",
    "\n",
    "#Plot the mean voltage\n",
    "plt.plot(x, volt_dat_mu, color='black', label='Mean Voltage')\n",
    "plt.title('Uncertainty in Voltage Calculation', pad=15.0)\n",
    "plt.xlabel('Time [s]')\n",
    "plt.ylabel('Voltage [V]')\n",
    "plt.grid()\n",
    "\n",
    "#If the user has performed the sensitivity analysis, and calculated uncertainty using dV_dx (see generate_inp_params.py)...\n",
    "#...They can plot this data on top of the other uncertainty propagation data \n",
    "if (sys.argv[1] == 'True'):\n",
    "    std_V_dat = pd.read_csv('./std_V_dat.csv')\n",
    "\n",
    "    std_V = np.array(std_V_dat['std_V'][:])\n",
    "    low_b = volt_dat_mu-(2*std_V)\n",
    "    up_b = volt_dat_mu+(2*std_V)\n",
    "    plt.fill_between(x, low_b, up_b, alpha=0.2, label='95% Confidence (SDs)', color='C3')\n",
    "    plt.plot(x, low_b, alpha=0.5, color='C3')\n",
    "    plt.plot(x, up_b, alpha=0.5, color='C3')\n",
    "\n",
    "#Display the plot\n",
    "plt.legend()\n",
    "\n",
    "inp_dat.close()\n",
    "dat_mu.close()\n",
    "\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48cedd5b",
   "metadata": {},
   "source": [
    "### Comparing Estimates"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cca3e2d6",
   "metadata": {},
   "source": [
    "To plot both the Latin hypercube random sampling estimate of the uncertainty, and the first order sensitivity analysis approximation of the uncertainty you can execute the command `make vis_sens_uncer`. \n",
    "\n",
    "The cell below executes this command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abf4f9a1",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Creates pop up window that visualises the uncertainty quantification and sensitivity analysis\n",
    "#Read the input parameters from the original input file\n",
    "inp_dat = NC.Dataset(f\"uq_code/data_store_up/SPM_input_ori.nc\", \"r\", format=\"NETCDF4\")\n",
    "num_samps = inp_dat['no_samples'][:][0]\n",
    "\n",
    "#Get the results from the mean parameters and save to an array\n",
    "dat_mu = NC.Dataset(f\"uq_code/data_store_up/SP_output_0.nc\", \"r\", format=\"NETCDF4\")\n",
    "volt_dat_mu = np.array(dat_mu['volt'][:][:,0])\n",
    "\n",
    "#Get the number of time steps and the step length in order to create a time x axis in seconds\n",
    "t_steps = dat_mu['sim_steps'][:][0]\n",
    "dt = dat_mu['dt'][:][0]\n",
    "x = (np.linspace(0, t_steps, t_steps).astype(int))*dt\n",
    "\n",
    "#Create an array to store all the voltage data from each run\n",
    "volt_array = np.empty((t_steps, (num_samps-1)))\n",
    "\n",
    "#Open the files from the data base, extract the voltage and save to the array\n",
    "for i in range(1, num_samps, 1):\n",
    "    dat = NC.Dataset(f\"uq_code/data_store_up/SP_output_{i}.nc\", \"r\", format=\"NETCDF4\")\n",
    "    volt_array[:, (i-1)] = np.array(dat['volt'][:][:,0])\n",
    "    dat.close()\n",
    "\n",
    "#Calculate the 2.5th and 97.5th percentile of the voltage at each time step to get a 95% confidence\n",
    "ana = np.percentile(volt_array, [2.5, 97.5], axis=1)\n",
    "\n",
    "#Plot the range onto a graph\n",
    "plt.fill_between(x, ana[0, :], ana[1, :], alpha=0.2, label='95% Confidence (random samps)', color='C0')\n",
    "plt.plot(x, ana[0, :], color='C0', alpha=0.5)\n",
    "plt.plot(x, ana[1, :], color='C0', alpha=0.5)\n",
    "\n",
    "#Save the data to a data frame, and a .csv file for the user to be able to move to a different plotting software\n",
    "dat_fram = {'mean volt': volt_dat_mu,\n",
    "            '5th percentile': ana[0,:],\n",
    "            '95th percentile': ana[1,:]}\n",
    "df = pd.DataFrame(dat_fram, columns=['mean', '5th percentile', '95th percentile'])\n",
    "df.to_csv('voltage_confidence_up.csv', index=False)\n",
    "\n",
    "#Plot the mean voltage\n",
    "plt.plot(x, volt_dat_mu, color='black', label='Mean Voltage')\n",
    "plt.title('Uncertainty in Voltage Calculation', pad=15.0)\n",
    "plt.xlabel('Time [s]')\n",
    "plt.ylabel('Voltage [V]')\n",
    "plt.grid()\n",
    "\n",
    "#If the user has performed the sensitivity analysis, and calculated uncertainty using dV_dx (see generate_inp_params.py)...\n",
    "#...They can plot this data on top of the other uncertainty propagation data \n",
    "std_V_dat = pd.read_csv('uq_code/data_store_sens/std_V_dat.csv')\n",
    "\n",
    "std_V = np.array(std_V_dat['std_V'][:])\n",
    "low_b = volt_dat_mu-(2*std_V)\n",
    "up_b = volt_dat_mu+(2*std_V)\n",
    "plt.fill_between(x, low_b, up_b, alpha=0.2, label='95% Confidence (SDs)', color='C3')\n",
    "plt.plot(x, low_b, alpha=0.5, color='C3')\n",
    "plt.plot(x, up_b, alpha=0.5, color='C3')\n",
    "\n",
    "#Display the plot\n",
    "plt.legend()\n",
    "\n",
    "inp_dat.close()\n",
    "dat_mu.close()\n",
    "\n",
    "plt.show()\n",
    "\n"
   ]
  }
 ],
 "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
}