From 24212ae731d968380fe80a6fa90121bf1b94d4ff Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 9 May 2020 10:26:14 -0400 Subject: [PATCH 1/6] Require function as first argument to cram.deplete Done in order to remove hardcoded used of CRAM48 --- openmc/deplete/cram.py | 7 +++- openmc/deplete/integrators.py | 58 ++++++++++++++------------ tests/unit_tests/test_deplete_chain.py | 2 +- 3 files changed, 37 insertions(+), 30 deletions(-) diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py index 7dad820a9c..939ff5432c 100644 --- a/openmc/deplete/cram.py +++ b/openmc/deplete/cram.py @@ -20,11 +20,14 @@ __all__ = [ "Cram16Solver", "Cram48Solver", "IPFCramSolver"] -def deplete(chain, x, rates, dt, matrix_func=None): +def deplete(func, chain, x, rates, dt, matrix_func=None): """Deplete materials using given reaction rates for a specified time Parameters ---------- + func : callable + Function to use to get new compositions. Expected to have the + signature ``func(A, n0, t) -> n1`` chain : openmc.deplete.Chain Depletion chain x : list of numpy.ndarray @@ -63,7 +66,7 @@ def deplete(chain, x, rates, dt, matrix_func=None): # Use multiprocessing pool to distribute work with Pool() as pool: inputs = zip(matrices, x, repeat(dt)) - x_result = list(pool.starmap(CRAM48, inputs)) + x_result = list(pool.starmap(func, inputs)) return x_result diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index ddb2130379..ad0137221b 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -2,7 +2,7 @@ import copy from itertools import repeat from .abc import Integrator, SIIntegrator, OperatorResult -from .cram import timed_deplete +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 @@ -96,7 +96,8 @@ class PredictorIntegrator(Integrator): operator with predictor """ - proc_time, conc_end = timed_deplete(self.chain, conc, rates, dt) + proc_time, conc_end = timed_deplete( + CRAM48, self.chain, conc, rates, dt) return proc_time, [conc_end], [] @@ -187,12 +188,14 @@ class CECMIntegrator(Integrator): Eigenvalue and reaction rates from transport simulations """ # deplete across first half of inteval - time0, x_middle = timed_deplete(self.chain, conc, rates, dt / 2) + time0, x_middle = timed_deplete( + CRAM48, self.chain, 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(self.chain, conc, res_middle.rates, dt) + time1, x_end = timed_deplete( + CRAM48, self.chain, conc, res_middle.rates, dt) return time0 + time1, [x_middle, x_end], [res_middle] @@ -288,26 +291,26 @@ class CF4Integrator(Integrator): """ # Step 1: deplete with matrix 1/2*A(y0) time1, conc_eos1 = timed_deplete( - self.chain, bos_conc, bos_rates, dt, matrix_func=cf4_f1) + CRAM48, self.chain, 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( - self.chain, bos_conc, res1.rates, dt, matrix_func=cf4_f1) + CRAM48, self.chain, 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( - self.chain, conc_eos1, list_rates, dt, matrix_func=cf4_f2) + CRAM48, self.chain, 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( - self.chain, bos_conc, list_rates, dt, matrix_func=cf4_f3) + CRAM48, self.chain, bos_conc, list_rates, dt, matrix_func=cf4_f3) time5, conc_eos5 = timed_deplete( - self.chain, conc_inter, list_rates, dt, matrix_func=cf4_f4) + CRAM48, self.chain, conc_inter, list_rates, dt, matrix_func=cf4_f4) return (time1 + time2 + time3 + time4 + time5, [conc_eos1, conc_eos2, conc_eos3, conc_eos5], @@ -403,17 +406,18 @@ class CELIIntegrator(Integrator): simulation """ # deplete to end using BOS rates - proc_time, conc_ce = timed_deplete(self.chain, bos_conc, rates, dt) + proc_time, conc_ce = timed_deplete( + CRAM48, self.chain, 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( - self.chain, bos_conc, list_rates, dt, matrix_func=celi_f1) + CRAM48, self.chain, bos_conc, list_rates, dt, matrix_func=celi_f1) time_le2, conc_end = timed_deplete( - self.chain, conc_inter, list_rates, dt, matrix_func=celi_f2) + CRAM48, self.chain, conc_inter, list_rates, dt, matrix_func=celi_f2) return proc_time + time_le1 + time_le1, [conc_ce, conc_end], [res_ce] @@ -508,23 +512,23 @@ class EPCRK4Integrator(Integrator): # Step 1: deplete with matrix A(y0) / 2 time1, conc1 = timed_deplete( - self.chain, conc, rates, dt, matrix_func=rk4_f1) + CRAM48, self.chain, conc, rates, dt, matrix_func=rk4_f1) res1 = self.operator(conc1, power) # Step 2: deplete with matrix A(y1) / 2 time2, conc2 = timed_deplete( - self.chain, conc, res1.rates, dt, matrix_func=rk4_f1) + CRAM48, self.chain, conc, res1.rates, dt, matrix_func=rk4_f1) res2 = self.operator(conc2, power) # Step 3: deplete with matrix A(y2) time3, conc3 = timed_deplete( - self.chain, conc, res2.rates, dt) + CRAM48, self.chain, 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( - self.chain, conc, list_rates, dt, matrix_func=rk4_f4) + CRAM48, self.chain, conc, list_rates, dt, matrix_func=rk4_f4) return (time1 + time2 + time3 + time4, [conc1, conc2, conc3, conc4], [res1, res2, res3]) @@ -646,9 +650,9 @@ class LEQIIntegrator(Integrator): self._prev_rates, bos_res.rates, repeat(prev_dt), repeat(dt))) time1, conc_inter = timed_deplete( - self.chain, bos_conc, le_inputs, dt, matrix_func=leqi_f1) + CRAM48, self.chain, bos_conc, le_inputs, dt, matrix_func=leqi_f1) time2, conc_eos0 = timed_deplete( - self.chain, conc_inter, le_inputs, dt, matrix_func=leqi_f2) + CRAM48, self.chain, conc_inter, le_inputs, dt, matrix_func=leqi_f2) res_inter = self.operator(conc_eos0, power) @@ -657,9 +661,9 @@ class LEQIIntegrator(Integrator): repeat(prev_dt), repeat(dt))) time3, conc_inter = timed_deplete( - self.chain, bos_conc, qi_inputs, dt, matrix_func=leqi_f3) + CRAM48, self.chain, bos_conc, qi_inputs, dt, matrix_func=leqi_f3) time4, conc_eos1 = timed_deplete( - self.chain, conc_inter, qi_inputs, dt, matrix_func=leqi_f4) + CRAM48, self.chain, conc_inter, qi_inputs, dt, matrix_func=leqi_f4) # store updated rates self._prev_rates = copy.deepcopy(bos_res.rates) @@ -754,7 +758,7 @@ class SICELIIntegrator(SIIntegrator): simulations """ proc_time, eos_conc = timed_deplete( - self.chain, bos_conc, bos_rates, dt) + CRAM48, self.chain, bos_conc, bos_rates, dt) inter_conc = copy.deepcopy(eos_conc) # Begin iteration @@ -770,9 +774,9 @@ class SICELIIntegrator(SIIntegrator): list_rates = list(zip(bos_rates, res_bar.rates)) time1, inter_conc = timed_deplete( - self.chain, bos_conc, list_rates, dt, matrix_func=celi_f1) + CRAM48, self.chain, bos_conc, list_rates, dt, matrix_func=celi_f1) time2, inter_conc = timed_deplete( - self.chain, inter_conc, list_rates, dt, matrix_func=celi_f2) + CRAM48, self.chain, inter_conc, list_rates, dt, matrix_func=celi_f2) proc_time += time1 + time2 # end iteration @@ -880,9 +884,9 @@ class SILEQIIntegrator(SIIntegrator): inputs = list(zip(self._prev_rates, bos_rates, repeat(prev_dt), repeat(dt))) proc_time, inter_conc = timed_deplete( - self.chain, bos_conc, inputs, dt, matrix_func=leqi_f1) + CRAM48, self.chain, bos_conc, inputs, dt, matrix_func=leqi_f1) time1, eos_conc = timed_deplete( - self.chain, inter_conc, inputs, dt, matrix_func=leqi_f2) + CRAM48, self.chain, inter_conc, inputs, dt, matrix_func=leqi_f2) proc_time += time1 inter_conc = copy.deepcopy(eos_conc) @@ -900,9 +904,9 @@ class SILEQIIntegrator(SIIntegrator): inputs = list(zip(self._prev_rates, bos_rates, res_bar.rates, repeat(prev_dt), repeat(dt))) time1, inter_conc = timed_deplete( - self.chain, bos_conc, inputs, dt, matrix_func=leqi_f3) + CRAM48, self.chain, bos_conc, inputs, dt, matrix_func=leqi_f3) time2, inter_conc = timed_deplete( - self.chain, inter_conc, inputs, dt, matrix_func=leqi_f4) + CRAM48, self.chain, inter_conc, inputs, dt, matrix_func=leqi_f4) proc_time += time1 + time2 return proc_time, [eos_conc, inter_conc], [res_bar] diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 36791b0182..9506ad3175 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -414,7 +414,7 @@ def test_fission_yield_attribute(simple_chain): dummy_conc = [[1, 2]] * (len(empty_chain.fission_yields) + 1) with pytest.raises( ValueError, match="fission yield.*not equal.*compositions"): - cram.deplete(empty_chain, dummy_conc, None, 0.5) + cram.deplete(cram.CRAM48, empty_chain, dummy_conc, None, 0.5) def test_validate(simple_chain): From 46e98e4a01d71ebf67cb65e4630df540b657e197 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 9 May 2020 14:38:46 -0400 Subject: [PATCH 2/6] Move deplete function to dedicated pool module; removed timed_deplete --- docs/source/pythonapi/deplete.rst | 3 +- openmc/deplete/cram.py | 77 +------------------------- openmc/deplete/pool.py | 57 +++++++++++++++++++ tests/unit_tests/test_deplete_chain.py | 4 +- 4 files changed, 61 insertions(+), 80 deletions(-) create mode 100644 openmc/deplete/pool.py diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index b9970ffb11..b01de477ce 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -143,8 +143,7 @@ with :func:`cram.CRAM48` being the default. cram.CRAM16 cram.CRAM48 - cram.deplete - cram.timed_deplete + pool.deplete The following classes are used to help the :class:`openmc.deplete.Operator` compute quantities like effective fission yields, reaction rates, and diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py index 939ff5432c..6b6a909ec1 100644 --- a/openmc/deplete/cram.py +++ b/openmc/deplete/cram.py @@ -4,9 +4,6 @@ Implements two different forms of CRAM for use in openmc.deplete. """ import numbers -from itertools import repeat -from multiprocessing import Pool -import time import numpy as np import scipy.sparse as sp @@ -15,79 +12,7 @@ import scipy.sparse.linalg as sla from openmc.checkvalue import check_type, check_length from .abc import DepSystemSolver -__all__ = [ - "deplete", "timed_deplete", "CRAM16", "CRAM48", - "Cram16Solver", "Cram48Solver", "IPFCramSolver"] - - -def deplete(func, chain, x, rates, dt, matrix_func=None): - """Deplete materials using given reaction rates for a specified time - - Parameters - ---------- - func : callable - Function to use to get new compositions. Expected to have the - signature ``func(A, n0, t) -> n1`` - chain : openmc.deplete.Chain - Depletion chain - x : list of numpy.ndarray - Atom number vectors for each material - rates : openmc.deplete.ReactionRates - Reaction rates (from transport operator) - dt : float - Time in [s] to deplete for - maxtrix_func : Callable, optional - Function to form the depletion matrix after calling - ``matrix_func(chain, rates, fission_yields)``, where - ``fission_yields = {parent: {product: yield_frac}}`` - Expected to return the depletion matrix required by - :func:`CRAM48`. - - Returns - ------- - x_result : list of numpy.ndarray - Updated atom number vectors for each material - """ - - fission_yields = chain.fission_yields - if len(fission_yields) == 1: - fission_yields = repeat(fission_yields[0]) - elif len(fission_yields) != len(x): - raise ValueError( - "Number of material fission yield distributions {} is not equal " - "to the number of compositions {}".format(len(fission_yields), - len(x))) - - if matrix_func is None: - matrices = map(chain.form_matrix, rates, fission_yields) - else: - matrices = map(matrix_func, repeat(chain), rates, fission_yields) - - # Use multiprocessing pool to distribute work - with Pool() as pool: - inputs = zip(matrices, x, repeat(dt)) - x_result = list(pool.starmap(func, inputs)) - - return x_result - - -def timed_deplete(*args, **kwargs): - """Wrapper over :func:`deplete` that also returns process time - - All arguments and keyword arguments are passed onto - :func:`deplete` directly. - - Returns - ------- - proc_time: float - Process time required to return from deplete - results: list of numpy arrays - Output from :func:`deplete` call - """ - - start = time.time() - results = deplete(*args, **kwargs) - return time.time() - start, results +__all__ = ["CRAM16", "CRAM48", "Cram16Solver", "Cram48Solver", "IPFCramSolver"] class IPFCramSolver(DepSystemSolver): diff --git a/openmc/deplete/pool.py b/openmc/deplete/pool.py new file mode 100644 index 0000000000..ef92569b73 --- /dev/null +++ b/openmc/deplete/pool.py @@ -0,0 +1,57 @@ +"""Dedicated module containing depletion function + +Provided to avoid some circular imports +""" +from itertools import repeat +from multiprocessing import Pool + + +def deplete(func, chain, x, rates, dt, matrix_func=None): + """Deplete materials using given reaction rates for a specified time + + Parameters + ---------- + func : callable + Function to use to get new compositions. Expected to have the + signature ``func(A, n0, t) -> n1`` + chain : openmc.deplete.Chain + Depletion chain + x : list of numpy.ndarray + Atom number vectors for each material + rates : openmc.deplete.ReactionRates + Reaction rates (from transport operator) + dt : float + Time in [s] to deplete for + maxtrix_func : Callable, optional + Function to form the depletion matrix after calling + ``matrix_func(chain, rates, fission_yields)``, where + ``fission_yields = {parent: {product: yield_frac}}`` + Expected to return the depletion matrix required by + :func:`CRAM48`. + + Returns + ------- + x_result : list of numpy.ndarray + Updated atom number vectors for each material + """ + + fission_yields = chain.fission_yields + if len(fission_yields) == 1: + fission_yields = repeat(fission_yields[0]) + elif len(fission_yields) != len(x): + raise ValueError( + "Number of material fission yield distributions {} is not equal " + "to the number of compositions {}".format( + len(fission_yields), len(x))) + + if matrix_func is None: + matrices = map(chain.form_matrix, rates, fission_yields) + else: + matrices = map(matrix_func, repeat(chain), rates, fission_yields) + + # Use multiprocessing pool to distribute work + with Pool() as pool: + inputs = zip(matrices, x, repeat(dt)) + x_result = list(pool.starmap(func, inputs)) + + return x_result diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 9506ad3175..7162be6de9 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -7,7 +7,7 @@ import os from pathlib import Path import numpy as np -from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram +from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram, pool import pytest from tests import cdtemp @@ -414,7 +414,7 @@ def test_fission_yield_attribute(simple_chain): dummy_conc = [[1, 2]] * (len(empty_chain.fission_yields) + 1) with pytest.raises( ValueError, match="fission yield.*not equal.*compositions"): - cram.deplete(cram.CRAM48, empty_chain, dummy_conc, None, 0.5) + pool.deplete(cram.CRAM48, empty_chain, dummy_conc, None, 0.5) def test_validate(simple_chain): From 8d2d11bf49ac0c67faf69ebdd5b35ef8210f2c8b Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 9 May 2020 14:54:03 -0400 Subject: [PATCH 3/6] Configure depletion solver with Integrator argument, attribute Users can pass a string or callable function to the Integrator construction to configure the depletion solver, defaulting to CRAM48. Supported strings are "cram16" and "cram48". Callables are passed to the new solver attribute, and must accept three input arguments for A, n0, and dt. Functions are expected to return an array of compositions for this material. Imports of CRAM48 and CRAM16 are delayed until inside the class since the cram module imports from abc and can lead to a circular import. --- openmc/deplete/abc.py | 145 ++++++++++++++-- openmc/deplete/integrators.py | 313 ++++++++++++++++++++++++++++------ 2 files changed, 389 insertions(+), 69 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index ddc48bc4de..cbccd05ee4 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 ad0137221b..020291555b 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] From a253bcccb17e037376094173e1e284dbc944b0f4 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 9 May 2020 15:20:35 -0400 Subject: [PATCH 4/6] Add tests for configurable integrator solver --- tests/unit_tests/test_deplete_integrator.py | 28 ++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index 6c09b3feac..e060b091d1 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -17,7 +17,8 @@ import pytest from openmc.deplete import ( ReactionRates, Results, ResultsList, comm, OperatorResult, PredictorIntegrator, CECMIntegrator, CF4Integrator, CELIIntegrator, - EPCRK4Integrator, LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator) + EPCRK4Integrator, LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator, + cram) from tests import dummy_operator @@ -147,6 +148,23 @@ def test_bad_integrator_inputs(): with pytest.raises(ValueError, match="n_steps"): SICELIIntegrator(op, timesteps, [1], n_steps=0) + with pytest.raises(ValueError, match="Solver failure"): + PredictorIntegrator(op, timesteps, power=1, solver="failure") + + with pytest.raises(TypeError, match=".*callable.*NoneType"): + PredictorIntegrator(op, timesteps, power=1, solver=None) + + with pytest.raises(ValueError, match=".*arguments"): + PredictorIntegrator(op, timesteps, power=1, solver=mock_bad_solver_nargs) + + +def mock_good_solver(A, n, t): + pass + + +def mock_bad_solver_nargs(A, n): + pass + @pytest.mark.parametrize("scheme", dummy_operator.SCHEMES) def test_integrator(run_in_tmpdir, scheme): @@ -174,6 +192,14 @@ def test_integrator(run_in_tmpdir, scheme): assert dep_time.shape == (2, ) assert all(dep_time > 0) + integrator = bundle.solver(operator, [0.75], 1, solver=cram.CRAM48) + assert integrator.solver is cram.CRAM48 + + integrator.solver = mock_good_solver + assert integrator.solver is mock_good_solver + + lfunc = lambda A, n, t: mock_good_solver(A, n, t) + @pytest.mark.parametrize("integrator", INTEGRATORS) def test_timesteps(integrator): From 7fdfe0ec0b59db10ab2eca59bc3de1931aaecfb5 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 14 May 2020 08:11:01 -0400 Subject: [PATCH 5/6] Documentation fixes for configurable CRAM solver Recommendations and typos found through review with @paulromano --- openmc/deplete/abc.py | 16 ++++---- openmc/deplete/integrators.py | 76 +++++++++++++++++++---------------- openmc/deplete/pool.py | 5 ++- 3 files changed, 52 insertions(+), 45 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index cbccd05ee4..8846906b44 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -626,7 +626,7 @@ class Integrator(ABC): 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 + If a string, must be the name of the solver responsible for solving the Bateman equations. Current options are: * ``cram16`` - 16th order IPF CRAM @@ -649,9 +649,9 @@ class Integrator(ABC): 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: + :math:`\frac{\partial}{\partial t}\vec{n} = 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 @@ -925,7 +925,7 @@ class SIIntegrator(Integrator): 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 + If a string, must be the name of the solver responsible for solving the Bateman equations. Current options are: * ``cram16`` - 16th order IPF CRAM @@ -950,9 +950,9 @@ class SIIntegrator(Integrator): 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: + :math:`\frac{\partial}{\partial t}\vec{n} = 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 diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index 020291555b..1fe6369003 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -55,7 +55,7 @@ class PredictorIntegrator(Integrator): .. versionadded:: 0.12 solver : str or callable, optional - If a string, must be the tame of the solver responsible for + If a string, must be the name of the solver responsible for solving the Bateman equations. Current options are: * ``cram16`` - 16th order IPF CRAM @@ -78,9 +78,9 @@ class PredictorIntegrator(Integrator): 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: + :math:`\frac{\partial}{\partial t}\vec{n} = 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 @@ -174,7 +174,7 @@ class CECMIntegrator(Integrator): .. versionadded:: 0.12 solver : str or callable, optional - If a string, must be the tame of the solver responsible for + If a string, must be the name of the solver responsible for solving the Bateman equations. Current options are: * ``cram16`` - 16th order IPF CRAM @@ -197,9 +197,9 @@ class CECMIntegrator(Integrator): 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: + :math:`\frac{\partial}{\partial t}\vec{n} = 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 @@ -301,7 +301,7 @@ class CF4Integrator(Integrator): .. versionadded:: 0.12 solver : str or callable, optional - If a string, must be the tame of the solver responsible for + If a string, must be the name of the solver responsible for solving the Bateman equations. Current options are: * ``cram16`` - 16th order IPF CRAM @@ -324,9 +324,9 @@ class CF4Integrator(Integrator): 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: + :math:`\frac{\partial}{\partial t}\vec{n} = 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 @@ -445,7 +445,7 @@ class CELIIntegrator(Integrator): .. versionadded:: 0.12 solver : str or callable, optional - If a string, must be the tame of the solver responsible for + If a string, must be the name of the solver responsible for solving the Bateman equations. Current options are: * ``cram16`` - 16th order IPF CRAM @@ -468,9 +468,9 @@ class CELIIntegrator(Integrator): 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: + :math:`\frac{\partial}{\partial t}\vec{n} = 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 @@ -576,7 +576,7 @@ class EPCRK4Integrator(Integrator): .. versionadded:: 0.12 solver : str or callable, optional - If a string, must be the tame of the solver responsible for + If a string, must be the name of the solver responsible for solving the Bateman equations. Current options are: * ``cram16`` - 16th order IPF CRAM @@ -598,7 +598,7 @@ class EPCRK4Integrator(Integrator): 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 + If a string, must be the name of the solver responsible for solving the Bateman equations. Current options are: * ``cram16`` - 16th order IPF CRAM @@ -721,7 +721,7 @@ class LEQIIntegrator(Integrator): .. versionadded:: 0.12 solver : str or callable, optional - If a string, must be the tame of the solver responsible for + If a string, must be the name of the solver responsible for solving the Bateman equations. Current options are: * ``cram16`` - 16th order IPF CRAM @@ -744,9 +744,9 @@ class LEQIIntegrator(Integrator): 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: + :math:`\frac{\partial}{\partial t}\vec{n} = 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 @@ -872,7 +872,7 @@ class SICELIIntegrator(SIIntegrator): 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 + If a string, must be the name of the solver responsible for solving the Bateman equations. Current options are: * ``cram16`` - 16th order IPF CRAM @@ -897,9 +897,9 @@ class SICELIIntegrator(SIIntegrator): 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: + :math:`\frac{\partial}{\partial t}\vec{n} = 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 @@ -1009,7 +1009,7 @@ class SILEQIIntegrator(SIIntegrator): 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 + If a string, must be the name of the solver responsible for solving the Bateman equations. Current options are: * ``cram16`` - 16th order IPF CRAM @@ -1032,18 +1032,24 @@ 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: + solver : callable + Function that will solve the Bateman equations + :math:`\frac{\partial}{\partial t}\vec{n} = 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 - * ``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`. + * ``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 diff --git a/openmc/deplete/pool.py b/openmc/deplete/pool.py index ef92569b73..d9d3597751 100644 --- a/openmc/deplete/pool.py +++ b/openmc/deplete/pool.py @@ -22,17 +22,18 @@ def deplete(func, chain, x, rates, dt, matrix_func=None): Reaction rates (from transport operator) dt : float Time in [s] to deplete for - maxtrix_func : Callable, optional + maxtrix_func : callable, optional Function to form the depletion matrix after calling ``matrix_func(chain, rates, fission_yields)``, where ``fission_yields = {parent: {product: yield_frac}}`` Expected to return the depletion matrix required by - :func:`CRAM48`. + ``func`` Returns ------- x_result : list of numpy.ndarray Updated atom number vectors for each material + """ fission_yields = chain.fission_yields From 133bffcb28f32e845033688d1a87e0d31a449951 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 14 May 2020 08:16:08 -0400 Subject: [PATCH 6/6] Expand integrator.solver test: lambda, named solver --- tests/unit_tests/test_deplete_integrator.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index e060b091d1..0560cb9f73 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -195,10 +195,15 @@ def test_integrator(run_in_tmpdir, scheme): integrator = bundle.solver(operator, [0.75], 1, solver=cram.CRAM48) assert integrator.solver is cram.CRAM48 + integrator = bundle.solver(operator, [0.75], 1, solver="cram16") + assert integrator.solver is cram.CRAM16 + integrator.solver = mock_good_solver assert integrator.solver is mock_good_solver lfunc = lambda A, n, t: mock_good_solver(A, n, t) + integrator.solver = lfunc + assert integrator.solver is lfunc @pytest.mark.parametrize("integrator", INTEGRATORS)