In [ ]:
%load_ext autoreload
%autoreload 2
In [ ]:
from collections import defaultdict
from IPython.display import Image
import numpy as np
import matplotlib.pyplot as plt
from pyodesys.symbolic import ScaledSys
from chempy import ReactionSystem
from chempy.kinetics.ode import get_odesys
from chempy.kinetics.analysis import plot_reaction_contributions, dominant_reactions_graph
from chempy.kinetics._native import get_native
%matplotlib inline
In [ ]:
import os
os.environ['OMP_NUM_THREADS'] = '1'
print(os.environ['OMP_NUM_THREADS'])
In [ ]:
# Rate constants choosen rather arbitrarily for demonstration purposes
system_str = """
2 Fe+2 + H2O2 -> 2 Fe+3 + 2 OH- ; 42
2 Fe+3 + H2O2 -> 2 Fe+2 + O2 + 2 H+ ; 17
H+ + OH- -> H2O ; %(kprot_H2O)s
H2O -> H+ + OH- ; %(kprot_H2O)s * %(Kw)s
Fe+3 + 2 H2O -> FeOOH(s) + 3 H+ ; 1
FeOOH(s) + 3 H+ -> Fe+3 + 2 H2O ; 2.5
H2O2 -> H+ + HO2- ; %(kprot_H2O2)s * %(Ka_H2O2)s
H+ + HO2- -> H2O2 ; %(kprot_H2O2)s
H2O2 + OH- -> HO2- + H2O ; %(kdepr_H2O2)s
HO2- + H2O -> H2O2 + OH- ; %(kdepr_H2O2)s / %(Ka_H2O2)s * %(Kw)s
""" % dict(Kw='1e-14', kprot_H2O='1e10', Ka_H2O2='10**-11.65', kprot_H2O2='5e10', kdepr_H2O2='1.3e10')
In [ ]:
rsys = ReactionSystem.from_string(system_str, name='Iron(II)/Iron(III) speciation') # "[H2O]" = 1.0 (actually 55.4 at RT)
odesys, extra = get_odesys(rsys, SymbolicSys=ScaledSys, dep_scaling=1e6)
rsys
In [ ]:
tout = sorted(np.concatenate((np.linspace(0, 1003), np.logspace(-8, 3))))
c0 = defaultdict(float, {'Fe+2': 0.05, 'H2O2': 0.1, 'H2O': 1.0, 'H+': 1e-7, 'OH-': 1e-7})
res = odesys.integrate(tout, c0, atol=1e-12, rtol=1e-14, integrator='gsl')
In [ ]:
def plot_result(result):
fig, axes = plt.subplots(1, 4, figsize=(16, 6))
result.plot(names=[k for k in rsys.substances if k != 'H2O'], ax=axes[0])
axes[0].legend(loc='best', prop={'size': 10}); _ = axes[0].set_xlabel('Time'); _ = axes[0].set_ylabel('Concentration')
result.plot(names=[k for k in rsys.substances if k != 'H2O'], ax=axes[1], xscale='log', yscale='log')
axes[1].legend(loc='best', prop={'size': 10}); _ = axes[1].set_xlabel('Time'); _ = axes[1].set_ylabel('Concentration')
result.plot(deriv=True, names=[k for k in rsys.substances if k != 'H2O'], ax=axes[2], xscale='log');
axes[2].set_yscale('symlog', linthreshy=1e-9)
axes[2].legend(loc='best', prop={'size': 10}); _ = axes[2].set_xlabel('Time'); _ = axes[2].set_ylabel('dC/dt')
result.plot_invariant_violations(ax=axes[3], xscale='log')
axes[3].legend(loc='best', prop={'size': 10}); _ = axes[3].set_xlabel('Time'); _ = axes[3].set_ylabel('Component conservation')
plt.tight_layout()
In [ ]:
plot_result(res)
In [ ]:
substance_keys = 'Fe+2 Fe+3 FeOOH(s) H2O2'.split()
selection = res.xout > 10
fig, axes = plt.subplots(len(substance_keys), 2, figsize=(16, 4*len(substance_keys)))
plot_reaction_contributions(res, rsys, extra['rate_exprs_cb'], substance_keys, axes=axes[:, 0])
plot_reaction_contributions(res, rsys, extra['rate_exprs_cb'],
substance_keys, axes=axes[:, 1], relative=True, xscale='linear', yscale='linear',
combine_equilibria=True, selection=selection)
plt.tight_layout()
In [ ]:
nativesys = get_native(rsys, odesys, 'cvode')
In [ ]:
nativeres = nativesys.integrate(tout[-1], c0, atol=1e-10, rtol=1e-10, nsteps=25000,
return_on_root=True, first_step=1e-40) #error_outside_bounds=True)
plot_result(nativeres)
In [ ]:
from _ipynb_progressbar import log_progress
native = get_native(rsys, odesys, 'cvode', steady_state_root=True)
def plot_tss_conv(factor, tols, ax, idx):
success = []
failure = []
c = 'krbgcmy'
for tol in log_progress(tols):
res = native.integrate(
1e11, c0, atol=tol, rtol=tol, nsteps=25000,
return_on_root=True, dx_min=1e-17, return_on_error=True, first_step=1e-40, special_settings=[factor])
if res.info['success']:
success.append((tol, res.xout[-1]))
else:
failure.append((tol, res.xout[-1]))
if success:
ax.loglog(*zip(*success), marker='.', color=c[idx % len(c)], label=factor)
if failure:
ax.loglog(*zip(*failure), marker='x', color=c[idx % len(c)], label=None if success else fac01tor)
In [ ]:
#fig, ax = plt.subplots(figsize=(14, 6))
#tols = np.logspace(-10, -6, 50)
#for idx, factor in enumerate([1e4, 1.1e4, 1e5, 1e6, 1e7, 1e8, 1e9]):
# plot_tss_conv(factor, tols, ax, idx)
#ax.legend()
In [ ]:
dominant_reactions_graph(res.yout[-1, :], extra['rate_exprs_cb'], rsys, 'H2O2', combine_equilibria=True,
fname='H2O2.png', output_dir='.', include_inactive=False, tex=False)
Image('H2O2.png')
In [ ]: