Kinetic modeling – sinusoidal temperature

In this notebook we will look at the solution of a kinetic problem where the temperature varies sinusoidally.

In [1]:
import math
from collections import defaultdict
from chempy import ReactionSystem
from chempy.kinetics.ode import get_odesys
from chempy.kinetics.rates import SinTemp
%matplotlib inline
In [2]:
rsys = ReactionSystem.from_string("""
2 HNO2 -> H2O + NO + NO2; EyringParam(dH=85e3, dS=10)
2 NO2 -> N2O4; EyringParam(dH=70e3, dS=20)
"""
)
In [3]:
st = SinTemp(unique_keys='Tbase Tamp Tangvel Tphase'.split())
In [4]:
odesys, extra = get_odesys(rsys, include_params=False, substitutions={'temperature': st})
In [5]:
init_conc = defaultdict(lambda: 0, HNO2=1, H2O=55)
params = dict(
    Tbase=300,
    Tamp=10,
    Tangvel=2*math.pi/10,
    Tphase=-math.pi/2
)
duration = 60
In [6]:
def integrate_and_plot(system):
    result = system.integrate(duration, init_conc, params, integrator='cvode', nsteps=2000)
    result.plot(names='NO HNO2 N2O4 NO2'.split(), title_info=2)
    print({k: v for k, v in sorted(result.info.items()) if not k.startswith('internal')})
In [7]:
integrate_and_plot(odesys)
{'atol': [1e-08], 'fpes': array([], dtype=int64), 'mode': 'adaptive', 'n_dls_jac_evals': 8, 'n_dls_rhs_evals': 0, 'n_err_test_fails': 23, 'n_lin_solv_setups': 57, 'n_nonlin_solv_conv_fails': 0, 'n_nonlin_solv_iters': 504, 'n_rhs_evals': 508, 'n_root_evals': 0, 'n_steps': 394, 'nfev': 508, 'njev': 8, 'njvev': 0, 'orders': array([], dtype=int64), 'rtol': 1e-08, 'steps': array([], dtype=float64), 'success': True, 'time_cpu': 0.186671, 'time_jac': 0.002946671, 'time_jtimes': 0.0, 'time_prec': 0.0, 'time_quads': 0.0, 'time_rhs': 0.17997687699999995, 'time_roots': 0.0, 'time_wall': 0.186610605}

Since we are using pyodesys we can reformulate our ODE-system as an autonomous one. This can help for stiff systems since time will then be explicitly represented in the Jacobian matrix.

In [8]:
autsys = odesys.as_autonomous()
assert len(autsys.dep) == len(odesys.dep) + 1
In [9]:
autsys.exprs
Out[9]:
(1.0*y_1**2*(69368672484.0789*p_0 + 69368672484.0789*p_1*sin(p_2*x + p_3))*exp(-10223.1386430792/(p_0 + p_1*sin(p_2*x + p_3))),
 -2.0*y_1**2*(69368672484.0789*p_0 + 69368672484.0789*p_1*sin(p_2*x + p_3))*exp(-10223.1386430792/(p_0 + p_1*sin(p_2*x + p_3))),
 1.0*y_4**2*(230939911607.726*p_0 + 230939911607.726*p_1*sin(p_2*x + p_3))*exp(-8419.05535312405/(p_0 + p_1*sin(p_2*x + p_3))),
 1.0*y_1**2*(69368672484.0789*p_0 + 69368672484.0789*p_1*sin(p_2*x + p_3))*exp(-10223.1386430792/(p_0 + p_1*sin(p_2*x + p_3))),
 1.0*y_1**2*(69368672484.0789*p_0 + 69368672484.0789*p_1*sin(p_2*x + p_3))*exp(-10223.1386430792/(p_0 + p_1*sin(p_2*x + p_3))) - 2.0*y_4**2*(230939911607.726*p_0 + 230939911607.726*p_1*sin(p_2*x + p_3))*exp(-8419.05535312405/(p_0 + p_1*sin(p_2*x + p_3))),
 1)
In [10]:
integrate_and_plot(autsys)
{'atol': [1e-08], 'fpes': array([], dtype=int64), 'mode': 'adaptive', 'n_dls_jac_evals': 8, 'n_dls_rhs_evals': 0, 'n_err_test_fails': 13, 'n_lin_solv_setups': 44, 'n_nonlin_solv_conv_fails': 0, 'n_nonlin_solv_iters': 497, 'n_rhs_evals': 501, 'n_root_evals': 0, 'n_steps': 419, 'nfev': 501, 'njev': 8, 'njvev': 0, 'orders': array([], dtype=int64), 'rtol': 1e-08, 'steps': array([], dtype=float64), 'success': True, 'time_cpu': 0.1864, 'time_jac': 0.007753973999999999, 'time_jtimes': 0.0, 'time_prec': 0.0, 'time_quads': 0.0, 'time_rhs': 0.17253364099999996, 'time_roots': 0.0, 'time_wall': 0.184317639}