TUTORIAL_insoluble_monolayers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Insoluble monolayers\n",
    "\n",
    "Experiments can be carried out on insoluble monolayers of material. This notebook explains how one would create and fit a model to data from such experiments. The dataset used is from Woden et al. (2021) (https://doi.org/10.5194/acp-21-1325-2021).\n",
    "\n",
    "This dataset is from a monolayer of deuterated oleic acid on an aqueous subphase. We will assume that none of the products evaporate from the surface and keep them lumped together as a \"products\" component for simplicity. \n",
    "\n",
    "The reaction scheme is consistent with other MultilayerPy tutorials: `oleic acid + ozone --> products`"
   ]
  },
  {
   "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",
    "# check the version\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-SUB in this case) and geometry (film)\n",
    "mod_type = ModelType('km-sub','film')\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, only component 3 (products) is a product\n",
    "# a tuple with a single value inside is defined (value,)\n",
    "product_tuple_list = [(3,)]\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 a gas\n",
    "O3 = ModelComponent(2,reaction_scheme,gas=True,name='Ozone') \n",
    "\n",
    "# products\n",
    "prod = ModelComponent(3,reaction_scheme, name='Reaction products')\n",
    "\n",
    "# collect into a dictionary\n",
    "model_components_dict = {'1':OA,\n",
    "                        '2':O3,\n",
    "                        '3':prod}"
   ]
  },
  {
   "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import ModelBuilder class\n",
    "from multilayerpy.build import ModelBuilder\n",
    "\n",
    "# create the model object, ignore [1,2,3] etc at the end\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": "markdown",
   "metadata": {},
   "source": [
    "## Making an insoluble monolayer\n",
    "\n",
    "An insoluble monolayer can be made by setting the bulk diffusion coefficient of that component = 0.0 cm2 s-1. This means that that component does not diffuse at all in the bulk. We will nominally create 5 model bulk layers and a bulk film thickness of 1 µm to satisfy the model building process. However, there is essentially no exchange of material between the bulk and surface layers. Henry's law coefficient is also set to 0.0.\n",
    "\n",
    "Ozone concentration was 323 ppb (323 x 2.46e10 cm-3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "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",
    "# make the parameter dictionary\n",
    "# SETTING BULK DIFFUSION AND HENRY'S LAW PARAMETERS TO 0.0\n",
    "param_dict = {'delta_3':Parameter(1e-7),  # cm\n",
    "              'alpha_s_0_2':Parameter(1e-3,vary=True,bounds=(1e-4,1e-2)),\n",
    "              'delta_2':Parameter(0.4e-7),  # cm\n",
    "              'Db_2':Parameter(0.0),  # cm2 s-1\n",
    "              'delta_1':Parameter(0.8e-7),  # cm\n",
    "              'Db_1':Parameter(0.0),  # cm2 s-1\n",
    "              'Db_3':Parameter(0.0),  # cm2 s-1\n",
    "              'k_1_2':Parameter(0.0),  # cm3 s-1\n",
    "              'H_2':Parameter(0.0),  # mol cm-3 atm-1\n",
    "              'Xgs_2': Parameter(323.0 * 2.46e10),  # cm-3\n",
    "              'Td_2': Parameter(1e-7,vary=True,bounds=(1e-8,1e-5)),  # s\n",
    "              'w_2':Parameter(3.6e4),  # cm s-1\n",
    "              'T':Parameter(294.0),  # K\n",
    "              'k_1_2_surf':Parameter(2.2e-10,vary=True,bounds=(1e-11,1e-9))}  # cm2 s-1\n",
    "\n",
    "# import the data\n",
    "from multilayerpy.simulate import Data\n",
    "\n",
    "raw_data = np.genfromtxt('woden_etal_acp_data_uncert.csv',delimiter=',')\n",
    "\n",
    "data = Data(raw_data)\n",
    "\n",
    "# make the simulate object with the model and parameter dictionary\n",
    "sim = Simulate(model,param_dict,data=data)\n",
    "\n",
    "# define required parameters\n",
    "n_layers = 5 # nominal 5 bulk layers\n",
    "rp = 2e-4 # nominal film thickness in cm\n",
    "time_span = [0,800] # 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\n",
    "bulk_conc_dict = {'1':0.0,'2':0.0,'3':0.0} # key=model component number, value=bulk conc, units: cm-3\n",
    "\n",
    "# initial surf conc of oleic acid calculated from neutron reflectometry model fit\n",
    "surf_conc_dict = {'1':262571428571428.56,'2':0.0,'3':0.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(norm=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now optimise the model\n",
    "\n",
    "from time import time\n",
    "\n",
    "# import the optimize module and Optimizer object\n",
    "import multilayerpy.optimize\n",
    "from multilayerpy.optimize import Optimizer\n",
    "\n",
    "fitter = Optimizer(sim)\n",
    "\n",
    "res = fitter.fit(method='least_squares');\n",
    "\n",
    "\n",
    "fig = sim.plot(norm=True)"
   ]
  }
 ],
 "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
}