In [ ]:
from collections import defaultdict
import numpy as np
from IPython.display import Latex
import matplotlib.pyplot as plt
from pyodesys.symbolic import SymbolicSys
from chempy import Substance, ReactionSystem
from chempy.kinetics.ode import get_odesys
from chempy.units import SI_base_registry, default_units as u
from chempy.util.graph import rsys2graph
%matplotlib inline
In [ ]:
rsys = ReactionSystem.from_string("A -> B; 'k'", substance_factory=lambda k: Substance(k))
rsys
In [ ]:
odesys, extra = get_odesys(rsys, include_params=False)
odesys.names, odesys.param_names
We can change the expressions of the ODE system manually to account for source and sink terms from the flux:
In [ ]:
t, c1, c2, IA, IB, f, k, fc_A, fc_B = map(odesys.be.Symbol, 't c1 c2 I_A I_B f k phi_A phi_B'.split())
newp = f, fc_A, fc_B
c_feed = {'A': fc_A, 'B': fc_B}
cstr = SymbolicSys(
[
(dep, expr - f*dep + f*c_feed[name]) for name, dep, expr
in zip(odesys.names, odesys.dep, odesys.exprs)
],
params=list(odesys.params) + list(newp),
names=odesys.names,
param_names=list(odesys.param_names) + [p.name for p in newp],
dep_by_name=True,
par_by_name=True,
)
Latex('$' + r'\\ '.join(map(lambda x: ':'.join(x), zip(map(lambda x: r'\frac{d%s}{dt}' % x, cstr.names),
map(cstr.be.latex, cstr.exprs)))) + '$')
In [ ]:
cstr.param_names
In [ ]:
init_c, pars = {'A': .15, 'B': .1}, {'k': 0.8, 'f': 0.3, 'phi_A': .7, 'phi_B': .1}
res = cstr.integrate(10, init_c, pars, integrator='cvode')
In [ ]:
res.plot()
We can derive an analytic solution to the ODE system:
In [ ]:
k = cstr.params[0]
e = cstr.be.exp
exprs = [
fc_A*f/(f + k) + c1 * e(-t*(f + k)),
(fc_A*k + fc_B*(f + k))/(f + k) - c1*e(-f*t)*(e(-t*k) - 1) + c2*e(-f*t)
]
cstr.be.init_printing()
exprs
In [ ]:
exprs0 = [expr.subs(t, 0) for expr in exprs]
exprs0
In [ ]:
sol = cstr.be.solve([expr - c0 for expr, c0 in zip(exprs0, (IA, IB))], (c1, c2))
sol
In [ ]:
exprs2 = [expr.subs(sol) for expr in exprs]
exprs2
In [ ]:
IA
In [ ]:
import sympy as sp
cses, expr_cse = sp.cse([expr.subs({fc_A: sp.Symbol('fr'), fc_B: sp.Symbol('fp'), f: sp.Symbol('fv'),
IA: sp.Symbol('r'), IB: sp.Symbol('p')}) for expr in exprs2])
s = '\n'.join(['%s = %s' % (lhs, rhs) for lhs, rhs in cses] + [str(tuple(expr_cse))])
print(s.replace('exp', 'be.exp').replace('\n(', '\nreturn ('))
In [ ]:
exprs2_0 = [expr.subs(t, 0).simplify() for expr in exprs2]
exprs2_0
In [ ]:
_cb = cstr.be.lambdify([t, IA, IB, k, f, fc_A, fc_B], exprs2)
def analytic(x, c0, params):
return _cb(x, c0['A'], c0['B'], params['k'], params['f'], params['phi_A'], params['phi_B'])
In [ ]:
def get_ref(result, parameters=None):
drctn = -1 if result.odesys.names[0] == 'B' else 1
return np.array(analytic(
result.xout,
{k: result.named_dep(k)[0] for k in result.odesys.names},
parameters or {k: result.named_param(k) for k in result.odesys.param_names})).T[:, ::drctn]
In [ ]:
yref = get_ref(res)
yref.shape
Plotting the error (by comparison to the analytic solution):
In [ ]:
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
res.plot(ax=axes[0])
res.plot(ax=axes[0], y=yref)
res.plot(ax=axes[1], y=res.yout - yref)
In [ ]:
cstr2, extra2 = get_odesys(rsys, include_params=False, cstr=True)
cstr2.exprs
Note how we only needed to pass cstr=True to get_odesys to get the right expressions.
In [ ]:
cstr2.param_names
In [ ]:
renamed_pars = {'k': pars['k'], 'fc_A': pars['phi_A'], 'fc_B': pars['phi_B'], 'feedratio': pars['f']}
res2 = cstr2.integrate(10, init_c, renamed_pars, integrator='cvode')
ref2 = get_ref(res2, pars)
res2.plot(y=res2.yout - ref2)
The analytic solution of a unary reaction in a CSTR is also available in ChemPy:
In [ ]:
from chempy.kinetics.integrated import unary_irrev_cstr
help(unary_irrev_cstr)
The symbols of the feedratio and feed concentrations are available in the second output of get_odesys (the extra dictionary):
In [ ]:
fr, fc = extra2['cstr_fr_fc']
In [ ]:
def get_ref2(result):
drctn = -1 if result.odesys.names[0] == 'B' else 1
return np.array(unary_irrev_cstr(
result.xout,
result.named_param('k'),
result.named_dep('A')[0],
result.named_dep('B')[0],
result.named_param(fc['A']),
result.named_param(fc['B']),
result.named_param(fr))).T[:, ::drctn]
In [ ]:
res2.plot(y=res2.yout - get_ref2(res2))
In [ ]:
assert np.allclose(res2.yout, get_ref2(res2))
In [ ]: