1# -*- coding: utf-8 -*- 2 3from collections import defaultdict, OrderedDict 4from itertools import permutations 5import math 6 7import pytest 8 9try: 10 import numpy as np 11except ImportError: 12 np = None 13 14 15from chempy import Equilibrium, Reaction, ReactionSystem, Substance 16from chempy.thermodynamics.expressions import MassActionEq 17from chempy.units import ( 18 SI_base_registry, 19 get_derived_unit, 20 allclose, 21 units_library, 22 linspace, 23 to_unitless, 24 default_constants as const, 25 default_units as u, 26) 27from chempy.util._expr import Expr 28from chempy.util.testing import requires 29from .test_rates import _get_SpecialFraction_rsys 30from ..arrhenius import ArrheniusParam 31from ..rates import Arrhenius, MassAction, Radiolytic, RampedTemp 32from .._rates import ShiftedTPoly 33from ..ode import ( 34 get_odesys, 35 chained_parameter_variation, 36 _mk_dedim, 37 _create_odesys as create_odesys, 38) 39from ..integrated import dimerization_irrev, binary_rev 40 41 42@requires("numpy", "pyodesys") 43def test_get_odesys_1(): 44 k = 0.2 45 a = Substance("A") 46 b = Substance("B") 47 r = Reaction({"A": 1}, {"B": 1}, param=k) 48 rsys = ReactionSystem([r], [a, b]) 49 assert sorted(rsys.substances.keys()) == ["A", "B"] 50 odesys = get_odesys(rsys, include_params=True)[0] 51 c0 = { 52 "A": 1.0, 53 "B": 3.0, 54 } 55 t = np.linspace(0.0, 10.0) 56 xout, yout, info = odesys.integrate(t, c0) 57 yref = np.zeros((t.size, 2)) 58 yref[:, 0] = np.exp(-k * t) 59 yref[:, 1] = 4 - np.exp(-k * t) 60 assert np.allclose(yout, yref) 61 62 63@requires("numpy", "pyodesys", "sympy") 64def test_get_odesys__rate_exprs_cb(): 65 k = 0.2 66 a = Substance("A") 67 b = Substance("B") 68 r = Reaction({"A": 1}, {"B": 1}, param=k) 69 rsys = ReactionSystem([r], [a, b]) 70 assert sorted(rsys.substances.keys()) == ["A", "B"] 71 odesys, extra = get_odesys(rsys) 72 c0 = {"A": 1.0, "B": 3.0} 73 t = np.linspace(0.0, 10.0) 74 res = odesys.integrate(t, c0) 75 yref = np.zeros((t.size, 2)) 76 yref[:, 0] = np.exp(-k * t) 77 yref[:, 1] = 4 - np.exp(-k * t) 78 assert np.allclose(res.yout, yref) 79 rate = extra["rate_exprs_cb"](res.xout, res.yout, res.params) 80 assert np.allclose(rate[:, 0], k * yref[:, 0]) 81 82 83@requires("numpy", "pyodesys") 84def test_get_odesys_2(): 85 g = Radiolytic([3.14]) 86 a = Substance("A") 87 b = Substance("B") 88 r = Reaction({"A": 1}, {"B": 1}, param=g) 89 rsys = ReactionSystem([r], [a, b]) 90 odesys = get_odesys(rsys, include_params=True)[0] 91 c0 = { 92 "A": 1.0, 93 "B": 3.0, 94 } 95 t = np.linspace(0.0, 0.1) 96 xout, yout, info = odesys.integrate( 97 t, rsys.as_per_substance_array(c0), {"doserate": 2.72, "density": 0.998} 98 ) 99 yref = np.zeros((t.size, 2)) 100 k = 3.14 * 2.72 * 0.998 101 yref[:, 0] = 1 - k * t 102 yref[:, 1] = 3 + k * t 103 assert np.allclose(yout, yref) 104 105 106@requires(units_library, "pyodesys") 107def test_get_odesys_3(): 108 M = u.molar 109 s = u.second 110 mol = u.mol 111 m = u.metre 112 substances = list(map(Substance, "H2O H+ OH-".split())) 113 dissociation = Reaction({"H2O": 1}, {"H+": 1, "OH-": 1}, 2.47e-5 / s) 114 recombination = Reaction({"H+": 1, "OH-": 1}, {"H2O": 1}, 1.37e11 / M / s) 115 rsys = ReactionSystem([dissociation, recombination], substances) 116 odesys = get_odesys( 117 rsys, include_params=True, unit_registry=SI_base_registry, output_conc_unit=M 118 )[0] 119 c0 = {"H2O": 55.4 * M, "H+": 1e-7 * M, "OH-": 1e-4 * mol / m ** 3} 120 x, y, p = odesys.to_arrays( 121 -42 * u.second, rsys.as_per_substance_array(c0, unit=M), () 122 ) 123 fout = odesys.f_cb(x, y, p) 124 125 time_unit = get_derived_unit(SI_base_registry, "time") 126 conc_unit = get_derived_unit(SI_base_registry, "concentration") 127 128 r1 = to_unitless(55.4 * 2.47e-5 * M / s, conc_unit / time_unit) 129 r2 = to_unitless(1e-14 * 1.37e11 * M / s, conc_unit / time_unit) 130 assert np.all(abs(fout[:, 0] - r2 + r1)) < 1e-10 131 assert np.all(abs(fout[:, 1] - r1 + r2)) < 1e-10 132 assert np.all(abs(fout[:, 2] - r1 + r2)) < 1e-10 133 134 135@requires(units_library, "pyodesys") 136def test_get_odesys__with_units(): 137 a = Substance("A") 138 b = Substance("B") 139 molar = u.molar 140 second = u.second 141 r = Reaction({"A": 2}, {"B": 1}, param=1e-3 / molar / second) 142 rsys = ReactionSystem([r], [a, b]) 143 odesys = get_odesys(rsys, include_params=True, unit_registry=SI_base_registry)[0] 144 c0 = {"A": 13 * u.mol / u.metre ** 3, "B": 0.2 * u.molar} 145 conc_unit = get_derived_unit(SI_base_registry, "concentration") 146 t = np.linspace(0, 10) * u.hour 147 xout, yout, info = odesys.integrate( 148 t, rsys.as_per_substance_array(c0, unit=conc_unit), atol=1e-10, rtol=1e-12 149 ) 150 151 t_unitless = to_unitless(xout, u.second) 152 Aref = dimerization_irrev(t_unitless, 1e-6, 13.0) 153 # Aref = 1/(1/13 + 2*1e-6*t_unitless) 154 yref = np.zeros((xout.size, 2)) 155 yref[:, 0] = Aref 156 yref[:, 1] = 200 + (13 - Aref) / 2 157 assert allclose(yout, yref * conc_unit) 158 159 160@requires("numpy", "pyodesys") 161def test_SpecialFraction(): 162 k, kprime = 3.142, 2.718 163 rsys = _get_SpecialFraction_rsys(k, kprime) 164 165 odesys = get_odesys(rsys, include_params=True)[0] 166 c0 = {"H2": 13, "Br2": 17, "HBr": 19} 167 r = k * c0["H2"] * c0["Br2"] ** (3 / 2) / (c0["Br2"] + kprime * c0["HBr"]) 168 ref = rsys.as_per_substance_array({"H2": -r, "Br2": -r, "HBr": 2 * r}) 169 res = odesys.f_cb(0, rsys.as_per_substance_array(c0)) 170 assert np.allclose(res, ref) 171 172 173@requires(units_library, "pyodesys") 174def test_SpecialFraction_with_units(): 175 k, kprime = 3.142 * u.s ** -1 * u.molar ** -0.5, 2.718 176 rsys = _get_SpecialFraction_rsys(k, kprime) 177 odesys = get_odesys(rsys, include_params=True, unit_registry=SI_base_registry)[0] 178 c0 = {"H2": 13 * u.molar, "Br2": 16 * u.molar, "HBr": 19 * u.molar} 179 r = k * c0["H2"] * c0["Br2"] ** (3 / 2) / (c0["Br2"] + kprime * c0["HBr"]) 180 conc_unit = u.mol / u.metre ** 3 181 rate_unit = conc_unit / u.second 182 ref = rsys.as_per_substance_array( 183 {"H2": -r, "Br2": -r, "HBr": 2 * r}, unit=rate_unit 184 ) 185 res = odesys.f_cb(0, rsys.as_per_substance_array(c0, unit=conc_unit)) 186 assert allclose(to_unitless(ref, rate_unit), res) 187 188 189@requires("pyodesys") 190def test_ode_with_global_parameters(): 191 ratex = MassAction(Arrhenius([1e10, 40e3 / 8.3145])) 192 rxn = Reaction({"A": 1}, {"B": 1}, ratex) 193 rsys = ReactionSystem([rxn], "A B") 194 odesys, extra = get_odesys(rsys, include_params=False) 195 param_keys, unique_keys, p_units = map( 196 extra.get, "param_keys unique p_units".split() 197 ) 198 conc = {"A": 3, "B": 5} 199 x, y, p = odesys.to_arrays(-37, conc, {"temperature": 298.15}) 200 fout = odesys.f_cb(x, y, p) 201 ref = 3 * 1e10 * np.exp(-40e3 / 8.3145 / 298.15) 202 assert np.all(abs((fout[:, 0] + ref) / ref) < 1e-14) 203 assert np.all(abs((fout[:, 1] - ref) / ref) < 1e-14) 204 205 206@requires("pyodesys") 207def test_get_ode__ArrheniusParam(): 208 rxn = Reaction({"A": 1}, {"B": 1}, None) 209 rxn.param = ArrheniusParam(1e10, 40e3) 210 rsys = ReactionSystem([rxn], "A B") 211 odesys = get_odesys(rsys, include_params=True)[0] 212 conc = {"A": 3, "B": 5} 213 x, y, p = odesys.to_arrays(-37, conc, {"temperature": 200}) 214 fout = odesys.f_cb(x, y, p) 215 ref = 3 * 1e10 * np.exp(-40e3 / 8.314472 / 200) 216 assert np.all(abs((fout[:, 0] + ref) / ref) < 1e-14) 217 assert np.all(abs((fout[:, 1] - ref) / ref) < 1e-14) 218 219 220@requires("pyodesys") 221def test_get_ode__Radiolytic(): 222 rad = Radiolytic([2.4e-7]) 223 rxn = Reaction({"A": 4, "B": 1}, {"C": 3, "D": 2}, rad) 224 rsys = ReactionSystem([rxn], "A B C D") 225 odesys = get_odesys(rsys, include_params=True)[0] 226 c = {"A": 3, "B": 5, "C": 11, "D": 13} 227 x, y, p = odesys.to_arrays(-37, c, {"doserate": 0.4, "density": 0.998}) 228 fout = odesys.f_cb(x, y, p) 229 r = 2.4e-7 * 0.4 * 0.998 230 ref = [-4 * r, -r, 3 * r, 2 * r] 231 assert np.all(abs((fout - ref) / ref) < 1e-14) 232 233 234@requires("pyodesys", units_library) 235def test_get_ode__Radiolytic__units(): 236 rad = Radiolytic([2.4e-7 * u.mol / u.joule]) 237 rxn = Reaction({"A": 4, "B": 1}, {"C": 3, "D": 2}, rad) 238 rsys = ReactionSystem([rxn], "A B C D") 239 odesys = get_odesys(rsys, include_params=True, unit_registry=SI_base_registry)[0] 240 conc = {"A": 3 * u.molar, "B": 5 * u.molar, "C": 11 * u.molar, "D": 13 * u.molar} 241 x, y, p = odesys.to_arrays( 242 -37 * u.second, 243 conc, 244 { 245 "doserate": 0.4 * u.gray / u.second, 246 "density": 0.998 * u.kg / u.decimetre ** 3, 247 }, 248 ) 249 fout = odesys.f_cb(x, y, p) # f_cb does not carry any units 250 r = 2.4e-7 * 0.4 * 0.998 * 1e3 # mol/m3 251 ref = [-4 * r, -r, 3 * r, 2 * r] 252 assert np.all(abs((fout - ref) / ref) < 1e-14) 253 254 255@requires("pyodesys", units_library) 256def test_get_ode__Radiolytic__units__multi(): 257 rad = Radiolytic([2.4e-7 * u.mol / u.joule]) 258 rxn = Reaction({"A": 4, "B": 1}, {"C": 3, "D": 2}, rad) 259 rsys = ReactionSystem([rxn], "A B C D") 260 odesys = get_odesys(rsys, include_params=True, unit_registry=SI_base_registry)[0] 261 conc = {"A": 3 * u.molar, "B": 5 * u.molar, "C": 11 * u.molar, "D": 13 * u.molar} 262 doserates = [dr * u.gray / u.second for dr in [0.1, 0.2, 0.3, 0.4]] 263 results = odesys.integrate( 264 37 * u.second, 265 conc, 266 {"doserate": doserates, "density": 0.998 * u.kg / u.decimetre ** 3}, 267 ) 268 assert len(results) == 4 269 for i, r in enumerate(results): 270 dr = r.params[odesys.param_names.index("doserate")] 271 assert dr.ndim == 0 or len(dr) == 1 272 assert dr == doserates[i] 273 274 275class Density(Expr): 276 """Arguments: rho0 drhodT T0""" 277 278 parameter_keys = ("temperature",) 279 kw = {} 280 281 def __call__(self, variables, backend=None): 282 rho0, drhodT, T0 = self.all_args(variables) 283 return rho0 + drhodT * (variables["temperature"] - T0) 284 285 286@requires("pyodesys") 287def test_get_ode__Radiolytic__substitutions(): 288 rad = Radiolytic([2.4e-7]) 289 rxn = Reaction({"A": 4, "B": 1}, {"C": 3, "D": 2}, rad) 290 rsys = ReactionSystem([rxn], "A B C D") 291 substance_rho = Density([1, -1e-3, 273.15]) 292 odesys = get_odesys( 293 rsys, include_params=True, substitutions={"density": substance_rho} 294 )[0] 295 conc = {"A": 3, "B": 5, "C": 11, "D": 13} 296 state = {"doserate": 0.4, "temperature": 298.15} 297 x, y, p = odesys.to_arrays(-37, conc, state) 298 fout = odesys.f_cb(x, y, p) 299 r = 2.4e-7 * 0.4 * substance_rho({"temperature": 298.15}) 300 ref = [-4 * r, -r, 3 * r, 2 * r] 301 assert np.all(abs((fout - ref) / ref) < 1e-14) 302 303 304@requires("pyodesys", units_library) 305def test_get_ode__Radiolytic__substitutions__units(): 306 rad = Radiolytic([2.4e-7 * u.mol / u.joule]) 307 rxn = Reaction({"A": 4, "B": 1}, {"C": 3, "D": 2}, rad) 308 rsys = ReactionSystem([rxn], "A B C D") 309 g_dm3 = u.gram / u.decimetre ** 3 310 kg_dm3 = u.kg / u.decimetre ** 3 311 substance_rho = Density([1 * kg_dm3, -1 * g_dm3 / u.kelvin, 273.15 * u.kelvin]) 312 odesys = get_odesys( 313 rsys, 314 include_params=True, 315 unit_registry=SI_base_registry, 316 substitutions={"density": substance_rho}, 317 )[0] 318 conc = {"A": 3 * u.molar, "B": 5 * u.molar, "C": 11 * u.molar, "D": 13 * u.molar} 319 x, y, p = odesys.to_arrays( 320 -37 * u.second, 321 conc, 322 {"doserate": 0.4 * u.gray / u.second, "temperature": 298.15 * u.kelvin}, 323 ) 324 fout = odesys.f_cb(x, y, p) 325 r = 2.4e-7 * 0.4 * 0.975 * 1e3 # mol/m3/s 326 ref = [-4 * r, -r, 3 * r, 2 * r] 327 assert np.all(abs((fout - ref) / ref) < 1e-14) 328 329 330@requires("pyodesys", units_library) 331def test_get_ode__TPoly(): 332 rate = MassAction( 333 ShiftedTPoly([273.15 * u.K, 10 / u.molar / u.s, 2 / u.molar / u.s / u.K]) 334 ) 335 rxn = Reaction({"A": 1, "B": 1}, {"C": 3, "D": 2}, rate, {"A": 3}) 336 rsys = ReactionSystem([rxn], "A B C D") 337 odesys = get_odesys(rsys, unit_registry=SI_base_registry)[0] 338 conc = {"A": 3 * u.molar, "B": 5 * u.molar, "C": 11 * u.molar, "D": 13 * u.molar} 339 x, y, p = odesys.to_arrays(-37 * u.second, conc, {"temperature": 298.15 * u.kelvin}) 340 fout = odesys.f_cb(x, y, p) 341 r = 3 * 5 * (10 + 2 * 25) * 1000 # mol/m3/s 342 ref = [-4 * r, -r, 3 * r, 2 * r] 343 assert np.all(abs((fout - ref) / ref) < 1e-14) 344 345 346@requires("pyodesys", units_library) 347def test_get_odesys__time_dep_rate(): 348 class RampedRate(Expr): 349 argument_names = ("rate_constant", "ramping_rate") 350 351 def __call__(self, variables, reaction, backend=math): 352 rate_constant, ramping_rate = self.all_args(variables, backend=backend) 353 return rate_constant * ramping_rate * variables["time"] 354 355 rate = MassAction(RampedRate([7, 2])) 356 rxn = Reaction({"A": 1}, {"B": 3}, rate) 357 rsys = ReactionSystem([rxn], "A B") 358 odesys = get_odesys(rsys)[0] 359 conc = {"A": 3, "B": 11} 360 x, y, p = odesys.to_arrays([5, 13, 17], conc, ()) 361 fout = odesys.f_cb(x, y, p) 362 r = 2 * 7 * 3 363 ref = np.array([[-r * 5, -r * 13, -r * 17], [r * 5 * 3, r * 13 * 3, r * 17 * 3]]).T 364 assert np.allclose(fout, ref) 365 366 367@requires("pyodesys", units_library) 368def test_get_odesys__time_dep_temperature(): 369 import sympy as sp 370 371 def refA(t, A0, A, Ea_over_R, T0, dTdt): 372 T = T0 + dTdt * t 373 d_Ei = sp.Ei(-Ea_over_R / T0).n(100).round(90) - sp.Ei(-Ea_over_R / T).n( 374 100 375 ).round(90) 376 d_Texp = T0 * sp.exp(-Ea_over_R / T0) - T * sp.exp(-Ea_over_R / T) 377 return A0 * sp.exp(A / dTdt * (Ea_over_R * d_Ei + d_Texp)).n(30) 378 379 params = A0, A, Ea_over_R, T0, dTdt = [13, 1e10, 56e3 / 8, 273, 2] 380 B0 = 11 381 rate = MassAction(Arrhenius([A, Ea_over_R])) 382 rxn = Reaction({"A": 1}, {"B": 3}, rate) 383 rsys = ReactionSystem([rxn], "A B") 384 rt = RampedTemp([T0, dTdt], ("init_temp", "ramp_rate")) 385 odesys, extra = get_odesys(rsys, False, substitutions={"temperature": rt}) 386 all_pk, unique, p_units = map(extra.get, "param_keys unique p_units".split()) 387 conc = {"A": A0, "B": B0} 388 tout = [2, 5, 10] 389 390 for ramp_rate in [2, 3, 4]: 391 unique["ramp_rate"] = ramp_rate 392 xout, yout, info = odesys.integrate(10, conc, unique, atol=1e-10, rtol=1e-12) 393 params[-1] = ramp_rate 394 Aref = np.array([float(refA(t, *params)) for t in xout]) 395 # Aref = 1/(1/13 + 2*1e-6*t_unitless) 396 yref = np.zeros((xout.size, 2)) 397 yref[:, 0] = Aref 398 yref[:, 1] = B0 + 3 * (A0 - Aref) 399 assert allclose(yout, yref) 400 401 unique["ramp_rate"] = 2 402 x, y, p = odesys.to_arrays(tout, conc, unique) 403 fout = odesys.f_cb(x, y, p) 404 405 def r(t): 406 return A * np.exp(-Ea_over_R / (T0 + dTdt * t)) * A0 # refA(t, *params) 407 408 ref = np.array([[-r(2), -r(5), -r(10)], [3 * r(2), 3 * r(5), 3 * r(10)]]).T 409 assert np.allclose(fout, ref) 410 411 412@requires("numpy", "pyodesys") 413def test_get_odesys__late_binding(): 414 def _gibbs(args, T, R, backend, **kwargs): 415 H, S = args 416 return backend.exp(-(H - T * S) / (R * T)) 417 418 def _eyring(args, T, R, k_B, h, backend, **kwargs): 419 H, S = args 420 return k_B / h * T * backend.exp(-(H - T * S) / (R * T)) 421 422 gibbs_pk = ("temperature", "molar_gas_constant") 423 eyring_pk = gibbs_pk + ("Boltzmann_constant", "Planck_constant") 424 425 GibbsEC = MassActionEq.from_callback( 426 _gibbs, argument_names=("H", "S"), parameter_keys=gibbs_pk 427 ) 428 EyringMA = MassAction.from_callback( 429 _eyring, argument_names=("H", "S"), parameter_keys=eyring_pk 430 ) 431 432 uk_equil = ("He_assoc", "Se_assoc") 433 beta = GibbsEC(unique_keys=uk_equil) # equilibrium parameters 434 435 uk_kinet = ("Ha_assoc", "Sa_assoc") 436 bimol_barrier = EyringMA(unique_keys=uk_kinet) # activation parameters 437 438 eq = Equilibrium({"Fe+3", "SCN-"}, {"FeSCN+2"}, beta) 439 rsys = ReactionSystem(eq.as_reactions(kf=bimol_barrier)) 440 odesys, extra = get_odesys(rsys, include_params=False) 441 pk, unique, p_units = map(extra.get, "param_keys unique p_units".split()) 442 assert sorted(unique) == sorted(uk_equil + uk_kinet) 443 assert sorted(pk) == sorted(eyring_pk) 444 445 446@requires("numpy", "pyodesys") 447def test_get_odesys__ScaledSys(): 448 from pyodesys.symbolic import ScaledSys 449 450 k = 0.2 451 a = Substance("A") 452 b = Substance("B") 453 r = Reaction({"A": 1}, {"B": 1}, param=k) 454 rsys = ReactionSystem([r], [a, b]) 455 assert sorted(rsys.substances.keys()) == ["A", "B"] 456 odesys = get_odesys(rsys, include_params=True, SymbolicSys=ScaledSys)[0] 457 c0 = { 458 "A": 1.0, 459 "B": 3.0, 460 } 461 t = np.linspace(0.0, 10.0) 462 xout, yout, info = odesys.integrate(t, c0) 463 yref = np.zeros((t.size, 2)) 464 yref[:, 0] = np.exp(-k * t) 465 yref[:, 1] = 4 - np.exp(-k * t) 466 assert np.allclose(yout, yref) 467 468 469@requires("numpy", "pyodesys", "sympy") 470def test_get_odesys__max_euler_step_cb(): 471 rsys = ReactionSystem.from_string( 472 "\n".join(["H2O -> H+ + OH-; 1e-4", "OH- + H+ -> H2O; 1e10"]) 473 ) 474 odesys, extra = get_odesys(rsys) 475 r1 = 1.01e-4 476 r2 = 6e-4 477 dH2Odt = r2 - r1 478 euler_ref = 2e-7 / dH2Odt 479 assert ( 480 abs( 481 extra["max_euler_step_cb"](0, {"H2O": 1.01, "H+": 2e-7, "OH-": 3e-7}) 482 - euler_ref 483 ) 484 / euler_ref 485 < 1e-8 486 ) 487 488 489@requires("numpy", "pyodesys", "sympy") 490@pytest.mark.parametrize("substances", permutations(["H2O", "H+", "OH-"])) 491def test_get_odesys__linear_dependencies__preferred(substances): 492 rsys = ReactionSystem.from_string( 493 "\n".join(["H2O -> H+ + OH-; 1e-4", "OH- + H+ -> H2O; 1e10"]), substances 494 ) 495 assert isinstance(rsys.substances, OrderedDict) 496 odesys, extra = get_odesys(rsys) 497 498 af_H2O_H = extra["linear_dependencies"](["H+", "H2O"]) 499 import sympy 500 501 y0 = {k: sympy.Symbol(k + "0") for k in rsys.substances} 502 af_H2O_H( 503 None, {odesys[k]: v for k, v in y0.items()}, None, sympy 504 ) # ensure idempotent 505 exprs_H2O_H = af_H2O_H(None, {odesys[k]: v for k, v in y0.items()}, None, sympy) 506 ref_H2O_H = { 507 "H2O": y0["H2O"] + y0["OH-"] - odesys["OH-"], # oxygen 508 "H+": 2 * y0["H2O"] 509 + y0["H+"] 510 + y0["OH-"] 511 - odesys["OH-"] 512 - 2 * (y0["H2O"] + y0["OH-"] - odesys["OH-"]), # hydrogen 513 } 514 for k, v in ref_H2O_H.items(): 515 assert (exprs_H2O_H[odesys[k]] - v) == 0 516 517 518@requires("numpy", "pyodesys", "sympy", "pycvodes") 519@pytest.mark.parametrize( 520 "preferred", [None, ["H+", "OH-"], ["H2O", "H+"], ["H2O", "OH-"]] 521) 522def test_get_odesys__linear_dependencies__PartiallySolvedSystem(preferred): 523 import sympy 524 from pyodesys.symbolic import PartiallySolvedSystem 525 526 rsys = ReactionSystem.from_string( 527 "\n".join(["H2O -> H+ + OH-; 1e-4", "OH- + H+ -> H2O; 1e10"]) 528 ) 529 odesys, extra = get_odesys(rsys) 530 c0 = {"H2O": 0, "H+": 2e-7, "OH-": 3e-7} 531 h0max = extra["max_euler_step_cb"](0, c0) 532 analytic_factory = extra["linear_dependencies"]() 533 y0 = {k: sympy.Symbol(k + "0") for k in rsys.substances} 534 analytic_factory(None, {odesys[k]: v for k, v in y0.items()}, None, sympy) 535 psys = PartiallySolvedSystem(odesys, analytic_factory) 536 xout, yout, info = psys.integrate( 537 1, c0, atol=1e-12, rtol=1e-10, first_step=h0max * 1e-12, integrator="cvode" 538 ) 539 c_reac = c0["H+"], c0["OH-"] 540 H2O_ref = binary_rev(xout, 1e10, 1e-4, c0["H2O"], max(c_reac), min(c_reac)) 541 assert np.allclose(yout[:, psys.names.index("H2O")], H2O_ref) 542 assert np.allclose(yout[:, psys.names.index("H+")], c0["H+"] + c0["H2O"] - H2O_ref) 543 assert np.allclose( 544 yout[:, psys.names.index("OH-")], c0["OH-"] + c0["H2O"] - H2O_ref 545 ) 546 547 548@requires("numpy", "pyodesys", "sympy", "pycvodes") 549def test_get_odesys__Equilibrium_as_reactions(): 550 from chempy import Equilibrium, ReactionSystem 551 552 eq = Equilibrium({"Fe+3", "SCN-"}, {"FeSCN+2"}, 10 ** 2) 553 substances = "Fe+3 SCN- FeSCN+2".split() 554 rsys = ReactionSystem(eq.as_reactions(kf=3.0), substances) 555 odesys, extra = get_odesys(rsys) 556 init_conc = {"Fe+3": 1.0, "SCN-": 0.3, "FeSCN+2": 0} 557 tout, Cout, info = odesys.integrate( 558 5, init_conc, integrator="cvode", atol=1e-11, rtol=1e-12 559 ) 560 cmplx_ref = binary_rev( 561 tout, 3, 3.0 / 100, init_conc["FeSCN+2"], init_conc["Fe+3"], init_conc["SCN-"] 562 ) 563 assert np.allclose(Cout[:, 2], cmplx_ref) 564 565 566@requires("numpy", "pyodesys", "sympy", "pycvodes") 567def test_get_odesys__Expr_as_param(): 568 def _eyring_pe(args, T, backend=math, **kwargs): 569 (freq,) = args 570 return freq * T 571 572 EyringPreExp = Expr.from_callback( 573 _eyring_pe, argument_names=("freq",), parameter_keys=("temperature",) 574 ) 575 576 def _k(args, T, backend=math, **kwargs): 577 A, H, S = args 578 return A * backend.exp(-(H - T * S) / (8.314511 * T)) 579 580 EyringMA = MassAction.from_callback( 581 _k, parameter_keys=("temperature",), argument_names=("Aa", "Ha", "Sa") 582 ) 583 kb_h = 2.08e10 584 rxn = Reaction({"A"}, {"B"}, EyringMA(unique_keys=("A_u", "H_u", "S_u"))) 585 rsys = ReactionSystem([rxn], ["A", "B"]) 586 odesys, extra = get_odesys( 587 rsys, include_params=False, substitutions={"A_u": EyringPreExp(kb_h)} 588 ) 589 y0 = defaultdict(float, {"A": 7.0}) 590 rt = 293.15 591 xout, yout, info = odesys.integrate( 592 5, 593 y0, 594 {"H_u": 117e3, "S_u": 150, "temperature": rt}, 595 integrator="cvode", 596 atol=1e-12, 597 rtol=1e-10, 598 nsteps=1000, 599 ) 600 kref = kb_h * rt * np.exp(-(117e3 - rt * 150) / (8.314511 * rt)) 601 ref = y0["A"] * np.exp(-kref * xout) 602 assert np.allclose(yout[:, 0], ref) 603 assert np.allclose(yout[:, 1], y0["A"] - ref) 604 605 606@requires("numpy", "pyodesys", "sympy", "pycvodes") 607def test_get_odesys__Expr_as_param__unique_as_param(): 608 def _eyring_pe_coupled(args, T, S, backend=math, **kwargs): 609 (freq,) = args 610 return freq * T / S 611 612 EyringPreExpCoupled = Expr.from_callback( 613 _eyring_pe_coupled, 614 argument_names=("freq",), 615 parameter_keys=("temperature", "S_u"), 616 ) 617 618 def _k(args, T, backend=math, **kwargs): 619 A, H, S = args 620 return A * backend.exp(-(H - T * S) / (8.314511 * T)) 621 622 EyringMA = MassAction.from_callback( 623 _k, parameter_keys=("temperature",), argument_names=("Aa", "Ha", "Sa") 624 ) 625 kb_h = 2.08e10 626 rxn = Reaction({"A"}, {"B"}, EyringMA(unique_keys=("A_u", "H_u", "S_u"))) 627 rsys = ReactionSystem([rxn], ["A", "B"]) 628 odesys2, extra2 = get_odesys( 629 rsys, include_params=False, substitutions={"A_u": EyringPreExpCoupled(kb_h)} 630 ) 631 y0 = defaultdict(float, {"A": 7.0}) 632 rt = 293.15 633 xout2, yout2, info2 = odesys2.integrate( 634 5, 635 y0, 636 {"H_u": 107e3, "S_u": 150, "temperature": rt}, 637 integrator="cvode", 638 atol=1e-12, 639 rtol=1e-10, 640 nsteps=1000, 641 ) 642 kref2 = kb_h * rt * np.exp(-(107e3 - rt * 150) / (8.314511 * rt)) / 150 643 ref2 = y0["A"] * np.exp(-kref2 * xout2) 644 assert np.allclose(yout2[:, 0], ref2) 645 assert np.allclose(yout2[:, 1], y0["A"] - ref2) 646 647 648@requires("pyodesys", "pycvodes") 649def test_chained_parameter_variation(): 650 ratex = MassAction(Arrhenius([1e10, 63e3 / 8.3145])) 651 rxn = Reaction({"A": 1}, {"B": 1}, ratex) 652 rsys = ReactionSystem([rxn], "A B") 653 odesys, extra = get_odesys(rsys, include_params=False) 654 param_keys, unique_keys, p_units = map( 655 extra.get, "param_keys unique p_units".split() 656 ) 657 conc = {"A": 3.17, "B": 5.03} 658 Ts = (294, 304, 317) 659 times = [3.1, 2.1, 5.3] 660 kw = dict(integrator="cvode", atol=1e-12, rtol=1e-13, first_step=1e-14) 661 tout, cout, info = chained_parameter_variation( 662 odesys, times, conc, {"temperature": Ts}, {}, integrate_kwargs=kw 663 ) 664 assert len(info["nfev"]) == 3 665 assert info["nfev"][0] > 2 666 assert info["nfev"][1] > 2 667 assert info["nfev"][2] > 2 668 assert np.all(np.diff(tout) > 0) 669 tout1 = tout[tout <= times[0]] 670 tout23 = tout[tout > times[0]] 671 tout2 = tout23[tout23 <= times[0] + times[1]] 672 tout3 = tout23[tout23 > times[0] + times[1]] 673 674 def _ref(y0, x, T, x0): 675 k = 1e10 * np.exp(-63e3 / 8.3145 / T) 676 return y0 * np.exp(-k * (x - x0)) 677 678 Aref1 = _ref(conc["A"], tout1, Ts[0], tout1[0]) 679 Bref1 = conc["B"] + conc["A"] - Aref1 680 681 Aref2 = _ref(Aref1[-1], tout2, Ts[1], tout1[-1]) 682 Bref2 = Bref1[-1] + Aref1[-1] - Aref2 683 684 Aref3 = _ref(Aref2[-1], tout3, Ts[2], tout2[-1]) 685 Bref3 = Bref2[-1] + Aref2[-1] - Aref3 686 687 cref = np.concatenate( 688 [ 689 np.vstack((a, b)).T 690 for a, b in [(Aref1, Bref1), (Aref2, Bref2), (Aref3, Bref3)] 691 ] 692 ) 693 forgive = 27 * 1.1 694 assert np.allclose(cref, cout, atol=kw["atol"] * forgive, rtol=kw["rtol"] * forgive) 695 696 697def _check_cstr(odesys, fr, fc, extra_pars=None): 698 tout, c0 = np.linspace(0, 0.13, 7), {"H2O2": 2, "O2": 4, "H2O": 3} 699 params = {fr: 13, fc["H2O2"]: 11, fc["O2"]: 43, fc["H2O"]: 45} 700 params.update(extra_pars or {}) 701 res = odesys.integrate(tout, c0, params) 702 from chempy.kinetics.integrated import binary_irrev_cstr 703 704 def get_analytic(result, k, n): 705 ref = binary_irrev_cstr( 706 result.xout, 707 5, 708 result.named_dep("H2O2")[0], 709 result.named_dep(k)[0], 710 result.named_param(fc["H2O2"]), 711 result.named_param(fc[k]), 712 result.named_param(fr), 713 n, 714 ) 715 return np.array(ref).T 716 717 ref_O2 = get_analytic(res, "O2", 1) 718 ref_H2O = get_analytic(res, "H2O", 2) 719 assert np.allclose(res.named_dep("H2O2"), ref_O2[:, 0]) 720 assert np.allclose(res.named_dep("H2O2"), ref_H2O[:, 0]) 721 assert np.allclose(res.named_dep("O2"), ref_O2[:, 1]) 722 assert np.allclose(res.named_dep("H2O"), ref_H2O[:, 1]) 723 724 725@requires("pyodesys", "scipy", "sym") 726def test_get_odesys__cstr(): 727 rsys = ReactionSystem.from_string("2 H2O2 -> O2 + 2 H2O; 5") 728 odesys, extra = get_odesys(rsys, cstr=True) 729 fr, fc = extra["cstr_fr_fc"] 730 _check_cstr(odesys, fr, fc) 731 732 733@requires("pyodesys", "scipy", "sym") 734def test_create_odesys__cstr(): 735 rsys = ReactionSystem.from_string("2 H2O2 -> O2 + 2 H2O; 'k2'") 736 fr, fc = "feedratio", OrderedDict([(sk, "fc_%s" % sk) for sk in rsys.substances]) 737 odesys, extra = create_odesys(rsys, rates_kw=dict(cstr_fr_fc=(fr, fc))) 738 _check_cstr(odesys, fr, fc, extra_pars=dict(k2=5)) 739 740 741@requires("pygslodeiv2", "sym", units_library) 742def test_get_odesys_rsys_with_units(): 743 rsys = ReactionSystem.from_string( 744 """ 745 A -> B; 0.096/s 746 B + C -> P; 4e3/M/s 747 """, 748 substance_factory=Substance, 749 ) 750 with pytest.raises(Exception): 751 get_odesys( 752 rsys 753 ) # not a strict test, SI_base_registry could be made the default 754 755 odesys, extra = get_odesys(rsys, unit_registry=SI_base_registry) 756 tend = 10 757 tend_units = tend * u.s 758 c0 = {"A": 1e-6, "B": 0, "C": 1, "P": 0} 759 c0_units = {k: v * u.molar for k, v in c0.items()} 760 result1 = odesys.integrate(tend_units, c0_units, integrator="gsl") 761 assert result1.info["success"] 762 763 with pytest.raises(Exception): 764 odesys.integrate(tend, c0, integrator="gsl") 765 766 767@requires("pyodeint", "sym", units_library) 768def test_get_odesys_rsys_with_units__named_params(): 769 rsys = ReactionSystem.from_string( 770 """ 771 A -> B; 'k1' 772 B + C -> P; 'k2' 773 """, 774 substance_factory=Substance, 775 ) 776 odesys, extra = get_odesys( 777 rsys, include_params=False, unit_registry=SI_base_registry 778 ) 779 tend = 10 780 tend_units = tend * u.s 781 c0 = {"A": 1e-6, "B": 0, "C": 1, "P": 0} 782 p = {"k1": 3, "k2": 4} 783 p_units = {"k1": 3 / u.s, "k2": 4 / u.M / u.s} 784 c0_units = {k: v * u.molar for k, v in c0.items()} 785 result1 = odesys.integrate(tend_units, c0_units, p_units, integrator="odeint") 786 assert result1.info["success"] 787 788 with pytest.raises(Exception): 789 odesys.integrate(tend, c0, p, integrator="odeint") 790 791 792@requires("pycvodes", "sym", units_library) 793def test_get_odesys__Eyring(): 794 R = 8.314472 795 T_K = 300 796 dH = 80e3 797 dS = 10 798 rsys1 = ReactionSystem.from_string( 799 """ 800 NOBr -> NO + Br; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol) 801 """.format( 802 dH=dH, dS=dS 803 ), 804 substances="NOBr NO Br".split(), 805 ) 806 kref = 20836643994.118652 * T_K * np.exp(-(dH - T_K * dS) / (R * T_K)) 807 NOBr0_M = 0.7 808 init_cond = dict(NOBr=NOBr0_M * u.M, NO=0 * u.M, Br=0 * u.M) 809 t = 5 * u.second 810 params = dict(temperature=T_K * u.K) 811 812 def check(rsys): 813 odesys, extra = get_odesys( 814 rsys, unit_registry=SI_base_registry, constants=const 815 ) 816 res = odesys.integrate(t, init_cond, params, integrator="cvode") 817 NOBr_ref = NOBr0_M * np.exp(-kref * to_unitless(res.xout, u.second)) 818 ref = np.array([NOBr_ref] + [NOBr0_M - NOBr_ref] * 2).T 819 cmp = to_unitless(res.yout, u.M) 820 assert np.allclose(cmp, ref) 821 822 check(rsys1) 823 rsys2 = ReactionSystem.from_string( 824 """ 825 NOBr -> NO + Br; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol])) 826 """.format( 827 dH=dH, dS=dS 828 ), 829 substances="NOBr NO Br".split(), 830 ) 831 check(rsys2) 832 833 834@requires("pycvodes", "sym", units_library) 835def test_get_odesys__Eyring_2nd_order(): 836 R = 8.314472 837 T_K = 300 838 dH = 80e3 839 dS = 10 840 rsys1b = ReactionSystem.from_string( 841 """ 842 NO + Br -> NOBr; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol) 843 """.format( 844 dH=dH, dS=dS 845 ) 846 ) 847 c0 = 1 # mol/dm3 === 1000 mol/m3 848 kbref = 20836643994.118652 * T_K * np.exp(-(dH - T_K * dS) / (R * T_K)) / c0 849 NO0_M = 1.5 850 Br0_M = 0.7 851 init_cond = dict(NOBr=0 * u.M, NO=NO0_M * u.M, Br=Br0_M * u.M) 852 t = 5 * u.second 853 params = dict(temperature=T_K * u.K) 854 855 def analytic_b(t): 856 U, V = NO0_M, Br0_M 857 d = U - V 858 return (U * (1 - np.exp(-kbref * t * d))) / (U / V - np.exp(-kbref * t * d)) 859 860 def check(rsys): 861 odesys, extra = get_odesys( 862 rsys, unit_registry=SI_base_registry, constants=const 863 ) 864 res = odesys.integrate(t, init_cond, params, integrator="cvode") 865 t_sec = to_unitless(res.xout, u.second) 866 NOBr_ref = analytic_b(t_sec) 867 cmp = to_unitless(res.yout, u.M) 868 ref = np.empty_like(cmp) 869 ref[:, odesys.names.index("NOBr")] = NOBr_ref 870 ref[:, odesys.names.index("Br")] = Br0_M - NOBr_ref 871 ref[:, odesys.names.index("NO")] = NO0_M - NOBr_ref 872 assert np.allclose(cmp, ref) 873 874 check(rsys1b) 875 rsys2b = ReactionSystem.from_string( 876 """ 877 NO + Br -> NOBr; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol])) 878 """.format( 879 dH=dH, dS=dS 880 ) 881 ) 882 check(rsys2b) 883 884 885@requires("pycvodes", "sym", "scipy", units_library) 886def test_get_odesys__Eyring_1st_order_linearly_ramped_temperature(): 887 from scipy.special import expi 888 889 def analytic_unit0(t, T0, dH, dS): 890 R = 8.314472 891 kB = 1.3806504e-23 892 h = 6.62606896e-34 893 A = kB / h * np.exp(dS / R) 894 B = dH / R 895 return np.exp( 896 A 897 * ( 898 (-(B ** 2) * np.exp(B / T0) * expi(-B / T0) - T0 * (B - T0)) 899 * np.exp(-B / T0) 900 + ( 901 B ** 2 * np.exp(B / (t + T0)) * expi(-B / (t + T0)) 902 - (t + T0) * (-B + t + T0) 903 ) 904 * np.exp(-B / (t + T0)) 905 ) 906 / 2 907 ) 908 909 T_K = 290 910 dH = 80e3 911 dS = 10 912 rsys1 = ReactionSystem.from_string( 913 """ 914 NOBr -> NO + Br; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol) 915 """.format( 916 dH=dH, dS=dS 917 ) 918 ) 919 920 NOBr0_M = 0.7 921 init_cond = dict(NOBr=NOBr0_M * u.M, NO=0 * u.M, Br=0 * u.M) 922 t = 20 * u.second 923 924 def check(rsys): 925 odes, extra = get_odesys( 926 rsys, 927 unit_registry=SI_base_registry, 928 constants=const, 929 substitutions={"temperature": RampedTemp([T_K * u.K, 1 * u.K / u.s])}, 930 ) 931 for odesys in [odes, odes.as_autonomous()]: 932 res = odesys.integrate(t, init_cond, integrator="cvode") 933 t_sec = to_unitless(res.xout, u.second) 934 NOBr_ref = NOBr0_M * analytic_unit0(t_sec, T_K, dH, dS) 935 cmp = to_unitless(res.yout, u.M) 936 ref = np.empty_like(cmp) 937 ref[:, odesys.names.index("NOBr")] = NOBr_ref 938 ref[:, odesys.names.index("Br")] = NOBr0_M - NOBr_ref 939 ref[:, odesys.names.index("NO")] = NOBr0_M - NOBr_ref 940 assert np.allclose(cmp, ref) 941 942 check(rsys1) 943 rsys2 = ReactionSystem.from_string( 944 """ 945 NOBr -> NO + Br; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol])) 946 """.format( 947 dH=dH, dS=dS 948 ) 949 ) 950 check(rsys2) 951 952 953@requires("pycvodes", "sym", "scipy", units_library) 954def test_get_odesys__Eyring_2nd_order_linearly_ramped_temperature(): 955 from scipy.special import expi 956 957 def analytic_unit0(t, k, m, dH, dS): 958 R = 8.314472 959 kB = 1.3806504e-23 960 h = 6.62606896e-34 961 A = kB / h * np.exp(dS / R) 962 B = dH / R 963 return ( 964 k 965 * np.exp(B * (k * t + 2 * m) / (m * (k * t + m))) 966 / ( 967 A 968 * ( 969 -(B ** 2) * np.exp(B / (k * t + m)) * expi(-B / (k * t + m)) 970 - B * k * t 971 - B * m 972 + k ** 2 * t ** 2 973 + 2 * k * m * t 974 + m ** 2 975 ) 976 * np.exp(B / m) 977 + ( 978 A * B ** 2 * np.exp(B / m) * expi(-B / m) 979 - A * m * (-B + m) 980 + k * np.exp(B / m) 981 ) 982 * np.exp(B / (k * t + m)) 983 ) 984 ) 985 986 T_K = 290 987 dTdt_Ks = 3 988 dH = 80e3 989 dS = 10 990 rsys1 = ReactionSystem.from_string( 991 """ 992 2 NO2 -> N2O4; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol) 993 """.format( 994 dH=dH, dS=dS 995 ) 996 ) 997 998 NO2_M = 1.0 999 init_cond = dict(NO2=NO2_M * u.M, N2O4=0 * u.M) 1000 t = 20 * u.second 1001 1002 def check(rsys): 1003 odes, extra = get_odesys( 1004 rsys, 1005 unit_registry=SI_base_registry, 1006 constants=const, 1007 substitutions={"temperature": RampedTemp([T_K * u.K, dTdt_Ks * u.K / u.s])}, 1008 ) 1009 for odesys in [odes, odes.as_autonomous()]: 1010 res = odesys.integrate(t, init_cond, integrator="cvode") 1011 t_sec = to_unitless(res.xout, u.second) 1012 NO2_ref = analytic_unit0(t_sec, dTdt_Ks, T_K, dH, dS) 1013 cmp = to_unitless(res.yout, u.M) 1014 ref = np.empty_like(cmp) 1015 ref[:, odesys.names.index("NO2")] = NO2_ref 1016 ref[:, odesys.names.index("N2O4")] = (NO2_M - NO2_ref) / 2 1017 assert np.allclose(cmp, ref) 1018 1019 check(rsys1) 1020 rsys2 = ReactionSystem.from_string( 1021 """ 1022 2 NO2 -> N2O4; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol])) 1023 """.format( 1024 dH=dH, dS=dS 1025 ) 1026 ) 1027 check(rsys2) 1028 1029 1030@requires("pycvodes", "sym", units_library) 1031def test_get_odesys__Eyring_2nd_order_reversible(): 1032 R = 8.314472 1033 T_K = 273.15 + 20 # 20 degree celsius 1034 kB = 1.3806504e-23 1035 h = 6.62606896e-34 1036 1037 dHf = 74e3 1038 dSf = R * np.log(h / kB / T_K * 1e16) 1039 1040 dHb = 79e3 1041 dSb = dSf - 23 1042 1043 rsys1 = ReactionSystem.from_string( 1044 """ 1045 Fe+3 + SCN- -> FeSCN+2; EyringParam(dH={dHf}*J/mol, dS={dSf}*J/K/mol) 1046 FeSCN+2 -> Fe+3 + SCN-; EyringParam(dH={dHb}*J/mol, dS={dSb}*J/K/mol) 1047 """.format( 1048 dHf=dHf, dSf=dSf, dHb=dHb, dSb=dSb 1049 ) 1050 ) 1051 kf_ref = 20836643994.118652 * T_K * np.exp(-(dHf - T_K * dSf) / (R * T_K)) 1052 kb_ref = 20836643994.118652 * T_K * np.exp(-(dHb - T_K * dSb) / (R * T_K)) 1053 Fe0 = 6e-3 1054 SCN0 = 2e-3 1055 init_cond = {"Fe+3": Fe0 * u.M, "SCN-": SCN0 * u.M, "FeSCN+2": 0 * u.M} 1056 t = 3 * u.second 1057 1058 def check(rsys, params): 1059 odes, extra = get_odesys( 1060 rsys, include_params=False, unit_registry=SI_base_registry, constants=const 1061 ) 1062 for odesys in [odes, odes.as_autonomous()]: 1063 res = odesys.integrate(t, init_cond, params, integrator="cvode") 1064 t_sec = to_unitless(res.xout, u.second) 1065 FeSCN_ref = binary_rev(t_sec, kf_ref, kb_ref, 0, Fe0, SCN0) 1066 cmp = to_unitless(res.yout, u.M) 1067 ref = np.empty_like(cmp) 1068 ref[:, odesys.names.index("FeSCN+2")] = FeSCN_ref 1069 ref[:, odesys.names.index("Fe+3")] = Fe0 - FeSCN_ref 1070 ref[:, odesys.names.index("SCN-")] = SCN0 - FeSCN_ref 1071 assert np.allclose(cmp, ref) 1072 1073 check(rsys1, {"temperature": T_K * u.K}) 1074 rsys2 = ReactionSystem.from_string( 1075 """ 1076 Fe+3 + SCN- -> FeSCN+2; MassAction(EyringHS([{dHf}*J/mol, {dSf}*J/K/mol])) 1077 FeSCN+2 -> Fe+3 + SCN-; MassAction(EyringHS([{dHb}*J/mol, {dSb}*J/K/mol])) 1078 """.format( 1079 dHf=dHf, dSf=dSf, dHb=dHb, dSb=dSb 1080 ) 1081 ) 1082 check(rsys2, {"temperature": T_K * u.K}) 1083 rsys3 = ReactionSystem.from_string( 1084 """ 1085 Fe+3 + SCN- -> FeSCN+2; MassAction(EyringHS.fk('dHf', 'dSf')) 1086 FeSCN+2 -> Fe+3 + SCN-; MassAction(EyringHS.fk('dHb', 'dSb')) 1087 """ 1088 ) 1089 check( 1090 rsys3, 1091 dict( 1092 temperature=T_K * u.K, 1093 dHf=dHf * u.J / u.mol, 1094 dSf=dSf * u.J / u.mol / u.K, 1095 dHb=dHb * u.J / u.mol, 1096 dSb=dSb * u.J / u.mol / u.K, 1097 ), 1098 ) 1099 1100 1101@requires("numpy", "pyodesys", "sympy", "pycvodes") 1102def test_create_odesys(): 1103 rsys = ReactionSystem.from_string( 1104 """ 1105 A -> B; 'k1' 1106 B + C -> P; 'k2' 1107 """, 1108 substance_factory=Substance, 1109 ) 1110 1111 odesys, odesys_extra = create_odesys(rsys, unit_registry=SI_base_registry) 1112 1113 tend_ul = 10 1114 init_conc_ul = {"A": 1e-6, "B": 0, "C": 1} 1115 params_ul = dict(k1=3, k2=4) 1116 1117 tend = tend_ul * u.s 1118 params = {"k1": params_ul["k1"] / u.s, "k2": params_ul["k2"] / u.M / u.s} 1119 init_conc = {k: v * u.molar for k, v in init_conc_ul.items()} 1120 1121 validation = odesys_extra["validate"](dict(init_conc, **params)) 1122 (P,) = validation["not_seen"] 1123 assert P == "P" 1124 1125 ref_rates = { 1126 "A": -params["k1"] * init_conc["A"], 1127 "P": params["k2"] * init_conc["B"] * init_conc["C"], 1128 } 1129 ref_rates["B"] = -ref_rates["A"] - ref_rates["P"] 1130 ref_rates["C"] = -ref_rates["P"] 1131 assert validation["rates"] == ref_rates 1132 1133 result1, result1_extra = odesys_extra["unit_aware_solve"]( 1134 tend, defaultdict(lambda: 0 * u.molar, init_conc), params, integrator="cvode" 1135 ) 1136 assert result1.info["success"] 1137 1138 result2 = odesys.integrate( 1139 tend_ul, defaultdict(lambda: 0, init_conc_ul), params_ul, integrator="cvode" 1140 ) 1141 assert np.allclose(result2.yout[-1, :], to_unitless(result1.yout[-1, :], u.molar)) 1142 1143 1144@requires("pycvodes", "sym", units_library) 1145def test_create_odesys__Radiolytic(): 1146 rsys1 = ReactionSystem.from_string( 1147 """ 1148 -> e-(aq); Radiolytic.fk('g_emaq') 1149 """, 1150 checks=(), 1151 ) 1152 ic1 = {"e-(aq)": 0.0} 1153 t1 = 5 1154 p1 = dict(g_emaq=42.0, doserate=17.0, density=5.0) 1155 odesys1, odesys_extra = create_odesys(rsys1) 1156 result1 = odesys1.integrate(t1, ic1, p1) 1157 yref1 = result1.xout * p1["g_emaq"] * p1["doserate"] * p1["density"] 1158 assert np.allclose(yref1, result1.yout.squeeze()) 1159 1160 1161@requires("pycvodes", "sym", units_library) 1162def test_create_odesys__validate__catalyst(): 1163 rsys1 = ReactionSystem.from_string( 1164 """ 1165 H2O2 + Pt -> 2 OH + Pt; 'k_decomp' 1166 """ 1167 ) 1168 ic1 = defaultdict(lambda: 0 * u.molar, {"H2O2": 3.0 * u.molar, "Pt": 0.5 * u.molar}) 1169 t1 = linspace(0 * u.s, 0.3 * u.s, 7) 1170 p1 = dict(k_decomp=42 / u.second / u.molar) 1171 odesys1, odesys_extra = create_odesys(rsys1) 1172 validation = odesys_extra["validate"](dict(ic1, **p1)) 1173 assert validation["not_seen"] == {"OH"} 1174 1175 dedim_ctx = _mk_dedim(SI_base_registry) 1176 (t, c, _p), dedim_extra = dedim_ctx["dedim_tcp"]( 1177 t1, [ic1[k] for k in odesys1.names], p1 1178 ) 1179 result1 = odesys1.integrate(t, c, _p) 1180 tout = result1.xout * dedim_extra["unit_time"] 1181 cout = result1.yout * dedim_extra["unit_conc"] 1182 yref1 = ic1["H2O2"] * np.exp(-tout * ic1["Pt"] * p1["k_decomp"]) 1183 assert allclose(yref1, cout[:, odesys1.names.index("H2O2")], rtol=1e-6) 1184 1185 1186@requires("pyodesys", units_library) 1187def test_create_odesys__ShiftedTPoly(): 1188 rxn = Reaction({"A": 1, "B": 1}, {"C": 3, "D": 2}, "k_bi", {"A": 3}) 1189 rsys = ReactionSystem([rxn], "A B C D") 1190 1191 _k0, _k1, T0C = 10, 2, 273.15 1192 rate = MassAction( 1193 ShiftedTPoly([T0C * u.K, _k0 / u.molar / u.s, _k1 / u.molar / u.s / u.K]) 1194 ) 1195 T_C = 25 1196 T = (T0C + T_C) * u.kelvin 1197 p1 = rate.rate_coeff({"temperature": T}) 1198 assert allclose(p1, (_k0 + _k1 * T_C) / u.molar / u.s) 1199 1200 odesys, odesys_extra = create_odesys(rsys) 1201 ics = {"A": 3 * u.molar, "B": 5 * u.molar, "C": 11 * u.molar, "D": 13 * u.molar} 1202 pars = dict(k_bi=p1) 1203 validation = odesys_extra["validate"](dict(ics, **pars)) 1204 assert set(map(str, validation["not_seen"])) == {"C", "D"} 1205 1206 dedim_ctx = _mk_dedim(SI_base_registry) 1207 (t, c, _p), dedim_extra = dedim_ctx["dedim_tcp"]( 1208 -37 * u.s, [ics[k] for k in odesys.names], pars 1209 ) 1210 fout = odesys.f_cb(t, c, [_p[pk] for pk in odesys.param_names]) 1211 r = 3 * 5 * (_k0 + _k1 * 25) * 1000 # mol/m3/s 1212 ref = [-4 * r, -r, 3 * r, 2 * r] 1213 assert np.all(abs((fout - ref) / ref) < 1e-14) 1214 1215 odesys.integrate(t, c, _p) 1216