diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index ddc48bc4d..cbccd05ee 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -4,15 +4,16 @@ This module contains the Operator class, which is then passed to an integrator to run a full depletion simulation. """ -from collections import namedtuple -from collections import defaultdict -from collections.abc import Iterable +from collections import namedtuple, defaultdict +from collections.abc import Iterable, Callable import os from pathlib import Path from abc import ABC, abstractmethod from copy import deepcopy from warnings import warn from numbers import Real, Integral +from inspect import signature +import time from numpy import nonzero, empty, asarray from uncertainties import ufloat @@ -23,6 +24,7 @@ from openmc.checkvalue import check_type, check_greater_than from .results import Results from .chain import Chain from .results_list import ResultsList +from .pool import deplete __all__ = [ @@ -595,7 +597,7 @@ class TalliedFissionYieldHelper(FissionYieldHelper): class Integrator(ABC): - """Abstract class for solving the time-integration for depletion + r"""Abstract class for solving the time-integration for depletion Parameters ---------- @@ -623,6 +625,17 @@ class Integrator(ABC): seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates that the values are given in burnup (MW-d of energy deposited per kilogram of initial heavy metal). + solver : str or callable, optional + If a string, must be the tame of the solver responsible for + solving the Bateman equations. Current options are: + + * ``cram16`` - 16th order IPF CRAM + * ``cram48`` - 48th order IPF CRAM [default] + + If a function or other callable, must adhere to the requirements in + :attr:`solver`. + + .. versionadded:: 0.12 Attributes ---------- @@ -634,10 +647,27 @@ class Integrator(ABC): Size of each depletion interval in [s] power : iterable of float Power of the reactor in [W] for each interval in :attr:`timesteps` + solver : callable + Function that will solve the Bateman equations + :math:`\vec{n}_{i+1} = A_i\vec{n}_i` with a step size :math:`t_i`. Can + be configured using the ``solver`` argument. User-supplied functions + are expected to have the following signature: + ``solver(A, n0, t) -> n1`` where + + * ``A`` is a :class:`scipy.sparse.csr_matrix` making up the + depletion matrix + * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions + for a given material in atoms/cm3 + * ``t`` is a float of the time step size in seconds, and + * ``n1`` is a :class:`numpy.ndarray` of compositions at the + next time step. Expected to be of the same shape as ``n0`` + + .. versionadded:: 0.12 + """ def __init__(self, operator, timesteps, power=None, power_density=None, - timestep_units='s'): + timestep_units='s', solver="cram48"): # Check number of stages previously used if operator.prev_res is not None: res = operator.prev_res[-1] @@ -677,24 +707,24 @@ class Integrator(ABC): # Determine number of seconds for each timestep seconds = [] - for time, unit, watts in zip(times, units, power): + for timestep, unit, watts in zip(times, units, power): # Make sure values passed make sense - check_type('timestep', time, Real) - check_greater_than('timestep', time, 0.0, False) + check_type('timestep', timestep, Real) + check_greater_than('timestep', timestep, 0.0, False) check_type('timestep units', unit, str) check_type('power', watts, Real) check_greater_than('power', watts, 0.0, True) if unit in ('s', 'sec'): - seconds.append(time) + seconds.append(timestep) elif unit in ('min', 'minute'): - seconds.append(time*_SECONDS_PER_MINUTE) + seconds.append(timestep*_SECONDS_PER_MINUTE) elif unit in ('h', 'hr', 'hour'): - seconds.append(time*_SECONDS_PER_HOUR) + seconds.append(timestep*_SECONDS_PER_HOUR) elif unit in ('d', 'day'): - seconds.append(time*_SECONDS_PER_DAY) + seconds.append(timestep*_SECONDS_PER_DAY) elif unit.lower() == 'mwd/kg': - watt_days_per_kg = 1e6*time + watt_days_per_kg = 1e6*timestep kilograms = 1e-3*operator.heavy_metal days = watt_days_per_kg * kilograms / watts seconds.append(days*_SECONDS_PER_DAY) @@ -704,6 +734,59 @@ class Integrator(ABC): self.timesteps = asarray(seconds) self.power = asarray(power) + if isinstance(solver, str): + # Delay importing of cram module, which requires this file + if solver == "cram48": + from .cram import CRAM48 + self._solver = CRAM48 + elif solver == "cram16": + from .cram import CRAM16 + self._solver = CRAM16 + else: + raise ValueError( + "Solver {} not understood. Expected 'cram48' or " + "'cram16'".format(solver)) + else: + self.solver = solver + + @property + def solver(self): + return self._solver + + @solver.setter + def solver(self, func): + if not isinstance(func, Callable): + raise TypeError( + "Solver must be callable, not {}".format(type(func))) + try: + sig = signature(func) + except ValueError: + # Guard against callables that aren't introspectable, e.g. + # fortran functions wrapped by F2PY + warn("Could not determine arguments to {}. Proceeding " + "anyways".format(func)) + self._solver = func + return + + # Inspect arguments + if len(sig.parameters) != 3: + raise ValueError("Function {} does not support three arguments: " + "{!s}".format(func, sig)) + + for ix, param in enumerate(sig.parameters.values()): + if param.kind in {param.KEYWORD_ONLY, param.VAR_KEYWORD}: + raise ValueError( + "Keyword arguments like {} at position {} are not " + "allowed".format(ix, param)) + + self._solver = func + + def _timed_deplete(self, concs, rates, dt, matrix_func=None): + start = time.time() + results = deplete( + self._solver, self.chain, concs, rates, dt, matrix_func) + return time.time() - start, results + @abstractmethod def __call__(self, conc, rates, dt, power, i): """Perform the integration across one time step @@ -807,7 +890,7 @@ class Integrator(ABC): class SIIntegrator(Integrator): - """Abstract class for the Stochastic Implicit Euler integrators + r"""Abstract class for the Stochastic Implicit Euler integrators Does not provide a ``__call__`` method, but scales and resets the number of particles used in initial transport calculation @@ -841,6 +924,17 @@ class SIIntegrator(Integrator): n_steps : int, optional Number of stochastic iterations per depletion interval. Must be greater than zero. Default : 10 + solver : str or callable, optional + If a string, must be the tame of the solver responsible for + solving the Bateman equations. Current options are: + + * ``cram16`` - 16th order IPF CRAM + * ``cram48`` - 48th order IPF CRAM [default] + + If a function or other callable, must adhere to the requirements in + :attr:`solver`. + + .. versionadded:: 0.12 Attributes ---------- @@ -854,12 +948,31 @@ class SIIntegrator(Integrator): Power of the reactor in [W] for each interval in :attr:`timesteps` n_steps : int Number of stochastic iterations per depletion interval + solver : callable + Function that will solve the Bateman equations + :math:`\vec{n}_{i+1} = A_i\vec{n}_i` with a step size :math:`t_i`. Can + be configured using the ``solver`` argument. User-supplied functions + are expected to have the following signature: + ``solver(A, n0, t) -> n1`` where + + * ``A`` is a :class:`scipy.sparse.csr_matrix` making up the + depletion matrix + * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions + for a given material in atoms/cm3 + * ``t`` is a float of the time step size in seconds, and + * ``n1`` is a :class:`numpy.ndarray` of compositions at the + next time step. Expected to be of the same shape as ``n0`` + + .. versionadded:: 0.12 + """ def __init__(self, operator, timesteps, power=None, power_density=None, - timestep_units='s', n_steps=10): + timestep_units='s', n_steps=10, solver="cram48"): check_type("n_steps", n_steps, Integral) check_greater_than("n_steps", n_steps, 0) - super().__init__(operator, timesteps, power, power_density, timestep_units) + super().__init__( + operator, timesteps, power, power_density, timestep_units, + solver=solver) self.n_steps = n_steps def _get_bos_data_from_operator(self, step_index, step_power, bos_conc): diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index ad0137221..020291555 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -2,7 +2,6 @@ import copy from itertools import repeat from .abc import Integrator, SIIntegrator, OperatorResult -from .cram import CRAM48, timed_deplete from ._matrix_funcs import ( cf4_f1, cf4_f2, cf4_f3, cf4_f4, celi_f1, celi_f2, leqi_f1, leqi_f2, leqi_f3, leqi_f4, rk4_f1, rk4_f4 @@ -54,6 +53,17 @@ class PredictorIntegrator(Integrator): that the values are given in burnup (MW-d of energy deposited per kilogram of initial heavy metal). + .. versionadded:: 0.12 + solver : str or callable, optional + If a string, must be the tame of the solver responsible for + solving the Bateman equations. Current options are: + + * ``cram16`` - 16th order IPF CRAM + * ``cram48`` - 48th order IPF CRAM [default] + + If a function or other callable, must adhere to the requirements in + :attr:`solver`. + .. versionadded:: 0.12 Attributes @@ -66,6 +76,23 @@ class PredictorIntegrator(Integrator): Size of each depletion interval in [s] power : iterable of float Power of the reactor in [W] for each interval in :attr:`timesteps` + solver : callable + Function that will solve the Bateman equations + :math:`\vec{n}_{i+1} = A_i\vec{n}_i` with a step size :math:`t_i`. Can + be configured using the ``solver`` argument. User-supplied functions + are expected to have the following signature: + ``solver(A, n0, t) -> n1`` where + + * ``A`` is a :class:`scipy.sparse.csr_matrix` making up the + depletion matrix + * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions + for a given material in atoms/cm3 + * ``t`` is a float of the time step size in seconds, and + * ``n1`` is a :class:`numpy.ndarray` of compositions at the + next time step. Expected to be of the same shape as ``n0`` + + .. versionadded:: 0.12 + """ _num_stages = 1 @@ -96,8 +123,7 @@ class PredictorIntegrator(Integrator): operator with predictor """ - proc_time, conc_end = timed_deplete( - CRAM48, self.chain, conc, rates, dt) + proc_time, conc_end = self._timed_deplete(conc, rates, dt) return proc_time, [conc_end], [] @@ -146,6 +172,17 @@ class CECMIntegrator(Integrator): that the values are given in burnup (MW-d of energy deposited per kilogram of initial heavy metal). + .. versionadded:: 0.12 + solver : str or callable, optional + If a string, must be the tame of the solver responsible for + solving the Bateman equations. Current options are: + + * ``cram16`` - 16th order IPF CRAM + * ``cram48`` - 48th order IPF CRAM [default] + + If a function or other callable, must adhere to the requirements in + :attr:`solver`. + .. versionadded:: 0.12 Attributes @@ -158,6 +195,23 @@ class CECMIntegrator(Integrator): Size of each depletion interval in [s] power : iterable of float Power of the reactor in [W] for each interval in :attr:`timesteps` + solver : callable + Function that will solve the Bateman equations + :math:`\vec{n}_{i+1} = A_i\vec{n}_i` with a step size :math:`t_i`. Can + be configured using the ``solver`` argument. User-supplied functions + are expected to have the following signature: + ``solver(A, n0, t) -> n1`` where + + * ``A`` is a :class:`scipy.sparse.csr_matrix` making up the + depletion matrix + * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions + for a given material in atoms/cm3 + * ``t`` is a float of the time step size in seconds, and + * ``n1`` is a :class:`numpy.ndarray` of compositions at the + next time step. Expected to be of the same shape as ``n0`` + + .. versionadded:: 0.12 + """ _num_stages = 2 @@ -188,14 +242,12 @@ class CECMIntegrator(Integrator): Eigenvalue and reaction rates from transport simulations """ # deplete across first half of inteval - time0, x_middle = timed_deplete( - CRAM48, self.chain, conc, rates, dt / 2) + time0, x_middle = self._timed_deplete(conc, rates, dt / 2) res_middle = self.operator(x_middle, power) # deplete across entire interval with BOS concentrations, # MOS reaction rates - time1, x_end = timed_deplete( - CRAM48, self.chain, conc, res_middle.rates, dt) + time1, x_end = self._timed_deplete(conc, res_middle.rates, dt) return time0 + time1, [x_middle, x_end], [res_middle] @@ -247,6 +299,17 @@ class CF4Integrator(Integrator): that the values are given in burnup (MW-d of energy deposited per kilogram of initial heavy metal). + .. versionadded:: 0.12 + solver : str or callable, optional + If a string, must be the tame of the solver responsible for + solving the Bateman equations. Current options are: + + * ``cram16`` - 16th order IPF CRAM + * ``cram48`` - 48th order IPF CRAM [default] + + If a function or other callable, must adhere to the requirements in + :attr:`solver`. + .. versionadded:: 0.12 Attributes @@ -259,6 +322,23 @@ class CF4Integrator(Integrator): Size of each depletion interval in [s] power : iterable of float Power of the reactor in [W] for each interval in :attr:`timesteps` + solver : callable + Function that will solve the Bateman equations + :math:`\vec{n}_{i+1} = A_i\vec{n}_i` with a step size :math:`t_i`. Can + be configured using the ``solver`` argument. User-supplied functions + are expected to have the following signature: + ``solver(A, n0, t) -> n1`` where + + * ``A`` is a :class:`scipy.sparse.csr_matrix` making up the + depletion matrix + * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions + for a given material in atoms/cm3 + * ``t`` is a float of the time step size in seconds, and + * ``n1`` is a :class:`numpy.ndarray` of compositions at the + next time step. Expected to be of the same shape as ``n0`` + + .. versionadded:: 0.12 + """ _num_stages = 4 @@ -290,27 +370,27 @@ class CF4Integrator(Integrator): simulations """ # Step 1: deplete with matrix 1/2*A(y0) - time1, conc_eos1 = timed_deplete( - CRAM48, self.chain, bos_conc, bos_rates, dt, matrix_func=cf4_f1) + time1, conc_eos1 = self._timed_deplete( + bos_conc, bos_rates, dt, matrix_func=cf4_f1) res1 = self.operator(conc_eos1, power) # Step 2: deplete with matrix 1/2*A(y1) - time2, conc_eos2 = timed_deplete( - CRAM48, self.chain, bos_conc, res1.rates, dt, matrix_func=cf4_f1) + time2, conc_eos2 = self._timed_deplete( + bos_conc, res1.rates, dt, matrix_func=cf4_f1) res2 = self.operator(conc_eos2, power) # Step 3: deplete with matrix -1/2*A(y0)+A(y2) list_rates = list(zip(bos_rates, res2.rates)) - time3, conc_eos3 = timed_deplete( - CRAM48, self.chain, conc_eos1, list_rates, dt, matrix_func=cf4_f2) + time3, conc_eos3 = self._timed_deplete( + conc_eos1, list_rates, dt, matrix_func=cf4_f2) res3 = self.operator(conc_eos3, power) # Step 4: deplete with two matrix exponentials list_rates = list(zip(bos_rates, res1.rates, res2.rates, res3.rates)) - time4, conc_inter = timed_deplete( - CRAM48, self.chain, bos_conc, list_rates, dt, matrix_func=cf4_f3) - time5, conc_eos5 = timed_deplete( - CRAM48, self.chain, conc_inter, list_rates, dt, matrix_func=cf4_f4) + time4, conc_inter = self._timed_deplete( + bos_conc, list_rates, dt, matrix_func=cf4_f3) + time5, conc_eos5 = self._timed_deplete( + conc_inter, list_rates, dt, matrix_func=cf4_f4) return (time1 + time2 + time3 + time4 + time5, [conc_eos1, conc_eos2, conc_eos3, conc_eos5], @@ -363,6 +443,17 @@ class CELIIntegrator(Integrator): that the values are given in burnup (MW-d of energy deposited per kilogram of initial heavy metal). + .. versionadded:: 0.12 + solver : str or callable, optional + If a string, must be the tame of the solver responsible for + solving the Bateman equations. Current options are: + + * ``cram16`` - 16th order IPF CRAM + * ``cram48`` - 48th order IPF CRAM [default] + + If a function or other callable, must adhere to the requirements in + :attr:`solver`. + .. versionadded:: 0.12 Attributes @@ -375,6 +466,23 @@ class CELIIntegrator(Integrator): Size of each depletion interval in [s] power : iterable of float Power of the reactor in [W] for each interval in :attr:`timesteps` + solver : callable + Function that will solve the Bateman equations + :math:`\vec{n}_{i+1} = A_i\vec{n}_i` with a step size :math:`t_i`. Can + be configured using the ``solver`` argument. User-supplied functions + are expected to have the following signature: + ``solver(A, n0, t) -> n1`` where + + * ``A`` is a :class:`scipy.sparse.csr_matrix` making up the + depletion matrix + * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions + for a given material in atoms/cm3 + * ``t`` is a float of the time step size in seconds, and + * ``n1`` is a :class:`numpy.ndarray` of compositions at the + next time step. Expected to be of the same shape as ``n0`` + + .. versionadded:: 0.12 + """ _num_stages = 2 @@ -406,18 +514,17 @@ class CELIIntegrator(Integrator): simulation """ # deplete to end using BOS rates - proc_time, conc_ce = timed_deplete( - CRAM48, self.chain, bos_conc, rates, dt) + proc_time, conc_ce = self._timed_deplete(bos_conc, rates, dt) res_ce = self.operator(conc_ce, power) # deplete using two matrix exponentials list_rates = list(zip(rates, res_ce.rates)) - time_le1, conc_inter = timed_deplete( - CRAM48, self.chain, bos_conc, list_rates, dt, matrix_func=celi_f1) + time_le1, conc_inter = self._timed_deplete( + bos_conc, list_rates, dt, matrix_func=celi_f1) - time_le2, conc_end = timed_deplete( - CRAM48, self.chain, conc_inter, list_rates, dt, matrix_func=celi_f2) + time_le2, conc_end = self._timed_deplete( + conc_inter, list_rates, dt, matrix_func=celi_f2) return proc_time + time_le1 + time_le1, [conc_ce, conc_end], [res_ce] @@ -467,6 +574,17 @@ class EPCRK4Integrator(Integrator): that the values are given in burnup (MW-d of energy deposited per kilogram of initial heavy metal). + .. versionadded:: 0.12 + solver : str or callable, optional + If a string, must be the tame of the solver responsible for + solving the Bateman equations. Current options are: + + * ``cram16`` - 16th order IPF CRAM + * ``cram48`` - 48th order IPF CRAM [default] + + If a function or other callable, must adhere to the requirements in + :attr:`solver`. + .. versionadded:: 0.12 Attributes @@ -479,6 +597,18 @@ class EPCRK4Integrator(Integrator): Size of each depletion interval in [s] power : iterable of float Power of the reactor in [W] for each interval in :attr:`timesteps` + solver : str or callable, optional + If a string, must be the tame of the solver responsible for + solving the Bateman equations. Current options are: + + * ``cram16`` - 16th order IPF CRAM + * ``cram48`` - 48th order IPF CRAM [default] + + If a function or other callable, must adhere to the requirements in + :attr:`solver`. + + .. versionadded:: 0.12 + """ _num_stages = 4 @@ -511,24 +641,23 @@ class EPCRK4Integrator(Integrator): """ # Step 1: deplete with matrix A(y0) / 2 - time1, conc1 = timed_deplete( - CRAM48, self.chain, conc, rates, dt, matrix_func=rk4_f1) + time1, conc1 = self._timed_deplete( + conc, rates, dt, matrix_func=rk4_f1) res1 = self.operator(conc1, power) # Step 2: deplete with matrix A(y1) / 2 - time2, conc2 = timed_deplete( - CRAM48, self.chain, conc, res1.rates, dt, matrix_func=rk4_f1) + time2, conc2 = self._timed_deplete( + conc, res1.rates, dt, matrix_func=rk4_f1) res2 = self.operator(conc2, power) # Step 3: deplete with matrix A(y2) - time3, conc3 = timed_deplete( - CRAM48, self.chain, conc, res2.rates, dt) + time3, conc3 = self._timed_deplete(conc, res2.rates, dt) res3 = self.operator(conc3, power) # Step 4: deplete with matrix built from weighted rates list_rates = list(zip(rates, res1.rates, res2.rates, res3.rates)) - time4, conc4 = timed_deplete( - CRAM48, self.chain, conc, list_rates, dt, matrix_func=rk4_f4) + time4, conc4 = self._timed_deplete( + conc, list_rates, dt, matrix_func=rk4_f4) return (time1 + time2 + time3 + time4, [conc1, conc2, conc3, conc4], [res1, res2, res3]) @@ -590,6 +719,17 @@ class LEQIIntegrator(Integrator): that the values are given in burnup (MW-d of energy deposited per kilogram of initial heavy metal). + .. versionadded:: 0.12 + solver : str or callable, optional + If a string, must be the tame of the solver responsible for + solving the Bateman equations. Current options are: + + * ``cram16`` - 16th order IPF CRAM + * ``cram48`` - 48th order IPF CRAM [default] + + If a function or other callable, must adhere to the requirements in + :attr:`solver`. + .. versionadded:: 0.12 Attributes @@ -602,6 +742,23 @@ class LEQIIntegrator(Integrator): Size of each depletion interval in [s] power : iterable of float Power of the reactor in [W] for each interval in :attr:`timesteps` + solver : callable + Function that will solve the Bateman equations + :math:`\vec{n}_{i+1} = A_i\vec{n}_i` with a step size :math:`t_i`. Can + be configured using the ``solver`` argument. User-supplied functions + are expected to have the following signature: + ``solver(A, n0, t) -> n1`` where + + * ``A`` is a :class:`scipy.sparse.csr_matrix` making up the + depletion matrix + * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions + for a given material in atoms/cm3 + * ``t`` is a float of the time step size in seconds, and + * ``n1`` is a :class:`numpy.ndarray` of compositions at the + next time step. Expected to be of the same shape as ``n0`` + + .. versionadded:: 0.12 + """ _num_stages = 2 @@ -649,10 +806,10 @@ class LEQIIntegrator(Integrator): le_inputs = list(zip( self._prev_rates, bos_res.rates, repeat(prev_dt), repeat(dt))) - time1, conc_inter = timed_deplete( - CRAM48, self.chain, bos_conc, le_inputs, dt, matrix_func=leqi_f1) - time2, conc_eos0 = timed_deplete( - CRAM48, self.chain, conc_inter, le_inputs, dt, matrix_func=leqi_f2) + time1, conc_inter = self._timed_deplete( + bos_conc, le_inputs, dt, matrix_func=leqi_f1) + time2, conc_eos0 = self._timed_deplete( + conc_inter, le_inputs, dt, matrix_func=leqi_f2) res_inter = self.operator(conc_eos0, power) @@ -660,10 +817,10 @@ class LEQIIntegrator(Integrator): self._prev_rates, bos_res.rates, res_inter.rates, repeat(prev_dt), repeat(dt))) - time3, conc_inter = timed_deplete( - CRAM48, self.chain, bos_conc, qi_inputs, dt, matrix_func=leqi_f3) - time4, conc_eos1 = timed_deplete( - CRAM48, self.chain, conc_inter, qi_inputs, dt, matrix_func=leqi_f4) + time3, conc_inter = self._timed_deplete( + bos_conc, qi_inputs, dt, matrix_func=leqi_f3) + time4, conc_eos1 = self._timed_deplete( + conc_inter, qi_inputs, dt, matrix_func=leqi_f4) # store updated rates self._prev_rates = copy.deepcopy(bos_res.rates) @@ -714,6 +871,17 @@ class SICELIIntegrator(SIIntegrator): n_steps : int, optional Number of stochastic iterations per depletion interval. Must be greater than zero. Default : 10 + solver : str or callable, optional + If a string, must be the tame of the solver responsible for + solving the Bateman equations. Current options are: + + * ``cram16`` - 16th order IPF CRAM + * ``cram48`` - 48th order IPF CRAM [default] + + If a function or other callable, must adhere to the requirements in + :attr:`solver`. + + .. versionadded:: 0.12 Attributes ---------- @@ -727,6 +895,23 @@ class SICELIIntegrator(SIIntegrator): Power of the reactor in [W] for each interval in :attr:`timesteps` n_steps : int Number of stochastic iterations per depletion interval + solver : callable + Function that will solve the Bateman equations + :math:`\vec{n}_{i+1} = A_i\vec{n}_i` with a step size :math:`t_i`. Can + be configured using the ``solver`` argument. User-supplied functions + are expected to have the following signature: + ``solver(A, n0, t) -> n1`` where + + * ``A`` is a :class:`scipy.sparse.csr_matrix` making up the + depletion matrix + * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions + for a given material in atoms/cm3 + * ``t`` is a float of the time step size in seconds, and + * ``n1`` is a :class:`numpy.ndarray` of compositions at the + next time step. Expected to be of the same shape as ``n0`` + + .. versionadded:: 0.12 + """ _num_stages = 2 @@ -757,8 +942,7 @@ class SICELIIntegrator(SIIntegrator): Eigenvalue and reaction rates from intermediate transport simulations """ - proc_time, eos_conc = timed_deplete( - CRAM48, self.chain, bos_conc, bos_rates, dt) + proc_time, eos_conc = self._timed_deplete(bos_conc, bos_rates, dt) inter_conc = copy.deepcopy(eos_conc) # Begin iteration @@ -773,10 +957,10 @@ class SICELIIntegrator(SIIntegrator): res_bar = OperatorResult(k, rates) list_rates = list(zip(bos_rates, res_bar.rates)) - time1, inter_conc = timed_deplete( - CRAM48, self.chain, bos_conc, list_rates, dt, matrix_func=celi_f1) - time2, inter_conc = timed_deplete( - CRAM48, self.chain, inter_conc, list_rates, dt, matrix_func=celi_f2) + time1, inter_conc = self._timed_deplete( + bos_conc, list_rates, dt, matrix_func=celi_f1) + time2, inter_conc = self._timed_deplete( + inter_conc, list_rates, dt, matrix_func=celi_f2) proc_time += time1 + time2 # end iteration @@ -824,6 +1008,17 @@ class SILEQIIntegrator(SIIntegrator): n_steps : int, optional Number of stochastic iterations per depletion interval. Must be greater than zero. Default : 10 + solver : str or callable, optional + If a string, must be the tame of the solver responsible for + solving the Bateman equations. Current options are: + + * ``cram16`` - 16th order IPF CRAM + * ``cram48`` - 48th order IPF CRAM [default] + + If a function or other callable, must adhere to the requirements in + :attr:`solver`. + + .. versionadded:: 0.12 Attributes ---------- @@ -837,6 +1032,18 @@ class SILEQIIntegrator(SIIntegrator): Power of the reactor in [W] for each interval in :attr:`timesteps` n_steps : int Number of stochastic iterations per depletion interval + solver : str or callable, optional + If a string, must be the tame of the solver responsible for + solving the Bateman equations. Current options are: + + * ``cram16`` - 16th order IPF CRAM + * ``cram48`` - 48th order IPF CRAM [default] + + If a function or other callable, must adhere to the requirements in + :attr:`solver`. + + .. versionadded:: 0.12 + """ _num_stages = 2 @@ -883,10 +1090,10 @@ class SILEQIIntegrator(SIIntegrator): # Perform remaining LE/QI inputs = list(zip(self._prev_rates, bos_rates, repeat(prev_dt), repeat(dt))) - proc_time, inter_conc = timed_deplete( - CRAM48, self.chain, bos_conc, inputs, dt, matrix_func=leqi_f1) - time1, eos_conc = timed_deplete( - CRAM48, self.chain, inter_conc, inputs, dt, matrix_func=leqi_f2) + proc_time, inter_conc = self._timed_deplete( + bos_conc, inputs, dt, matrix_func=leqi_f1) + time1, eos_conc = self._timed_deplete( + inter_conc, inputs, dt, matrix_func=leqi_f2) proc_time += time1 inter_conc = copy.deepcopy(eos_conc) @@ -903,10 +1110,10 @@ class SILEQIIntegrator(SIIntegrator): inputs = list(zip(self._prev_rates, bos_rates, res_bar.rates, repeat(prev_dt), repeat(dt))) - time1, inter_conc = timed_deplete( - CRAM48, self.chain, bos_conc, inputs, dt, matrix_func=leqi_f3) - time2, inter_conc = timed_deplete( - CRAM48, self.chain, inter_conc, inputs, dt, matrix_func=leqi_f4) + time1, inter_conc = self._timed_deplete( + bos_conc, inputs, dt, matrix_func=leqi_f3) + time2, inter_conc = self._timed_deplete( + inter_conc, inputs, dt, matrix_func=leqi_f4) proc_time += time1 + time2 return proc_time, [eos_conc, inter_conc], [res_bar]