TUTORIAL_branching_ratio_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reaction stoichiometry (branching ratio) example\n",
    "\n",
    "In this notebook, a branching ratio will be applied to the reaction: `A + B --> C + D` to afford `A + B --> (0.75)C + (0.25)D` with the branching ratios in brackets. \n",
    "\n",
    "Model construction is dealt with in the crash course notebook. Here we will focus on incorporating a branching ratio into the reaction scheme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# importing the necessaries\n",
    "import numpy as np\n",
    "import multilayerpy\n",
    "import multilayerpy.build as build \n",
    "import multilayerpy.simulate as simulate\n",
    "import multilayerpy.optimize as optimize\n",
    "\n",
    "\n",
    "# useful to know what version of multilayerpy we are using\n",
    "print(multilayerpy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the reaction scheme with a branching ratio\n",
    "The first thing to do is define the reaction scheme. In this case, `A + B --> (0.75)C + (0.25)D`. There are 4 components. \n",
    "\n",
    "It is possible to define branching ratos for both reactants and products. In this case, we are only interested in the products. In addition to the `product_tuple_list` supplied to the `ReactionScheme` object, we will define a `product_stoich` (product stoichiometry) tuple list in a similar fashion:\n",
    "\n",
    "`product_stoich = [(0.75,0.25)]`\n",
    "\n",
    "This states that for reaction 1, the products 3 and 4 (`C` and `D`) have stoichiometric coefficients (branching ratios) of 0.75 and 0.25, respectively.\n",
    "\n",
    "This is then supplied to the `ReactionScheme` object as outlined in the code cell below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the ModelType class\n",
    "from multilayerpy.build import ModelType\n",
    "\n",
    "# import the ReactionScheme class\n",
    "from multilayerpy.build import ReactionScheme\n",
    "\n",
    "# define the model type (KM-SUB in this case) and geometry (spherical or film)\n",
    "mod_type = ModelType('km-sub','spherical')\n",
    "\n",
    "# build the reaction tuple list, in this case only 1 tuple in the list (for 1 reaction)\n",
    "reaction_tuple_list = [(1,2)]\n",
    "\n",
    "\n",
    "# build the product tuple list\n",
    "product_tuple_list = [(3,4)]\n",
    "\n",
    "# now add the product stoichiometry **********\n",
    "product_stoich = [(0.75,0.25)]\n",
    "\n",
    "# now construct the reaction scheme with the product stoichiometry\n",
    "# we can give it a name and define the nuber of components as below\n",
    "reaction_scheme = ReactionScheme(mod_type,\n",
    "                                 name='AB_CD',\n",
    "                                 reactants=reaction_tuple_list,\n",
    "                                 products=product_tuple_list,\n",
    "                                 product_stoich=product_stoich)\n",
    "\n",
    "# let's print out a representation of the reaction scheme\n",
    "reaction_scheme.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model components\n",
    "\n",
    "# import ModelComponent class\n",
    "from multilayerpy.build import ModelComponent\n",
    "\n",
    "# making model components\n",
    "\n",
    "# A\n",
    "A = ModelComponent(1,reaction_scheme,name='A')\n",
    "\n",
    "# B (gas)\n",
    "B = ModelComponent(2,reaction_scheme,gas=True,name='B') \n",
    "\n",
    "# C\n",
    "C = ModelComponent(3,reaction_scheme, name='C')\n",
    "\n",
    "# D\n",
    "D = ModelComponent(4,reaction_scheme, name='D')\n",
    "\n",
    "# collect into a dictionary\n",
    "model_components_dict = {'1':A,\n",
    "                        '2':B,\n",
    "                        '3':C,\n",
    "                        '4':D}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Diffusion\n",
    "\n",
    "# import DiffusionRegime class\n",
    "from multilayerpy.build import DiffusionRegime\n",
    "\n",
    "# making the diffusion dictionary\n",
    "diff_dict = None  \n",
    "\n",
    "# make diffusion regime\n",
    "diff_regime = DiffusionRegime(mod_type,model_components_dict,diff_dict=diff_dict)\n",
    "\n",
    "# call it to build diffusion code ready for the builder\n",
    "diff_regime()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the model\n",
    "\n",
    "# import ModelBuilder class\n",
    "from multilayerpy.build import ModelBuilder\n",
    "\n",
    "# create the model object\n",
    "model = ModelBuilder(reaction_scheme,model_components_dict,diff_regime)\n",
    "\n",
    "# build the model. Will save a file, won't include the date in the model filename\n",
    "model.build(date=False)\n",
    "\n",
    "# print out the parameters required for the model to run (useful for constructing the parameter dictionary later)\n",
    "print(model.req_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run simulation and plot\n",
    "\n",
    "# import the Simulate class\n",
    "from multilayerpy.simulate import Simulate\n",
    "\n",
    "# import the Parameter class\n",
    "from multilayerpy.build import Parameter\n",
    "\n",
    "# make the parameter dictionary\n",
    "param_dict = {'delta_3':Parameter(1e-7),  # cm\n",
    "              'alpha_s_0_2':Parameter(4.2e-4),  \n",
    "              'delta_2':Parameter(0.4e-7),  # cm\n",
    "              'Db_2':Parameter(1e-5),  # cm2 s-1\n",
    "              'delta_1':Parameter(0.8e-7),  # cm\n",
    "              'Db_1':Parameter(1e-10),  # cm2 s-1\n",
    "              'Db_3':Parameter(1e-10),  # cm2 s-1\n",
    "              'k_1_2':Parameter(1.7e-15),  # cm3 s-1\n",
    "              'H_2':Parameter(4.8e-4),  # mol cm-3 atm-1\n",
    "              'Xgs_2': Parameter(7.0e13),  # cm-3\n",
    "              'Td_2': Parameter(1e-2),  # s\n",
    "              'w_2':Parameter(3.6e4),  # cm s-1\n",
    "              'T':Parameter(298.0),  # K\n",
    "              'k_1_2_surf':Parameter(6.0e-12),  # cm2 s-1\n",
    "             'delta_4':Parameter(1e-7),  # cm\n",
    "             'Db_4':Parameter(1e-10)}  # cm2 s-1\n",
    "\n",
    "\n",
    "# make the simulate object with the model and parameter dictionary\n",
    "sim = Simulate(model,param_dict)\n",
    "\n",
    "# define required parameters\n",
    "n_layers = 10\n",
    "rp = 0.2e-4 # radius in cm\n",
    "time_span = [0,40] # in s (times between which to run the model)\n",
    "n_time = 999 # number of timepoints to save to output\n",
    "\n",
    "#spherical V and A\n",
    "# use simulate.make_layers function\n",
    "V, A, layer_thick = simulate.make_layers(mod_type,n_layers,rp)\n",
    "\n",
    "# initial conc. of everything\n",
    "\n",
    "bulk_conc_dict = {'1':1.21e21,'2':0,'3':0,'4':0} # key=model component number, value=bulk conc, units: cm-3\n",
    "surf_conc_dict = {'1':9.68e13,'2':0,'3':0,'4':0} # key=model component number, value=surf conc, units: cm-2\n",
    "\n",
    "y0 = simulate.initial_concentrations(mod_type,bulk_conc_dict,surf_conc_dict,n_layers) \n",
    "    \n",
    "# now run the model\n",
    "output = sim.run(n_layers,rp,time_span,n_time,V,A,layer_thick,Y0=y0)\n",
    "\n",
    "%matplotlib inline\n",
    "# plot the model\n",
    "fig = sim.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "Have a play with the branching ratios in the reaction scheme. You will see different relative amounts of `C` and `D` in the model output. This can be done for many-reaction reaction schemes to increase complexity.\n",
    "\n",
    "You can also check out the model code outputted as `km-sub_AB_CD.py` and see where the branching ratios are applied in the differential equations defined for `C` and `D`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}