Merge pull request #96 from bjodah/surf

Homogenous treatment of multiphase systems

CHANGES.rst

@@ -2,6 +2,10 @@ v0.7.x (unreleased)
 ===================
 - Drop official support for Python 2.7
 
+v0.6.x (unreleased)
+===================
+- Upper limit on pyodesys (<0.12), chempy 0.7.x will use pyodesys>=0.12
+
 v0.6.7
 ======
 - Updated manuscript for JOSS.

chempy/kinetics/ode.py

@@ -20,7 +20,7 @@
 
 from ..units import (
     to_unitless, get_derived_unit, rescale, magnitude, unitless_in_registry,
-    default_units as u
+    default_unit_in_registry, default_units as u
 )
 from ..util.pyutil import deprecated
 from ..util._expr import Expr
@@ -451,7 +451,7 @@ def chained_parameter_variation(odesys, durations, init_conc, varied_params, def
 
 
 def _create_odesys(rsys, substance_symbols=None, parameter_symbols=None, pretty_replace=lambda x: x,
-                   backend=None, SymbolicSys=None, time_symbol=None):
+                   backend=None, SymbolicSys=None, time_symbol=None, unit_registry=None):
     """ This will be a simpler version of get_odesys without the unit handling code.
     The motivation is to reduce complexity (the code of get_odesys is long with multiple closures).
 
@@ -490,6 +490,13 @@ def _create_odesys(rsys, substance_symbols=None, parameter_symbols=None, pretty_
         if list(substance_symbols) != list(rsys.substances):
             raise ValueError("substance_symbols needs to have same (oredered) keys as rsys.substances")
 
+    if parameter_symbols is None:
+        keys = []
+        for rxn in rsys.rxns:
+            key, = rxn.param.unique_keys
+            keys.append(key)
+        parameter_symbols = OrderedDict([(key, backend.Symbol(key)) for key in keys])
+
     if not isinstance(parameter_symbols, OrderedDict):
         raise ValueError("parameter_symbols needs to be an OrderedDict")
 
@@ -515,14 +522,47 @@ def _create_odesys(rsys, substance_symbols=None, parameter_symbols=None, pretty_
         par_by_name=True
     )
 
+    validate = partial(_validate, rsys=rsys, symbols=symbols, odesys=odesys, backend=backend)
     return odesys, {
         'symbols': symbols,
-        'validate': partial(_validate, rsys=rsys, symbols=symbols,
-                            odesys=odesys, backend=backend),
+        'validate': validate,
+        'unit_aware_solve': _mk_unit_aware_solve(odesys, unit_registry, validate=validate) if unit_registry else None
     }
 
 
-def _validate(conditions, rsys, symbols, odesys, backend=None, transform=None):
+def _mk_dedim(unit_registry):
+    unit_time = get_derived_unit(unit_registry, 'time')
+    unit_conc = get_derived_unit(unit_registry, 'concentration')
+
+    def dedim_tcp(t, c, p, param_unit=lambda k, v: default_unit_in_registry(v, unit_registry)):
+        _t = to_unitless(t, unit_time)
+        _c = to_unitless(c, unit_conc)
+        _p, pu = {}, {}
+        for k, v in p.items():
+            pu[k] = param_unit(k, v)
+            _p[k] = to_unitless(v, pu[k])
+        return (_t, _c, _p), dict(unit_time=unit_time, unit_conc=unit_conc, param_units=pu)
+
+    return locals()
+
+
+def _mk_unit_aware_solve(odesys, unit_registry, validate):
+
+    dedim_ctx = _mk_dedim(unit_registry)
+
+    def solve(t, c, p, **kwargs):
+        for name in odesys.names:
+            c[name]  # to e.g. populate defaultdict
+        validate(dict(**c, **p))
+        tcp, dedim_extra = dedim_ctx['dedim_tcp'](t, c, p)
+        result = odesys.integrate(*tcp, **kwargs)
+        result.xout = result.xout * dedim_extra['unit_time']
+        result.yout = result.yout * dedim_extra['unit_conc']  # assumes only concentrations in y
+        return result, dedim_extra
+    return solve
+
+
+def _validate(conditions, rsys, symbols, odesys, backend=None, transform=None, ignore=('time',)):
     """ For use with create_odesys
 
     Parameters
@@ -547,7 +587,7 @@ def _validate(conditions, rsys, symbols, odesys, backend=None, transform=None):
 
     if transform is None:
         if backend.__name__ != 'sympy':
-            warnings.warn("Backend does not seem to be SymPy, please provide your own transform function.")
+            warnings.warn("Backend != SymPy, provide your own transform function.")
 
         def transform(arg):
             expr = backend.logcombine(arg, force=True)
@@ -566,6 +606,6 @@ def transform(arg):
         seen = [b or a in expr.free_symbols for b, a in zip(seen, args)]
     not_seen = [a for s, a in zip(seen, args) if not s]
     for k in conditions:
-        if k not in odesys.param_names and k not in odesys.names:
+        if k not in odesys.param_names and k not in odesys.names and k not in ignore:
             raise KeyError("Unknown param: %s" % k)
     return {'not_seen': not_seen, 'rates': rates}

chempy/kinetics/tests/test_ode.py

@@ -24,7 +24,7 @@
 from ..arrhenius import ArrheniusParam
 from ..rates import Arrhenius, MassAction, Radiolytic, RampedTemp
 from .._rates import ShiftedTPoly
-from ..ode import get_odesys, chained_parameter_variation
+from ..ode import get_odesys, chained_parameter_variation, _create_odesys as create_odesys
 from ..integrated import dimerization_irrev, binary_rev
 
 
@@ -910,3 +910,34 @@ def check(rsys, params):
     check(rsys3, dict(temperature=T_K*u.K,
                       dHf=dHf*u.J/u.mol, dSf=dSf*u.J/u.mol/u.K,
                       dHb=dHb*u.J/u.mol, dSb=dSb*u.J/u.mol/u.K))
+
+
+@requires('numpy', 'pyodesys', 'sympy', 'pycvodes')
+def test_create_odesys():
+    rsys = ReactionSystem.from_string("""
+    A -> B; 'k1'
+    B + C -> P; 'k2'
+    """, substance_factory=Substance)
+
+    odesys, odesys_extra = create_odesys(rsys, unit_registry=SI_base_registry)
+
+    tend = 10
+    tend_units = tend*u.s
+    c0 = {'A': 1e-6, 'B': 0, 'C': 1, 'P': 0}
+    p = dict(k1=3, k2=4)
+    p_units = {'k1': p['k1']/u.s, 'k2': p['k2']/u.M/u.s}
+    c0_units = {k: v*u.molar for k, v in c0.items()}
+    validation = odesys_extra['validate'](dict(**c0_units, **p_units))
+    P, = validation['not_seen']
+    assert P.name == 'P'
+    ref_rates = {'A': -p_units['k1']*c0_units['A'], 'P': p_units['k2']*c0_units['B']*c0_units['C']}
+    ref_rates['B'] = -ref_rates['A'] - ref_rates['P']
+    ref_rates['C'] = -ref_rates['P']
+    assert validation['rates'] == ref_rates
+
+    result1, result1_extra = odesys_extra['unit_aware_solve'](
+        tend_units, c0_units, p_units, integrator='cvode')
+    assert result1.info['success']
+
+    result2 = odesys.integrate(tend, c0, p, integrator='cvode')
+    assert np.allclose(result2.yout[-1, :], to_unitless(result1.yout[-1, :], u.molar))

chempy/util/parsing.py

@@ -326,13 +326,12 @@ def to_reaction(line, substance_keys, token, Cls, globals_=None, **kwargs):
     with running this on untrusted data.
 
     """
-    # TODO: add handling of units.
     if globals_ is None:
         import chempy
         from chempy.kinetics import rates
         from chempy.units import default_units
         globals_ = {k: getattr(rates, k) for k in dir(rates)}
-        globals_.update({'chempy': chempy, 'default_units': default_units})
+        globals_.update({'chempy': chempy})
         if default_units is not None:
             globals_.update({k: v for k, v in chempy.__dict__.items()
                              if not k.startswith('_')})

conda-recipe/meta.yaml

@@ -1,5 +1,5 @@
 {% set name = "chempy" %}
-{% set version = "0.6.0.dev0+git" %}
+{% set version = "0.6.8.dev0+git" %}
 
 package:
   name: {{ name|lower }}
@@ -19,16 +19,16 @@ requirements:
     - setuptools
   run:
     - python
-    - numpy >1.7
-    - scipy >=0.16.1
-    - matplotlib >=1.3.1
-    - sympy >=1.1.1
+    - numpy >=1.11.3
+    - scipy >=1.0.1
+    - matplotlib >=2.2.3
+    - sympy >=1.1.1,!=1.2
     - quantities >=0.12.1
     - pulp >=1.6.8
-    - pyneqsys >=0.5.1
-    - pyodesys >=0.11.16
+    - pyneqsys >=0.5.4
+    - pyodesys >=0.11.16,<0.12.0
     - pyparsing >=2.0.3
-    - sym >=0.3.0
+    - sym >=0.3.3
     - notebook
     - nbconvert
 

examples/_surface_reactions_homogeneous.ipynb

@@ -0,0 +1,140 @@
+{
+ "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.7.0+"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

setup.py

@@ -118,7 +118,7 @@ def _path_under_setup(*args):
     packages=[pkg_name] + submodules + tests,
     classifiers=classifiers,
     install_requires=[
-        'numpy>1.11.3', 'scipy>=1.0.1', 'matplotlib>=2.2.2',
+        'numpy>1.11.3', 'scipy>=1.0.1', 'matplotlib>=2.2.3',
         'sympy>=1.1.1,!=1.2', 'quantities>=0.12.1', 'pyneqsys>=0.5.4',
         'pyodesys>=0.11.16,!=0.12.0', 'pyparsing>=2.0.3', 'sym>=0.3.3', 'jupyter',
         'pulp>=1.6.8',