{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook shows how to treat surface reactions in reacting systems where diffusion is a much faster process in comparison to the chemical reactions (surface is treated as a \"homogenized\" concentration)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from collections import defaultdict\n",
    "from chempy import ReactionSystem, Species, Reaction\n",
    "from chempy.units import to_unitless, SI_base_registry, rescale, default_constants as const, default_units as u\n",
    "from chempy.kinetics.ode import (\n",
    "    _create_odesys as create_odesys,\n",
    ")\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ads_token = '(ads)'\n",
    "phases = {'(aq)': 0, ads_token: 1}\n",
    "ads_comp = -1\n",
    "vacancy = Species('vacancy(ads)', composition={ads_comp: 1},\n",
    "                  phase_idx=phases[ads_token], latex_name='vacancy(ads)')\n",
    "\n",
    "def substance_factory(name):\n",
    "    if name == vacancy.name:\n",
    "        return vacancy\n",
    "    s = Species.from_formula(name, phases=phases)\n",
    "    if name.endswith(ads_token):\n",
    "        s.composition[ads_comp] = 1\n",
    "    return s\n",
    "\n",
    "rsys = ReactionSystem.from_string(\"\"\"\n",
    "vacancy(ads) + H2O2 -> H2O2(ads); 'k_ads_H2O2'\n",
    "H2O2(ads) + vacancy(ads) -> 2 OH(ads); 'k_split_H2O2ads'\n",
    "OH(ads) + H2O2 -> H2O + HO2 + vacancy(ads); 'k_abst_OHads_H2O2'\n",
    "HO2 + HO2 -> O2 + H2O2; 'k_disprop_HO2'\n",
    "\"\"\", substance_factory=substance_factory)\n",
    "rsys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ureg = SI_base_registry.copy()\n",
    "ureg['length'] = u.dm\n",
    "odesys, odesys_extra = create_odesys(rsys, unit_registry=ureg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "powder_mass = 42*u.g\n",
    "specific_surf_area = 35*u.m**2/u.g  # e.g. from BET\n",
    "sample_volume = 250*u.cm3\n",
    "surf_volume_ratio = powder_mass*specific_surf_area/sample_volume\n",
    "\n",
    "adsorption_cross_section = 0.162*u.nm**2  # N2 at 77K on active carbon\n",
    "conc_vacancy = rescale(surf_volume_ratio/(adsorption_cross_section * const.Avogadro_constant), u.molar)\n",
    "conc_vacancy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c0 = defaultdict(lambda: 0*u.M, H2O2=.2*u.M, **{'vacancy(ads)': conc_vacancy})\n",
    "\n",
    "tend = 3600*u.s\n",
    "result, result_extra = odesys_extra['unit_aware_solve'](\n",
    "    (1e-7*u.s, tend), c0,\n",
    "    dict(\n",
    "        k_ads_H2O2=.1/u.M/u.s,\n",
    "        k_split_H2O2ads=.2/u.M/u.s,\n",
    "        k_abst_OHads_H2O2=.3/u.M/u.s,\n",
    "        k_disprop_HO2=.4/u.M/u.s\n",
    "    ), integrator='cvode'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(14, 4))\n",
    "result.plot(ax=axes[0])\n",
    "result.plot(ax=axes[1], xscale='log', yscale='log')\n",
    "for ax in axes:\n",
    "    #ax.set_ylabel('Concentration / M')\n",
    "    #ax.set_xlabel('Time / min')\n",
    "    ax.set_ylim(to_unitless([1e-6*u.molar, c0['H2O2']], result_extra['unit_conc']))\n",
    "    ax.set_xlim(to_unitless([1*u.s, tend], result_extra['unit_time']))"
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}