In [ ]:
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()}
In [ ]:
# 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 [ ]:
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 [ ]:
# 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 [ ]:
# The productions reactions have hardcoded rates corresponding to 1 Gy/s
rsys = ReactionSystem(prod_rxns + reactions, species)
rsys
In [ ]:
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 [ ]:
%debug
In [ ]:
uni
In [ ]:
bi
In [ ]:
odesys, extra = get_odesys(rsys)
odesys.exprs[:3] # take a look at the first three ODEs in the system
In [ ]:
j = odesys.get_jac()
j.shape
In [ ]:
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 [ ]:
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 [ ]:
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 [ ]:
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 [ ]:
integrate_and_plot(odesys, c0_dict, 'cvode', first_step=1e-12)
In [ ]:
integrate_and_plot(odesys, c0_dict, 'scipy')
One way to avoid negative concentrations is to solve for the logarithm of the concentration instead of the concentration directly:
In [ ]:
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 [ ]:
loglogsys.exprs[:3]
In [ ]:
integrate_and_plot(loglogsys, c0_dict, 'gsl', zero_conc=1e-26, first_step=1e-3, log10_t0=-13)
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 [ ]:
ssys, _ = get_odesys(rsys, SymbolicSys=ScaledSys, dep_scaling=1e16)
for speed we will use native compiled C++ code for the numerical evalatuations:
In [ ]:
from chempy.kinetics._native import get_native
nsys = get_native(rsys, ssys, 'cvode')
In [ ]:
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)
In [ ]:
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)
We can also access the generated C++ code (which can be handy if it needs to be run where Python is not available)
In [ ]:
print(''.join(open(next(filter(lambda s: s.endswith('.cpp'), nsys._native._written_files))).readlines()[:42]))
In [ ]:
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 [ ]:
from chempy.kinetics.analysis import plot_reaction_contributions
In [ ]:
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()