Coverage for chemreac/integrate : 88%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
# -*- coding: utf-8 -*- chemreac.integrate ==================
This module provides functions for integrating the system of ODEs which the ReactionDiffusion represent. The main class representing a numerical integration of the system of ODEs is :py:class:`Integration`.
If one does not want to hard code the choice of solver and solver parameters (e.g. tolerances), one may use :py:func:`run` which defers those choices to the user of the script through the use of environment variables. """
print_function, unicode_literals)
'atol': 1e-9, 'rtol': 1e-7, }
""" see integrate.
kwargs: method: linear multistep method: 'bdf' or 'adams'
"""
# Handle kwargs raise NotImplementedError( "Sundials integrator auto-selects banded for N>1") raise KeyError("Unkown kwargs: {}".format(kwargs))
# Run the integration np.asarray(tout).flatten(), **new_kwargs) except RuntimeError: yout = np.empty((len(tout), rd.N, rd.n), order='C')/0 # NaN success = False else:
'neval_f': rd.neval_f, 'neval_j': rd.neval_j, 'texec': texec, 'success': success }) info['nprec_setup'] = rd.nprec_setup info['nprec_solve'] = rd.nprec_solve info['njacvec_dot'] = rd.njacvec_dot
""" For demonstration purposes only, fixed step size give no error control and requires excessive work for accurate solutions. Unknown kwargs are simply ignored.
see integrate """ 'neval_f': 4*(tout.size-1), 'neval_j': 0, 'texec': texec, 'success': True, }
integrator_name='vode', dense_output=False, **kwargs): """ see integrate
Parameters ---------- tout: array-like at what times to report, e.g.: - np.linspace(t0, tend, nt+1) - np.logspace(np.log10(t0 + 1e-12), np.log10(tend), nt+1) integrator_name: string (default: vode) dense_output: bool (default: False) if True, tout is taken to be length 2 tuple (t0, tend)
Returns ======= yout: numpy array of shape (len(tout), rd.N, rd.n)
"""
fmtstr = "y0.size (={})not compatible with rd.n*rd.N (={})" raise ValueError(fmtstr.format(y0.size, rd.n*rd.N))
else: raise NotImplementedError
raise KeyError("Unkown kwargs: {}".format(kwargs))
# Create python callbacks with right signature
# Python function closure circumvents reallocation
# Currently SciPy <= v0.14 needs extra padding jout = rd.alloc_jout(banded=True, order='F', pad=rd.n) else: # SciPy >= v0.15 need no extra padding else: raise NotImplementedError
else:
raise ValueError("dense_output implies tout == (t0, tend)") # suppress warning printed by Fortran else:
'integrator_name': integrator_name, 'success': runner.successful(), 'texec': texec, 'neval_f': f.neval, 'neval_j': jac.neval, })
""" Algebraic sigmoid to avoid overflow/underflow of 'double exp(double)' """
""" Model kinetcs by integrating system of ODEs using user specified solver.
Parameters ---------- solver: string "sundials" or "scipy" where scipy uses VODE as the solver. rd: ReactionDiffusion instance C0: array initial concentrations (untransformed, i.e. linear) tout: array times for which to report solver results (untransformed) sigm_damp: bool or tuple of (lim: float, n: int) conditionally damp C0 with an algebraic sigmoid when rd.logy == True. s(x) = x/((x/lim)**n+1)**(1./n) if sigm==True then `lim` and `n` are the default of `sigm()` C0_is_log: bool If True: passed values in C0 are taken to be the natural logarithm of initial concentrations. If False and rd.logy == True: a very small number is added to C0 to avoid applying log to zero (see `tiny`). tiny: float added to C0 when rd.logy==True and C0_is_log==False. Note that if you explicitly want to avoid adding tiny you need to set it to zero (e.g. when manually setting any C0==0 to some epsilon). (default: None => numpy.finfo(np.float64).tiny)
**kwargs: mode: not supported by Sundials solver (current wrapper code auto selects banded for N>1 and uses dense mode for N==1) atol: float or sequence absolute tolerance of solution rtol: float relative tolerance of solution
Attributes ---------- Cout: array linear output concentrations yout: array output from solver: log(concentrations) if rd.logy == True info: dict Information from solver. Guaranteed to contain: - 'texec': execution time in seconds. - 'atol': float or array, absolute tolerance(s). - 'rtol': float, relative tolerance
Methods ------- _integrate() performs the integration, automatically called by __init__
"""
'sundials': _integrate_sundials, 'scipy': _integrate_scipy, 'rk4': _integrate_rk4, }
C0_is_log=False, tiny=None, **kwargs): raise KeyError("Unknown solver %s" % solver) rd.units, 'concentration')).flatten() rd.units, 'time'))
raise ValueError("Negative concentrations encountered in C0")
""" Performs the integration by calling the callback chosen by self.solver. If rd.logy == True, a transformation of self.C0 to log(C0) will be performed before running the integration (the same is done for self.tout / rd.logt == True).
After the integration is done the attributes `Cout`, `info` and `yout` are set. Cout is guaranteed to be linear concentrations (transformed from yout by calling exp if rd.logy==True) and yout is the unprocessed output from the solver. """
# Transform initial concentrations
y0 = sigm(C0, *self.sigm_damp) else: else: y0 = np.exp(sigm(C0)) y0 = np.exp(sigm(C0, *self.sigm_damp)) else: else:
# Transform time else:
# Run the integration self.rd, y0, t, **self.kwargs)
# Back-transform independent variable into linear time np.exp(self.internal_t) if self.rd.logt else self.internal_t, 'time')
# Back-transform integration output into linear concentration np.exp(self.yout) if self.rd.logy else self.yout, 'concentration')
""" ``run`` is provided for environment variable directed solver choice.
Set ``CHEMREAC_SOLVER`` to indicate what integrator to use (default: "scipy").
Set ``CHEMREAC_SOLVER_KWARGS`` to a string which can be eval'd to a python dictionary. e.g. "{'atol': 1e-4, 'rtol'=1e-7}" """ fmtstr = "CHEMREAC_SOLVER_KWARGS not evaluated to a dictinary: {}" raise TypeError(fmtstr.format(environ_kwargs)) os.getenv('CHEMREAC_SOLVER', 'scipy'), *args, **kwargs) |