TUTORIAL_KM-GAP_model_creation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Making a KM-GAP model\n",
    "\n",
    "The kinetic multi-layer model of gas-particle interations in aerosols and clouds (KM-GAP) is related to the KM-SUB models used in the rest of these tutorials. One of the main differences is that it can incorporate changes in the size of model layers as a function of the number of molecules in each layer. So particle size and film thickness can be followed during a reaction (https://doi.org/10.5194/acp-12-2777-2012). \n",
    "\n",
    "These models are computationally more expensive and need careful consideration. Here is how to create a KM-GAP model using the `oleic acid + ozone --> nonanal + products`. Note that in this case, we will include the volatile product nonanal in the model system and watch it evaporate. \n",
    "\n",
    "The model building steps are essentially the same as for KM-SUB, with a few additions. "
   ]
  },
  {
   "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",
    "# useful to know what version of multilayerpy we are using\n",
    "print(multilayerpy.__version__)"
   ]
  },
  {
   "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-GAP in this case) and geometry (spherical or film)\n",
    "mod_type = ModelType('km-gap','spherical')\n",
    "\n",
    "# build the reaction tuple list, in this case only 1 tuple in the list (for 1 reaction)\n",
    "# component 1 (oleic acid) reacts with component 2 (ozone)\n",
    "reaction_tuple_list = [(1,2)]\n",
    "\n",
    "# build the product tuple list, component 3 (products) and 4 (nonanal)\n",
    "product_tuple_list = [(3,4)]\n",
    "\n",
    "# now construct the reaction scheme\n",
    "# we can give it a name and define the nuber of components as below\n",
    "reaction_scheme = ReactionScheme(mod_type,name='Oleic acid ozonolysis',\n",
    "                                                   reactants=reaction_tuple_list,\n",
    "                                                products=product_tuple_list)\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": [
    "# import ModelComponent class\n",
    "from multilayerpy.build import ModelComponent\n",
    "\n",
    "# making model components\n",
    "\n",
    "# oleic acid\n",
    "OA = ModelComponent(1,reaction_scheme,name='Oleic acid')\n",
    "\n",
    "# ozone, declare that it is in the gas phase\n",
    "O3 = ModelComponent(2,reaction_scheme,gas=True,name='Ozone') \n",
    "\n",
    "# products\n",
    "prod = ModelComponent(3,reaction_scheme, name='Reaction products')\n",
    "\n",
    "# nonanal\n",
    "nonanal = ModelComponent(4,reaction_scheme, name='Nonanal')\n",
    "\n",
    "# collect into a dictionary\n",
    "model_components_dict = {'1':OA,\n",
    "                        '2':O3,\n",
    "                        '3':prod,\n",
    "                        '4':nonanal}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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": [
    "# 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, don'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\n",
    "print(model.req_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the Simulate class\n",
    "from multilayerpy.simulate import Simulate\n",
    "\n",
    "# import the Parameter class\n",
    "from multilayerpy.build import Parameter\n",
    "\n",
    "# build up the parameter dictionary - note the large amount of model parameters required for KM-GAP\n",
    "param_dict = {'w_2':Parameter(3.6e4),  # cm s-1\n",
    "              'delta_1':Parameter(0.8e-7),  # cm\n",
    "              'p_4':Parameter(0.4 / 760.0),  # Pa\n",
    "              'w_4':Parameter(590.0),  # cm s-1\n",
    "              'Zgs_2':Parameter(7.0e13),  # cm-3\n",
    "              'Zgs_3':Parameter(0.0),  # cm-3\n",
    "              'Td_1':Parameter(1e6),  # s\n",
    "              'Td_4':Parameter(1e-2,vary=True,bounds=(1e-6,1e-1)),  # s\n",
    "              'delta_4':Parameter(0.4e-7),  # cm\n",
    "              'Td_2':Parameter(1e-2),  # s\n",
    "              'alpha_s_0_2':Parameter(4.2e-4,vary=True,bounds=(1e-4,1.0)),  \n",
    "              'Zgs_1':Parameter(0.0),  # cm-3\n",
    "              'alpha_s_0_1':Parameter(1.0),  \n",
    "              'p_2':Parameter(1.0),  # Pa\n",
    "              'delta_3':Parameter(0.4e-7),  # cm\n",
    "              'p_1':Parameter(0.0),  # Pa\n",
    "              'Td_3':Parameter(1e6),  # s\n",
    "              'Db_3':Parameter(1e-10),  # cm2 s-1\n",
    "              'alpha_s_0_4':Parameter(1.0),  \n",
    "              'k_1_2':Parameter(1.7e-15),  # cm3 s-1\n",
    "              'Db_2':Parameter(1e-5),  # cm2 s-1\n",
    "              'p_3':Parameter(0.0),  # Pa\n",
    "              'Zgs_4':Parameter(0.0),  # cm-3 \n",
    "              'w_3':Parameter(590.0),  # cm s-1\n",
    "              'w_1':Parameter(0.0),  # cm s-1\n",
    "              'T':Parameter(298.0),  # K\n",
    "              'Db_4':Parameter(1e-10),  # cm2 s-1\n",
    "              'Db_1':Parameter(1e-10),  # cm2 s-1\n",
    "              'alpha_s_0_3':Parameter(1.0),  \n",
    "              'k_1_2_surf':Parameter(6e-12),  # cm2 s-1\n",
    "              'delta_2':Parameter(0.4e-7)}  # cm\n",
    "\n",
    "\n",
    "# make the simulate object with the model and parameter dictionary\n",
    "sim = Simulate(model,param_dict)\n",
    "\n",
    "# import the data to fit to later\n",
    "from multilayerpy.simulate import Data\n",
    "raw_data = np.genfromtxt('ziemann_data.csv',delimiter=',')\n",
    "\n",
    "data = Data(raw_data)\n",
    "\n",
    "sim.data = data\n",
    "\n",
    "# define required parameters\n",
    "n_layers = 10\n",
    "rp = 0.2e-4 # radius in cm\n",
    "time_span = [0,40] # in s\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 (in cm-3 for bulk concentrations and cm-2 for surface concentrations)\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",
    "static_surf_conc_dict = {'1':9.68e13,'2':0,'3':0,'4':0} # units: cm-2\n",
    "\n",
    "# make initial concentrations but with additional keyword arguements required for KM-GAP\n",
    "y0 = simulate.initial_concentrations(mod_type,bulk_conc_dict,surf_conc_dict,n_layers,\n",
    "                                     parameter_dict=param_dict,\n",
    "                                     static_surf_conc_dict = static_surf_conc_dict,\n",
    "                                    V=V,A=A) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now run the model\n",
    "output = sim.run(n_layers,rp,time_span,n_time,V,A,layer_thick,Y0=y0)\n",
    "\n",
    "\n",
    "%matplotlib inline\n",
    "# plot the model\n",
    "fig = sim.plot(norm=True)\n",
    "\n",
    "# plot bulk concentrations\n",
    "sim.plot_bulk_concs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the particle radius as a function of time\n",
    "rp_t = sim.rp_vs_t()\n",
    "\n",
    "# plot rp vs t\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(output.t,rp_t*1e4)\n",
    "plt.ylabel('rp / µm')\n",
    "plt.xlabel('Time / s')\n",
    "plt.title('Particle radius vs time')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# model optimisation to the data\n",
    "\n",
    "from multilayerpy.optimize import Optimizer\n",
    "\n",
    "fitter = Optimizer(sim)\n",
    "\n",
    "res = fitter.fit(method='differential_evolution',weighted=False) #set weighted to False as there are no uncertainties\n",
    "\n",
    "fig = sim.plot(norm=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.plot_bulk_concs()\n",
    "\n",
    "# get the particle radius as a function of time again\n",
    "rp_t = sim.rp_vs_t()\n",
    "\n",
    "# plot rp vs t\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(output.t,rp_t*1e4) # convert cm to µm\n",
    "plt.ylabel('rp / µm')\n",
    "plt.xlabel('Time / s')\n",
    "plt.title('Particle radius vs time')\n",
    "plt.show()"
   ]
  }
 ],
 "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": 5
}