Aqueous radiolysis

In this notebook we will look at the reactive species in water subjected to ionizing radiation. The reaction-set is simply has been taken from the litterature and we will not pay to much attention to it, but rather focus on the analysis of dominant reactions toward the end.

In [1]:
from __future__ import absolute_import, division, print_function
from collections import defaultdict
import numpy as np
import sympy
import matplotlib.pyplot as plt
import chempy
from chempy import Reaction, Substance, ReactionSystem
from chempy.kinetics.ode import get_odesys
from chempy.kinetics.analysis import plot_reaction_contributions
from chempy.printing.tables import UnimolecularTable, BimolecularTable
import pyodesys
from pyodesys.symbolic import ScaledSys
from pyneqsys.plotting import mpl_outside_legend
sympy.init_printing()
%matplotlib inline
{k: globals()[k].__version__ for k in 'sympy pyodesys chempy'.split()}
Out[1]:
{'sympy': '1.4', 'pyodesys': '0.12.5', 'chempy': '0.8.0.dev0+git'}
In [2]:
# Never mind the next row, it contains all the reaction data of aqueous radiolysis at room temperature
_species, _reactions = ['H', 'H+', 'H2', 'e-(aq)', 'HO2-', 'HO2', 'H2O2', 'HO3', 'O-', 'O2', 'O2-', 'O3', 'O3-', 'OH', 'OH-', 'H2O'], [({1: 1, 14: 1}, {15: 1}, 140000000000.0, {}), ({15: 1}, {1: 1, 14: 1}, 2.532901323662559e-05, {}), ({6: 1}, {1: 1, 4: 1}, 0.11193605692841689, {}), ({1: 1, 4: 1}, {6: 1}, 50000000000.0, {}), ({6: 1, 14: 1}, {4: 1, 15: 1}, 13000000000.0, {}), ({4: 1, 15: 0}, {6: 1, 14: 1}, 58202729.542927094, {15: 1}), ({3: 1, 15: 0}, {0: 1, 14: 1}, 19.0, {15: 1}), ({0: 1, 14: 1}, {3: 1, 15: 1}, 18000000.0, {}), ({0: 1}, {1: 1, 3: 1}, 3.905960400662016, {}), ({1: 1, 3: 1}, {0: 1}, 23000000000.0, {}), ({13: 1, 14: 1}, {8: 1, 15: 1}, 13000000000.0, {}), ({8: 1, 15: 0}, {13: 1, 14: 1}, 103500715.55425139, {15: 1}), ({13: 1}, {8: 1, 1: 1}, 0.12589254117941662, {}), ({8: 1, 1: 1}, {13: 1}, 100000000000.0, {}), ({5: 1}, {1: 1, 10: 1}, 1345767.401963457, {}), ({1: 1, 10: 1}, {5: 1}, 50000000000.0, {}), ({5: 1, 14: 1}, {10: 1, 15: 1}, 50000000000.0, {}), ({10: 1, 15: 0}, {5: 1, 14: 1}, 18.619585312728415, {15: 1}), ({3: 1, 13: 1}, {14: 1}, 30000000000.0, {}), ({3: 1, 6: 1}, {13: 1, 14: 1}, 14000000000.0, {}), ({10: 1, 3: 1, 15: 0}, {4: 1, 14: 1}, 13000000000.0, {15: 1}), ({3: 1, 5: 1}, {4: 1}, 20000000000.0, {}), ({9: 1, 3: 1}, {10: 1}, 22200000000.0, {}), ({3: 2, 15: 0}, {2: 1, 14: 2}, 5000000000.0, {15: 2}), ({0: 1, 3: 1, 15: 0}, {2: 1, 14: 1}, 25000000000.0, {15: 1}), ({3: 1, 4: 1}, {8: 1, 14: 1}, 3500000000.0, {}), ({8: 1, 3: 1, 15: 0}, {14: 2}, 22000000000.0, {15: 1}), ({3: 1, 12: 1, 15: 0}, {9: 1, 14: 2}, 16000000000.0, {15: 1}), ({3: 1, 11: 1}, {12: 1}, 36000000000.0, {}), ({0: 1, 15: 0}, {2: 1, 13: 1}, 11.0, {15: 1}), ({0: 1, 8: 1}, {14: 1}, 10000000000.0, {}), ({0: 1, 4: 1}, {13: 1, 14: 1}, 90000000.0, {}), ({0: 1, 12: 1}, {9: 1, 14: 1}, 10000000000.0, {}), ({0: 2}, {2: 1}, 7750000000.0, {}), ({0: 1, 13: 1}, {15: 1}, 7000000000.0, {}), ({0: 1, 6: 1}, {13: 1, 15: 1}, 90000000.0, {}), ({0: 1, 9: 1}, {5: 1}, 21000000000.0, {}), ({0: 1, 5: 1}, {6: 1}, 10000000000.0, {}), ({0: 1, 10: 1}, {4: 1}, 20000000000.0, {}), ({0: 1, 11: 1}, {7: 1}, 38000000000.0, {}), ({13: 2}, {6: 1}, 3600000000.0, {}), ({5: 1, 13: 1}, {9: 1, 15: 1}, 6000000000.0, {}), ({10: 1, 13: 1}, {9: 1, 14: 1}, 8200000000.0, {}), ({2: 1, 13: 1}, {0: 1, 15: 1}, 40000000.0, {}), ({13: 1, 6: 1}, {5: 1, 15: 1}, 30000000.0, {}), ({8: 1, 13: 1}, {4: 1}, 20000000000.0, {}), ({4: 1, 13: 1}, {5: 1, 14: 1}, 7500000000.0, {}), ({12: 1, 13: 1}, {11: 1, 14: 1}, 2550000000.0, {}), ({12: 1, 13: 1}, {1: 1, 10: 2}, 5950000000.0, {}), ({11: 1, 13: 1}, {9: 1, 5: 1}, 110000000.0, {}), ({10: 1, 5: 1}, {9: 1, 4: 1}, 80000000.0, {}), ({5: 2}, {9: 1, 6: 1}, 800000.0, {}), ({8: 1, 5: 1}, {9: 1, 14: 1}, 6000000000.0, {}), ({5: 1, 6: 1}, {9: 1, 13: 1, 15: 1}, 0.5, {}), ({4: 1, 5: 1}, {9: 1, 13: 1, 14: 1}, 0.5, {}), ({12: 1, 5: 1}, {9: 2, 14: 1}, 6000000000.0, {}), ({11: 1, 5: 1}, {9: 1, 7: 1}, 500000000.0, {}), ({10: 2, 15: 0}, {9: 1, 6: 1, 14: 2}, 100.0, {15: 2}), ({8: 1, 10: 1, 15: 0}, {9: 1, 14: 2}, 600000000.0, {15: 1}), ({10: 1, 6: 1}, {9: 1, 13: 1, 14: 1}, 0.13, {}), ({10: 1, 4: 1}, {8: 1, 9: 1, 14: 1}, 0.13, {}), ({10: 1, 12: 1, 15: 0}, {9: 2, 14: 2}, 10000.0, {15: 1}), ({10: 1, 11: 1}, {9: 1, 12: 1}, 1500000000.0, {}), ({8: 2, 15: 0}, {4: 1, 14: 1}, 1000000000.0, {15: 1}), ({8: 1, 9: 1}, {12: 1}, 3600000000.0, {}), ({8: 1, 2: 1}, {0: 1, 14: 1}, 80000000.0, {}), ({8: 1, 6: 1}, {10: 1, 15: 1}, 500000000.0, {}), ({8: 1, 4: 1}, {10: 1, 14: 1}, 400000000.0, {}), ({8: 1, 12: 1}, {10: 2}, 700000000.0, {}), ({8: 1, 11: 1}, {9: 1, 10: 1}, 5000000000.0, {}), ({12: 1}, {8: 1, 9: 1}, 300.0, {}), ({1: 1, 12: 1}, {9: 1, 13: 1}, 90000000000.0, {}), ({12: 1, 6: 1}, {9: 1, 10: 1, 15: 1}, 1600000.0, {}), ({12: 1, 4: 1}, {9: 1, 10: 1, 14: 1}, 890000.0, {}), ({2: 1, 12: 1}, {0: 1, 9: 1, 14: 1}, 250000.0, {}), ({7: 1}, {9: 1, 13: 1}, 110000.0, {}), ({13: 1, 7: 1}, {9: 1, 6: 1}, 5000000000.0, {}), ({7: 2}, {9: 2, 6: 1}, 5000000000.0, {}), ({10: 1, 7: 1}, {9: 2, 14: 1}, 10000000000.0, {}), ({7: 1}, {1: 1, 12: 1}, 328.097819129701, {}), ({1: 1, 12: 1}, {7: 1}, 52000000000.0, {}), ({11: 1, 14: 1}, {10: 1, 5: 1}, 70.0, {}), ({11: 1, 4: 1}, {9: 1, 10: 1, 13: 1}, 2800000.0, {})]
In [3]:
species = [Substance.from_formula(s) for s in _species]
reactions = [
    Reaction({_species[k]: v for k, v in reac.items()},
             {_species[k]: v for k, v in prod.items()}, param,
             {_species[k]: v for k, v in inact_reac.items()})
    for reac, prod, param, inact_reac in _reactions
]
In [4]:
# radiolytic yields for gamma radiolysis of neat water
C_H2O = 55.5
YIELD_CONV = 1.0364e-07 # mol * eV / (J * molecules)
prod_rxns = [
    Reaction({'H2O': 1}, {'H+': 1,  'OH-': 1},            0.5  * YIELD_CONV / C_H2O),
    Reaction({'H2O': 1}, {'H+': 1, 'e-(aq)': 1, 'OH': 1}, 2.6  * YIELD_CONV / C_H2O),
    Reaction({'H2O': 1}, {'H':  2, 'H2O2': 1},            0.66 * YIELD_CONV / C_H2O, {'H2O': 1}),
    Reaction({'H2O': 1}, {'H2': 1, 'H2O2': 1},            0.74 * YIELD_CONV / C_H2O, {'H2O': 1}),
    Reaction({'H2O': 1}, {'H2': 1,   'OH': 2},            0.1  * YIELD_CONV / C_H2O, {'H2O': 1}),
    Reaction({'H2O': 1}, {'H2': 3,  'HO2': 2},            0.04 * YIELD_CONV / C_H2O, {'H2O': 3}),
]
In [5]:
# The productions reactions have hardcoded rates corresponding to 1 Gy/s
rsys = ReactionSystem(prod_rxns + reactions, species)
rsys
Out[5]:
H2OH+ + OH- 9.3369⋅10-10
H2OH+ + OH + e-(aq) 4.8552⋅10-9
H2O + ( H2O)2 H + H2O2 1.2325⋅10-9
H2O + ( H2O)H2 + H2O2 1.3819⋅10-9
H2O + ( H2O)H2 + 2 OH 1.8674⋅10-10
H2O + ( 3 H2O)3 H2 + 2 HO2 7.4695⋅10-11
H+ + OH-H2O 1.4⋅1011
H2OH+ + OH- 2.5329⋅10-5
H2O2H+ + HO2- 0.11194
H+ + HO2-H2O2 5⋅1010
H2O2 + OH-H2O + HO2- 1.3⋅1010
HO2- + ( H2O)H2O2 + OH- 5.8203⋅107
e-(aq) + ( H2O)H + OH- 19
H + OH-H2O + e-(aq) 1.8⋅107
HH+ + e-(aq) 3.906
H+ + e-(aq)H 2.3⋅1010
OH + OH-H2O + O- 1.3⋅1010
O- + ( H2O)OH + OH- 1.035⋅108
OHH+ + O- 0.12589
H+ + O-OH 1011
HO2H+ + O2- 1.3458⋅106
H+ + O2-HO2 5⋅1010
HO2 + OH-H2O + O2- 5⋅1010
O2- + ( H2O)HO2 + OH- 18.62
OH + e-(aq)OH- 3⋅1010
H2O2 + e-(aq)OH + OH- 1.4⋅1010
O2- + e-(aq) + ( H2O)HO2- + OH- 1.3⋅1010
HO2 + e-(aq)HO2- 2⋅1010
O2 + e-(aq)O2- 2.22⋅1010
2 e-(aq) + ( 2 H2O)H2 + 2 OH- 5⋅109
H + e-(aq) + ( H2O)H2 + OH- 2.5⋅1010
HO2- + e-(aq)O- + OH- 3.5⋅109
O- + e-(aq) + ( H2O)2 OH- 2.2⋅1010
O3- + e-(aq) + ( H2O)O2 + 2 OH- 1.6⋅1010
O3 + e-(aq)O3- 3.6⋅1010
H + ( H2O)H2 + OH 11
H + O-OH- 1010
H + HO2-OH + OH- 9⋅107
H + O3-O2 + OH- 1010
2 HH2 7.75⋅109
H + OHH2O 7⋅109
H + H2O2H2O + OH 9⋅107
H + O2HO2 2.1⋅1010
H + HO2H2O2 1010
H + O2-HO2- 2⋅1010
H + O3HO3 3.8⋅1010
2 OHH2O2 3.6⋅109
HO2 + OHH2O + O2 6⋅109
O2- + OHO2 + OH- 8.2⋅109
H2 + OHH + H2O 4⋅107
H2O2 + OHH2O + HO2 3⋅107
O- + OHHO2- 2⋅1010
HO2- + OHHO2 + OH- 7.5⋅109
O3- + OHO3 + OH- 2.55⋅109
O3- + OHH+ + 2 O2- 5.95⋅109
O3 + OHHO2 + O2 1.1⋅108
HO2 + O2-HO2- + O2 8⋅107
2 HO2H2O2 + O2 8⋅105
HO2 + O-O2 + OH- 6⋅109
H2O2 + HO2H2O + O2 + OH 0.5
HO2 + HO2-O2 + OH + OH- 0.5
HO2 + O3-2 O2 + OH- 6⋅109
HO2 + O3HO3 + O2 5⋅108
2 O2- + ( 2 H2O)H2O2 + O2 + 2 OH- 100
O- + O2- + ( H2O)O2 + 2 OH- 6⋅108
H2O2 + O2-O2 + OH + OH- 0.13
HO2- + O2-O- + O2 + OH- 0.13
O2- + O3- + ( H2O)2 O2 + 2 OH- 10000
O2- + O3O2 + O3- 1.5⋅109
2 O- + ( H2O)HO2- + OH- 109
O- + O2O3- 3.6⋅109
H2 + O-H + OH- 8⋅107
H2O2 + O-H2O + O2- 5⋅108
HO2- + O-O2- + OH- 4⋅108
O- + O3-2 O2- 7⋅108
O- + O3O2 + O2- 5⋅109
O3-O- + O2 300
H+ + O3-O2 + OH 9⋅1010
H2O2 + O3-H2O + O2 + O2- 1.6⋅106
HO2- + O3-O2 + O2- + OH- 8.9⋅105
H2 + O3-H + O2 + OH- 2.5⋅105
HO3O2 + OH 1.1⋅105
HO3 + OHH2O2 + O2 5⋅109
2 HO3H2O2 + 2 O2 5⋅109
HO3 + O2-2 O2 + OH- 1010
HO3H+ + O3- 328.1
H+ + O3-HO3 5.2⋅1010
O3 + OH-HO2 + O2- 70
HO2- + O3O2 + O2- + OH 2.8⋅106
In [6]:
uni, not_uni = UnimolecularTable.from_ReactionSystem(rsys)
bi, not_bi = BimolecularTable.from_ReactionSystem(rsys)
assert not (not_uni & not_bi), "There are only uni- & bi-molecular reactions in this set"
In [7]:
%debug
ERROR:root:No traceback has been produced, nothing to debug.
In [10]:
odesys, extra = get_odesys(rsys)
odesys.exprs[:3]  # take a look at the first three ODEs in the system
Out[10]:
$\displaystyle \left( - 15500000000.0 y_{0}^{2} - 20000000000.0 y_{0} y_{10} - 38000000000.0 y_{0} y_{11} - 10000000000.0 y_{0} y_{12} - 7000000000.0 y_{0} y_{13} - 18000000.0 y_{0} y_{14} - 25000000000.0 y_{0} y_{3} - 90000000.0 y_{0} y_{4} - 10000000000.0 y_{0} y_{5} - 90000000.0 y_{0} y_{6} - 10000000000.0 y_{0} y_{8} - 21000000000.0 y_{0} y_{9} - 14.905960400662 y_{0} + 23000000000.0 y_{1} y_{3} + 250000.0 y_{12} y_{2} + 40000000.0 y_{13} y_{2} + 2.46495135135135 \cdot 10^{-9} y_{15} + 80000000.0 y_{2} y_{8} + 19.0 y_{3}, \ 3.90596040066202 y_{0} - 50000000000.0 y_{1} y_{10} - 142000000000.0 y_{1} y_{12} - 140000000000.0 y_{1} y_{14} - 23000000000.0 y_{1} y_{3} - 50000000000.0 y_{1} y_{4} - 100000000000.0 y_{1} y_{8} + 5950000000.0 y_{12} y_{13} + 0.125892541179417 y_{13} + 2.53348021375265 \cdot 10^{-5} y_{15} + 1345767.40196346 y_{5} + 0.111936056928417 y_{6} + 328.097819129701 y_{7}, \ 7750000000.0 y_{0}^{2} + 25000000000.0 y_{0} y_{3} + 11.0 y_{0} - 250000.0 y_{12} y_{2} - 40000000.0 y_{13} y_{2} + 1.79269189189189 \cdot 10^{-9} y_{15} - 80000000.0 y_{2} y_{8} + 5000000000.0 y_{3}^{2}\right)$
In [11]:
j = odesys.get_jac()
j.shape
Out[11]:
$\displaystyle \left( 16, \ 16\right)$
In [12]:
def integrate_and_plot(odesys, c0_dict, integrator, ax=None, zero_conc=0, log10_t0=-12, tout=None,
                       print_info=False, **kwargs):
    if tout is None:
        tout = np.logspace(log10_t0, 6)
    c0 = [c0_dict.get(k, zero_conc)+zero_conc for k in _species]
    result = odesys.integrate(tout, c0, integrator=integrator, **kwargs)
    if ax is None:
        ax = plt.subplot(1, 1, 1)
    result.plot(ax=ax, title_info=2)
    ax.set_xscale('log'); ax.set_yscale('log')
    mpl_outside_legend(ax, prop={'size': 9})
    if print_info:
        print({k: v for k, v in result.info.items() if not k.startswith('internal')})
    return result
In [13]:
c0_dict = defaultdict(float, {'H+': 1e-7, 'OH-': 1e-7, 'H2O': 55.5})
plt.figure(figsize=(16,6))
integrate_and_plot(odesys, c0_dict, integrator='cvode', first_step=1e-14, atol=1e-10, rtol=1e-10, autorestart=5)
plt.legend()
_ = plt.ylim([1e-16, 60])
In [14]:
def integrate_and_plot_scaled(rsys, dep_scaling, *args, **kwargs):
    integrate_and_plot(get_odesys(
            rsys, SymbolicSys=ScaledSys, dep_scaling=dep_scaling)[0], *args, **kwargs)

Let's see how different solvers behave when integrating this problem:

In [15]:
fig, axes = plt.subplots(2, 2, figsize=(14, 14))
solvers = 'cvode', 'odeint', 'gsl', 'scipy'
for idx, ax in enumerate(np.ravel(axes)):
    #if idx == 1:
    #    continue
    kw = {'method': 'vode'} if solvers[idx] == 'scipy' else {}
    integrate_and_plot_scaled(rsys, 1e3, c0_dict, solvers[idx], ax=ax, atol=1e-6, rtol=1e-6,
                              nsteps=4000, first_step=1e-10*extra['max_euler_step_cb'](0, c0_dict), **kw)
    ax.set_title(solvers[idx])
In [16]:
integrate_and_plot(odesys, c0_dict, 'cvode', first_step=1e-12)
Out[16]:
<pyodesys.results.Result at 0x7f28ebbf0dd8>
In [17]:
integrate_and_plot(odesys, c0_dict, 'scipy')
Out[17]:
<pyodesys.results.Result at 0x7f28eb667ef0>

One way to avoid negative concentrations is to solve for the logarithm of the concentration instead of the concentration directly:

In [18]:
from pyodesys.symbolic import symmetricsys
logexp = (sympy.log, sympy.exp)
LogLogSys = symmetricsys(logexp, logexp, exprs_process_cb=lambda exprs: [
        sympy.powsimp(expr.expand(), force=True) for expr in exprs])
loglogsys, _ = get_odesys(rsys, SymbolicSys=LogLogSys)
In [19]:
loglogsys.exprs[:3]
Out[19]:
$\displaystyle \left( - 14.905960400662 e^{x_{0}} - 15500000000.0 e^{x_{0} + y_{0}} - 20000000000.0 e^{x_{0} + y_{10}} - 38000000000.0 e^{x_{0} + y_{11}} - 10000000000.0 e^{x_{0} + y_{12}} - 7000000000.0 e^{x_{0} + y_{13}} - 18000000.0 e^{x_{0} + y_{14}} - 25000000000.0 e^{x_{0} + y_{3}} - 90000000.0 e^{x_{0} + y_{4}} - 10000000000.0 e^{x_{0} + y_{5}} - 90000000.0 e^{x_{0} + y_{6}} - 10000000000.0 e^{x_{0} + y_{8}} - 21000000000.0 e^{x_{0} + y_{9}} + 2.46495135135135 \cdot 10^{-9} e^{x_{0} - y_{0} + y_{15}} + 19.0 e^{x_{0} - y_{0} + y_{3}} + 23000000000.0 e^{x_{0} - y_{0} + y_{1} + y_{3}} + 250000.0 e^{x_{0} - y_{0} + y_{12} + y_{2}} + 40000000.0 e^{x_{0} - y_{0} + y_{13} + y_{2}} + 80000000.0 e^{x_{0} - y_{0} + y_{2} + y_{8}}, \ - 50000000000.0 e^{x_{0} + y_{10}} - 142000000000.0 e^{x_{0} + y_{12}} - 140000000000.0 e^{x_{0} + y_{14}} - 23000000000.0 e^{x_{0} + y_{3}} - 50000000000.0 e^{x_{0} + y_{4}} - 100000000000.0 e^{x_{0} + y_{8}} + 3.90596040066202 e^{x_{0} + y_{0} - y_{1}} + 0.125892541179417 e^{x_{0} - y_{1} + y_{13}} + 2.53348021375265 \cdot 10^{-5} e^{x_{0} - y_{1} + y_{15}} + 1345767.40196346 e^{x_{0} - y_{1} + y_{5}} + 0.111936056928417 e^{x_{0} - y_{1} + y_{6}} + 328.097819129701 e^{x_{0} - y_{1} + y_{7}} + 5950000000.0 e^{x_{0} - y_{1} + y_{12} + y_{13}}, \ - 250000.0 e^{x_{0} + y_{12}} - 40000000.0 e^{x_{0} + y_{13}} - 80000000.0 e^{x_{0} + y_{8}} + 11.0 e^{x_{0} + y_{0} - y_{2}} + 7750000000.0 e^{x_{0} + 2 y_{0} - y_{2}} + 1.79269189189189 \cdot 10^{-9} e^{x_{0} + y_{15} - y_{2}} + 5000000000.0 e^{x_{0} - y_{2} + 2 y_{3}} + 25000000000.0 e^{x_{0} + y_{0} - y_{2} + y_{3}}\right)$
In [20]:
integrate_and_plot(loglogsys, c0_dict, 'gsl', zero_conc=1e-26, first_step=1e-3, log10_t0=-13)
/usr/local/lib/python3.6/dist-packages/numpy/__init__.py:1: RuntimeWarning: overflow encountered in exp
  """
/usr/local/lib/python3.6/dist-packages/numpy/__init__.py:1: RuntimeWarning: invalid value encountered in double_scalars
  """
Out[20]:
<pyodesys.results.Result at 0x7f28eb7eb668>

It works, but at the cost of signigicant overhead (much larger number of function evaluations were needed). Another option is to scale the problem:

In [21]:
ssys, _ = get_odesys(rsys, SymbolicSys=ScaledSys, dep_scaling=1e16)

for speed we will use native compiled C++ code for the numerical evalatuations:

In [22]:
from chempy.kinetics._native import get_native
nsys = get_native(rsys, ssys, 'cvode')
In [23]:
integrate_and_plot(nsys, c0_dict, 'cvode', nsteps=96000,
                   error_outside_bounds=True, get_dx_max_factor=-1.0, autorestart=2, atol=1e-9, rtol=1e-10)
Out[23]:
<pyodesys.results.Result at 0x7f28e7f5dfd0>
In [24]:
integrate_and_plot(nsys, c0_dict, 'cvode', nsteps=96000, tout=(1e-16, 1e6),
                   error_outside_bounds=True, get_dx_max_factor=-1.0, autorestart=2, atol=1e-9, rtol=1e-10,
                   print_info=True, first_step=1e-16)
{'n_err_test_fails': 247, 'n_lin_solv_setups': 2427, 'n_rhs_evals': 11647, 'n_nonlin_solv_conv_fails': 512, 'njev': 1047, 'n_steps': 6679, 'nfev': 11647, 'njvev': 0, 'n_dls_rhs_evals': 0, 'n_dls_jac_evals': 1047, 'n_root_evals': 0, 'n_nonlin_solv_iters': 11646, 'time_jac': 0.00098874, 'time_cpu': 0.095827, 'time_quads': 0.0, 'time_prec': 0.0, 'time_jtimes': 0.0, 'time_rhs': 0.007206880000000014, 'time_wall': 0.09586960200000001, 'time_roots': 0.0, 'steps': array([], dtype=float64), 'orders': array([], dtype=int64), 'fpes': array([], dtype=int64), 'root_indices': [], 'mode': 'adaptive', 'success': True}
Out[24]:
<pyodesys.results.Result at 0x7f28e7f89710>

We can also access the generated C++ code (which can be handy if it needs to be run where Python is not available)

In [25]:
print(''.join(open(next(filter(lambda s: s.endswith('.cpp'), nsys._native._written_files))).readlines()[:42]))
// -*- mode: read-only -*-
// This file was generated using pyodesys-0.12.5 at 2019-10-07T00:39:11.421784
// -*- coding: utf-8 -*-

// User provided system description: None
// Names of dependent variables: ('H', 'H+', 'H2', 'e-(aq)', 'HO2-', 'HO2', 'H2O2', 'HO3', 'O-', 'O2', 'O2-', 'O3', 'O3-', 'OH', 'OH-', 'H2O')
// Names of parameters: ()

#include <algorithm>
#include <iostream>
#include <limits>
#include <cmath>
#include <vector>
#include "odesys_anyode_iterative.hpp"
#include <vector>
#include <type_traits>

namespace {  // anonymous namespace for user-defined helper functions
    std::vector<std::string> p_odesys_names = {"H", "H+", "H2", "e-(aq)", "HO2-", "HO2", "H2O2", "HO3", "O-", "O2", "O2-", "O3", "O3-", "OH", "OH-", "H2O"};
    
    template <typename T>
    constexpr T vecmin(T&& a){
        return std::forward<T>(a);
    }
    template <typename T1, typename T2>
    constexpr typename std::common_type<T1, T2>::type vecmin(T1&& a, T2&& b){
        return (a < b) ? std::forward<T1>(a) : std::forward<T2>(b);
    }
    template <typename T, typename... Ts>
    constexpr typename std::common_type<T, Ts...>::type vecmin(T&& a, Ts&&... args){
        return vecmin(std::forward<T>(a), vecmin(std::forward<Ts>(args)...));
    }
    std::vector<double> upper_conc_bounds(const double * const y){
        auto bounds = std::vector<double>(16);
        double cc[2];
        cc[0] = y[0] + y[1] + 2*y[2] + y[4] + y[5] + 2*y[6] + y[7] + y[13] + y[14] + 2*y[15];
        cc[1] = 2*y[4] + 2*y[5] + 2*y[6] + 3*y[7] + y[8] + 2*y[9] + 2*y[10] + 3*y[11] + 3*y[12] + y[13] + y[14] + y[15];
        bounds[0] = vecmin(cc[0]);
        bounds[1] = vecmin(cc[0]);
        bounds[2] = vecmin(cc[0]/2);
        bounds[3] = INFINITY;
        bounds[4] = vecmin(cc[0], cc[1]/2);

In [26]:
init_conc_H2O2 = c0_dict.copy()
init_conc_H2O2['H2O2'] = 0.01
odesys2 = ScaledSys.from_other(odesys, lower_bounds=0, dep_scaling=1e8, indep_scaling=1e-6)
res = integrate_and_plot(odesys2, init_conc_H2O2, integrator='cvode', nsteps=500,
                         first_step=1e-16, atol=1e-5, rtol=1e-10, autorestart=1)

In order to plot the most important reactions vs. time we can use a function provided by ChemPy:

In [27]:
from chempy.kinetics.analysis import plot_reaction_contributions
In [28]:
sks = ['H2O2', 'OH', 'HO2']
fig, axes = plt.subplots(1, len(sks), figsize=(16, 10))
selection = res.xout > 1e0
plot_reaction_contributions(res, rsys, extra['rate_exprs_cb'], selection=selection,
                            substance_keys=sks, axes=axes, combine_equilibria=True, total=True)
plt.tight_layout()