In [ ]:

    

from collections import defaultdict
import numpy as np
import matplotlib.pyplot as plt
from pyodesys.symbolic import symmetricsys, get_logexp
from pyodesys.native import native_sys
from _chem_kinet import from_b64_gz_pkl, radiolysis1
%matplotlib inline


    

In [ ]:

    

linsys, (params,) = from_b64_gz_pkl(radiolysis1, steady_state_root=True)


    

In [ ]:

    

def integrate_and_plot(odesys, tend=1e6, H2O2=(0, 1e-4, 1e-3), **kwargs):
    res = odesys.integrate(tend, defaultdict(
            float, {'H2O': 55.4e3, 'H+': 1e-4, 'OH-': 1e-4, 'O2': 1e-6, 'H2O2': H2O2}
        ), params, integrator='cvode', **kwargs)
    fig, axes = plt.subplots(3, len(res), figsize=(16, 15))
    axes = np.reshape(axes, (3, len(res)))
    plot_kw = dict(legend=dict(loc='center left', prop={'size': 8}), xscale='log')
    for r, axcol in zip(res, axes.T):
        r.plot(ax=axcol[0], yscale='log', title_info=2, **plot_kw)
        r.plot(ax=axcol[1], deriv=True, yscale='symlog;linthreshy=1e-9', **plot_kw)
        r.plot_invariant_violations(ax=axcol[2], **plot_kw)
        for ax in axcol:
            ax.set_xlim([1e-13, 1e-3])
            ax.set_xlim([1e-6, tend])


    

In [ ]:

    

integrate_and_plot(linsys, nsteps=5000, atol=1e-13, rtol=1e-13)


    

In [ ]:

    

integrate_and_plot(linsys, return_on_root=True, nsteps=5000, atol=1e-13, rtol=1e-13)


    

In [ ]:

    

native = native_sys['cvode'].from_other(linsys)
print(next(filter(lambda s: s.endswith('.cpp'), native._native._written_files)))


    

In [ ]:

    

integrate_and_plot(native, nsteps=5000, atol=1e-13, rtol=1e-13, return_on_root=True)


    

In [ ]:

    

integrate_and_plot(native, nsteps=1024000, atol=1e-11, rtol=1e-11, tend=1e10, H2O2=(0,))


    

In [ ]:

    

try:
    integrate_and_plot(native, nsteps=1024000, atol=1e-11, rtol=1e-11, tend=1e10, H2O2=(0,),
                       max_invariant_violation=2e-8, autorestart=1)
except Exception as e:
    print("Raises:\n%s" % str(e))
else:
    print("No exception")


    

In [ ]:

    

try:
    integrate_and_plot(native, nsteps=1024000, atol=1e-11, rtol=1e-11, tend=1e10, H2O2=(0,),
                       max_invariant_violation=-5e-9, autorestart=0)
except Exception as e:
    print("Raises:\n%s" % str(e))
else:
    print("No exception")


    

In [ ]:

    

# would neeed to tile in pre_process:
#LogLogSys = symmetricsys(get_logexp(1, 1e-24), get_logexp(1, 1e-12))
#tsys = LogLogSys.from_other(linsys)
#integrate_and_plot(tsys)


    

In [ ]: