diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index a72504ac92..5b26eb38aa 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -6,7 +6,7 @@ .. module:: openmc.deplete -Several functions are provided that implement different time-integration +Several classes are provided that implement different time-integration algorithms for depletion calculations, which are described in detail in Colin Josey's thesis, `Development and analysis of high order neutron transport-depletion coupling algorithms `_. @@ -25,7 +25,7 @@ transport-depletion coupling algorithms `_. SICELIIntegrator SILEQIIntegrator -Each of these functions expects a "transport operator" to be passed. An operator +Each of these classes expects a "transport operator" to be passed. An operator specific to OpenMC is available using the following class: .. autosummary:: @@ -65,6 +65,8 @@ for a depletion chain: DecayTuple Nuclide ReactionTuple + FissionYieldDistribution + FissionYield The following classes are used during a depletion simulation and store auxiliary data, such as number densities and reaction rates for each material. @@ -75,13 +77,25 @@ data, such as number densities and reaction rates for each material. :template: myclass.rst AtomNumber - ChainFissionHelper - DirectReactionRateHelper OperatorResult ReactionRates Results ResultsList +The following classes are used to help the :class:`openmc.deplete.Operator` +compute quantities like effective fission yields, reaction rates, and +total system energy. + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + helpers.AveragedFissionYieldHelper + helpers.ChainFissionHelper + helpers.ConstantFissionYieldHelper + helpers.DirectReactionRateHelper + helpers.FissionYieldCutoffHelper The following classes are abstract classes that can be used to extend the :mod:`openmc.deplete` capabilities: @@ -91,8 +105,10 @@ The following classes are abstract classes that can be used to extend the :nosignatures: :template: myclass.rst - ReactionRateHelper EnergyHelper + FissionYieldHelper + ReactionRateHelper + TalliedFissionYieldHelper TransportOperator Custom integrators can be developed by subclassing from the following abstract diff --git a/openmc/deplete/_matrix_funcs.py b/openmc/deplete/_matrix_funcs.py index 1cacbd2c06..a606d739f7 100644 --- a/openmc/deplete/_matrix_funcs.py +++ b/openmc/deplete/_matrix_funcs.py @@ -1,77 +1,78 @@ """Functions to form the special matrix for depletion""" -def celi_f1(chain, rates): - return (5 / 12 * chain.form_matrix(rates[0]) - + 1 / 12 * chain.form_matrix(rates[1])) +def celi_f1(chain, rates, fission_yields=None): + return (5 / 12 * chain.form_matrix(rates[0], fission_yields) + + 1 / 12 * chain.form_matrix(rates[1], fission_yields)) -def celi_f2(chain, rates): - return (1 / 12 * chain.form_matrix(rates[0]) - + 5 / 12 * chain.form_matrix(rates[1])) +def celi_f2(chain, rates, fission_yields=None): + return (1 / 12 * chain.form_matrix(rates[0], fission_yields) + + 5 / 12 * chain.form_matrix(rates[1], fission_yields)) -def cf4_f1(chain, rates): - return 1 / 2 * chain.form_matrix(rates) +def cf4_f1(chain, rates, fission_yields=None): + return 1 / 2 * chain.form_matrix(rates, fission_yields) -def cf4_f2(chain, rates): - return -1 / 2 * chain.form_matrix(rates[0]) + chain.form_matrix(rates[1]) +def cf4_f2(chain, rates, fission_yields=None): + return (-1 / 2 * chain.form_matrix(rates[0], fission_yields) + + chain.form_matrix(rates[1], fission_yields)) -def cf4_f3(chain, rates): - return (1 / 4 * chain.form_matrix(rates[0]) - + 1 / 6 * chain.form_matrix(rates[1]) - + 1 / 6 * chain.form_matrix(rates[2]) - - 1 / 12 * chain.form_matrix(rates[3])) +def cf4_f3(chain, rates, fission_yields=None): + return (1 / 4 * chain.form_matrix(rates[0], fission_yields) + + 1 / 6 * chain.form_matrix(rates[1], fission_yields) + + 1 / 6 * chain.form_matrix(rates[2], fission_yields) + - 1 / 12 * chain.form_matrix(rates[3], fission_yields)) -def cf4_f4(chain, rates): - return (-1 / 12 * chain.form_matrix(rates[0]) - + 1 / 6 * chain.form_matrix(rates[1]) - + 1 / 6 * chain.form_matrix(rates[2]) - + 1 / 4 * chain.form_matrix(rates[3])) +def cf4_f4(chain, rates, fission_yields=None): + return (-1 / 12 * chain.form_matrix(rates[0], fission_yields) + + 1 / 6 * chain.form_matrix(rates[1], fission_yields) + + 1 / 6 * chain.form_matrix(rates[2], fission_yields) + + 1 / 4 * chain.form_matrix(rates[3], fission_yields)) -def rk4_f1(chain, rates): - return 1 / 2 * chain.form_matrix(rates) +def rk4_f1(chain, rates, fission_yields=None): + return 1 / 2 * chain.form_matrix(rates, fission_yields) -def rk4_f4(chain, rates): - return (1 / 6 * chain.form_matrix(rates[0]) - + 1 / 3 * chain.form_matrix(rates[1]) - + 1 / 3 * chain.form_matrix(rates[2]) - + 1 / 6 * chain.form_matrix(rates[3])) +def rk4_f4(chain, rates, fission_yields=None): + return (1 / 6 * chain.form_matrix(rates[0], fission_yields) + + 1 / 3 * chain.form_matrix(rates[1], fission_yields) + + 1 / 3 * chain.form_matrix(rates[2], fission_yields) + + 1 / 6 * chain.form_matrix(rates[3], fission_yields)) -def leqi_f1(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) +def leqi_f1(chain, inputs, fission_yields): + f1 = chain.form_matrix(inputs[0], fission_yields) + f2 = chain.form_matrix(inputs[1], fission_yields) dt_l, dt = inputs[2], inputs[3] return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 -def leqi_f2(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) +def leqi_f2(chain, inputs, fission_yields=None): + f1 = chain.form_matrix(inputs[0], fission_yields) + f2 = chain.form_matrix(inputs[1], fission_yields) dt_l, dt = inputs[2], inputs[3] return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 -def leqi_f3(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - f3 = chain.form_matrix(inputs[2]) +def leqi_f3(chain, inputs, fission_yields=None): + f1 = chain.form_matrix(inputs[0], fission_yields) + f2 = chain.form_matrix(inputs[1], fission_yields) + f3 = chain.form_matrix(inputs[2], fission_yields) dt_l, dt = inputs[3], inputs[4] return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1 + (dt ** 2 + 6 * dt * dt_l + 5 * dt_l ** 2) / (12 * dt_l * (dt + dt_l)) * f2 + dt_l / (12 * (dt + dt_l)) * f3) -def leqi_f4(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - f3 = chain.form_matrix(inputs[2]) +def leqi_f4(chain, inputs, fission_yields=None): + f1 = chain.form_matrix(inputs[0], fission_yields) + f2 = chain.form_matrix(inputs[1], fission_yields) + f3 = chain.form_matrix(inputs[2], fission_yields) dt_l, dt = inputs[3], inputs[4] return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1 + (dt ** 2 + 2 * dt * dt_l + dt_l ** 2) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 0f1203aeb8..fb638e5c21 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -13,10 +13,11 @@ from copy import deepcopy from warnings import warn from numbers import Real, Integral -from numpy import nonzero, empty +from numpy import nonzero, empty, asarray from uncertainties import ufloat from openmc.data import DataLibrary, JOULE_PER_EV +from openmc.capi import MaterialFilter, Tally from openmc.checkvalue import check_type, check_greater_than from .results import Results from .chain import Chain @@ -364,6 +365,217 @@ class EnergyHelper(ABC): self._nuclides = nuclides +class FissionYieldHelper(ABC): + """Abstract class for processing energy dependent fission yields + + Parameters + ---------- + chain_nuclides : iterable of openmc.deplete.Nuclide + Nuclides tracked in the depletion chain. All nuclides are + not required to have fission yield data. + + Attributes + ---------- + constant_yields : dict of str to :class:`openmc.deplete.FissionYield` + Fission yields for all nuclides that only have one set of + fission yield data. Can be accessed as ``{parent: {product: yield}}`` + """ + + def __init__(self, chain_nuclides): + self._chain_nuclides = {} + self._constant_yields = {} + + # Get all nuclides with fission yield data + for nuc in chain_nuclides: + if nuc.yield_data is None: + continue + if len(nuc.yield_data) == 1: + self._constant_yields[nuc.name] = ( + nuc.yield_data[nuc.yield_energies[0]]) + elif len(nuc.yield_data) > 1: + self._chain_nuclides[nuc.name] = nuc + self._chain_set = set(self._chain_nuclides) | set(self._constant_yields) + + @property + def constant_yields(self): + return deepcopy(self._constant_yields) + + @abstractmethod + def weighted_yields(self, local_mat_index): + """Return fission yields for a specific material + + Parameters + ---------- + local_mat_index : int + Index for material tracked on this process that + exists in :attr:`local_mat_index` and fits within + the first axis in :attr:`results` + + Returns + ------- + library : dict + Dictionary of ``{parent: {product: fyield}}`` + """ + + @staticmethod + def unpack(): + """Unpack tally data prior to compute fission yields. + + Called after a :meth:`openmc.deplete.Operator.__call__` + routine during the normalization of reaction rates. + + Not necessary for all subclasses to implement, unless tallies + are used. + """ + + @staticmethod + def generate_tallies(materials, mat_indexes): + """Construct tallies necessary for computing fission yields + + Called during the operator set up phase prior to depleting. + Not necessary for subclasses to implement + + Parameters + ---------- + materials : iterable of C-API materials + Materials to be used in :class:`openmc.capi.MaterialFilter` + mat_indexes : iterable of int + Indices of tallied materials that will have their fission + yields computed by this helper. Necessary as the + :class:`openmc.deplete.Operator` that uses this helper + may only burn a subset of all materials when running + in parallel mode. + """ + + def update_tally_nuclides(self, nuclides): + """Return nuclides with non-zero densities and yield data + + Parameters + ---------- + nuclides : iterable of str + Nuclides with non-zero densities from the + :class:`openmc.deplete.Operator` + + Returns + ------- + nuclides : list of str + Union of nuclides that the :class:`openmc.deplete.Operator` + says have non-zero densities at this stage and those that + have yield data. Sorted by nuclide name + + """ + return sorted(self._chain_set & set(nuclides)) + + @classmethod + def from_operator(cls, operator, **kwargs): + """Create a new instance by pulling data from the operator + + All keyword arguments should be identical to their counterpart + in the main ``__init__`` method + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + Operator with a depletion chain + kwargs: optional + Additional keyword arguments to be used in constuction + """ + return cls(operator.chain.nuclides, **kwargs) + + +class TalliedFissionYieldHelper(FissionYieldHelper): + """Abstract class for computing fission yields with tallies + + Generates a basic fission rate tally in all burnable materials with + :meth:`generate_tallies`, and set nuclides to be tallied with + :meth:`update_tally_nuclides`. Subclasses will need to implement + :meth:`unpack` and :meth:`weighted_yields`. + + Parameters + ---------- + chain_nuclides : iterable of openmc.deplete.Nuclide + Nuclides tracked in the depletion chain. Not necessary + that all have yield data. + + Attributes + ---------- + constant_yields : dict of str to :class:`openmc.deplete.FissionYield` + Fission yields for all nuclides that only have one set of + fission yield data. Can be accessed as ``{parent: {product: yield}}`` + results : None or numpy.ndarray + Tally results shaped in a manner useful to this helper. + """ + + _upper_energy = 20.0e6 # upper energy for tallies + + def __init__(self, chain_nuclides): + super().__init__(chain_nuclides) + self._local_indexes = None + self._fission_rate_tally = None + self._tally_nucs = [] + self.results = None + + def generate_tallies(self, materials, mat_indexes): + """Construct the fission rate tally + + Parameters + ---------- + materials : iterable of :class:`openmc.capi.Material` + Materials to be used in :class:`openmc.capi.MaterialFilter` + mat_indexes : iterable of int + Indices of tallied materials that will have their fission + yields computed by this helper. Necessary as the + :class:`openmc.deplete.Operator` that uses this helper + may only burn a subset of all materials when running + in parallel mode. + """ + self._local_indexes = asarray(mat_indexes) + + # Tally group-wise fission reaction rates + self._fission_rate_tally = Tally() + self._fission_rate_tally.scores = ['fission'] + + self._fission_rate_tally.filters = [MaterialFilter(materials)] + + def update_tally_nuclides(self, nuclides): + """Tally nuclides with non-zero density and multiple yields + + Must be run after :meth:`generate_tallies`. + + Parameters + ---------- + nuclides : iterable of str + Potential nuclides to be tallied, such as those with + non-zero density at this stage. + + Returns + ------- + nuclides : list of str + Union of input nuclides and those that have multiple sets + of yield data. Sorted by nuclide name + + Raises + ------ + AttributeError + If tallies not generated + """ + assert self._fission_rate_tally is not None, ( + "Run generate_tallies first") + overlap = set(self._chain_nuclides).intersection(set(nuclides)) + nuclides = sorted(overlap) + self._tally_nucs = [self._chain_nuclides[n] for n in nuclides] + self._fission_rate_tally.nuclides = nuclides + return nuclides + + @abstractmethod + def unpack(self): + """Unpack tallies after a transport run. + + Abstract because each subclass will need to arrange its + tally data. + """ + + class Integrator(ABC): """Abstract class for solving the time-integration for depletion diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index b3e2dae29e..c2b5b75b31 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -9,12 +9,13 @@ from itertools import chain import math import re from collections import OrderedDict, defaultdict -from collections.abc import Mapping +from collections.abc import Mapping, Iterable from numbers import Real from warnings import warn from openmc.checkvalue import check_type, check_greater_than from openmc.data import gnd_name, zam +from .nuclide import FissionYieldDistribution # Try to use lxml if it is available. It preserves the order of attributes and # provides a pretty-printer by default. If not available, @@ -136,12 +137,22 @@ class Chain(object): Reactions that are tracked in the depletion chain nuclide_dict : OrderedDict of str to int Maps a nuclide name to an index in nuclides. + fission_yields : None or iterable of dict + List of effective fission yields for materials. Each dictionary + should be of the form ``{parent: {product: yield}}`` with + types ``{str: {str: float}}``, where ``yield`` is the fission product + yield for isotope ``parent`` producing isotope ``product``. + A single entry indicates yields are constant across all materials. + Otherwise, an entry can be added for each material to be burned. + Ordering should be identical to how the operator orders reaction + rates for burnable materials. """ def __init__(self): self.nuclides = [] self.reactions = [] self.nuclide_dict = OrderedDict() + self._fission_yields = None def __contains__(self, nuclide): return nuclide in self.nuclide_dict @@ -209,8 +220,7 @@ class Chain(object): for idx, parent in enumerate(sorted(decay_data, key=openmc.data.zam)): data = decay_data[parent] - nuclide = Nuclide() - nuclide.name = parent + nuclide = Nuclide(parent) chain.nuclides.append(nuclide) chain.nuclide_dict[parent] = idx @@ -225,7 +235,8 @@ class Chain(object): if mode.daughter in decay_data: target = mode.daughter else: - print('missing {} {} {}'.format(parent, ','.join(mode.modes), mode.daughter)) + print('missing {} {} {}'.format( + parent, ','.join(mode.modes), mode.daughter)) target = replace_missing(mode.daughter, decay_data) # Write branching ratio, taking care to ensure sum is unity @@ -279,15 +290,16 @@ class Chain(object): fpy = fpy_data[parent] if fpy.energies is not None: - nuclide.yield_energies = fpy.energies + yield_energies = fpy.energies else: - nuclide.yield_energies = [0.0] + yield_energies = [0.0] - for E, table in zip(nuclide.yield_energies, fpy.independent): + yield_data = {} + for E, table in zip(yield_energies, fpy.independent): yield_replace = 0.0 yields = defaultdict(float) for product, y in table.items(): - # Handle fission products that have no decay data available + # Handle fission products that have no decay data if product not in decay_data: daughter = replace_missing(product, decay_data) product = daughter @@ -297,10 +309,9 @@ class Chain(object): if yield_replace > 0.0: missing_fp.append((parent, E, yield_replace)) + yield_data[E] = yields - nuclide.yield_data[E] = [] - for k in sorted(yields, key=openmc.data.zam): - nuclide.yield_data[E].append((k, yields[k])) + nuclide.yield_data = FissionYieldDistribution(yield_data) # Display warnings if missing_daughter: @@ -387,23 +398,56 @@ class Chain(object): clean_indentation(root_elem) tree.write(str(filename), encoding='utf-8') - def form_matrix(self, rates): + def get_thermal_fission_yields(self): + """Return fission yields at lowest incident neutron energy + + Used as the default set of fission yields for :meth:`form_matrix` + if ``fission_yields`` are not provided + + Returns + ------- + fission_yields : dict + Dictionary of ``{parent: {product: f_yield}}`` + where ``parent`` and ``product`` are both string + names of nuclides with yield data and ``f_yield`` + is a float for the fission yield. + """ + out = {} + for nuc in self.nuclides: + if nuc.yield_data is None: + continue + yield_obj = nuc.yield_data[min(nuc.yield_energies)] + out[nuc.name] = dict(yield_obj) + return out + + def form_matrix(self, rates, fission_yields=None): """Forms depletion matrix. Parameters ---------- rates : numpy.ndarray 2D array indexed by (nuclide, reaction) + fission_yields : dict, optional + Option to use a custom set of fission yields. Expected + to be of the form ``{parent : {product : f_yield}}`` + with string nuclide names for ``parent`` and ``product``, + and ``f_yield`` as the respective fission yield Returns ------- scipy.sparse.csr_matrix Sparse matrix representing depletion. + See Also + -------- + :meth:`get_thermal_fission_yields` """ matrix = defaultdict(float) reactions = set() + if fission_yields is None: + fission_yields = self.get_thermal_fission_yields() + for i, nuc in enumerate(self.nuclides): if nuc.n_decay_modes != 0: @@ -448,11 +492,7 @@ class Chain(object): k = self.nuclide_dict[target] matrix[k, i] += path_rate * br else: - # Assume that we should always use thermal fission - # yields. At some point it would be nice to account - # for the energy-dependence.. - energy, data = sorted(nuc.yield_data.items())[0] - for product, y in data: + for product, y in fission_yields[nuc.name].items(): yield_val = y * path_rate if yield_val != 0.0: k = self.nuclide_dict[product] @@ -678,6 +718,20 @@ class Chain(object): parent.reactions.append(ReactionTuple( reaction, ground_tgt, rxn_Q, ground_br)) + @property + def fission_yields(self): + if self._fission_yields is None: + self._fission_yields = [self.get_thermal_fission_yields()] + return self._fission_yields + + @fission_yields.setter + def fission_yields(self, yields): + if yields is not None: + if isinstance(yields, Mapping): + yields = [yields] + check_type("fission_yields", yields, Iterable, Mapping) + self._fission_yields = yields + def validate(self, strict=True, quiet=False, tolerance=1e-4): """Search for possible inconsistencies diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py index 41a4e4c710..184c411bf3 100644 --- a/openmc/deplete/cram.py +++ b/openmc/deplete/cram.py @@ -30,7 +30,9 @@ def deplete(chain, x, rates, dt, matrix_func=None): dt : float Time in [s] to deplete for maxtrix_func : Callable, optional - Function of two variables: ``chain`` and ``rates``. + 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`. @@ -40,9 +42,19 @@ def deplete(chain, x, rates, dt, matrix_func=None): 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))) + # Use multiprocessing pool to distribute work with Pool() as pool: - iters = zip(repeat(chain), x, rates, repeat(dt), repeat(matrix_func)) + iters = zip(repeat(chain), x, rates, repeat(dt), + fission_yields, repeat(matrix_func)) x_result = list(pool.starmap(_cram_wrapper, iters)) return x_result @@ -67,7 +79,7 @@ def timed_deplete(*args, **kwargs): return time.time() - start, results -def _cram_wrapper(chain, n0, rates, dt, matrix_func=None): +def _cram_wrapper(chain, n0, rates, dt, fission_yields, matrix_func=None): """Wraps depletion matrix creation / CRAM solve for multiprocess execution Parameters @@ -82,6 +94,9 @@ def _cram_wrapper(chain, n0, rates, dt, matrix_func=None): Time to integrate to. maxtrix_func : function, optional Function to form the depletion matrix + fission_yields : dict + Single-energy fission yields of the form + ``{parent: {product: fission_yield}}`` Returns ------- @@ -90,9 +105,9 @@ def _cram_wrapper(chain, n0, rates, dt, matrix_func=None): """ if matrix_func is None: - A = chain.form_matrix(rates) + A = chain.form_matrix(rates, fission_yields) else: - A = matrix_func(chain, rates) + A = matrix_func(chain, rates, fission_yields) return CRAM48(A, n0, dt) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 66722c00a2..56644a3bba 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -1,12 +1,24 @@ """ Class for normalizing fission energy deposition """ +from copy import deepcopy from itertools import product +from numbers import Real +import bisect -from numpy import dot, zeros +from numpy import dot, zeros, newaxis -from openmc.capi import Tally, MaterialFilter -from .abc import ReactionRateHelper, EnergyHelper +from openmc.checkvalue import check_type, check_greater_than +from openmc.capi import ( + Tally, MaterialFilter, EnergyFilter, EnergyFunctionFilter) +from .abc import ( + ReactionRateHelper, EnergyHelper, FissionYieldHelper, + TalliedFissionYieldHelper) + +__all__ = ( + "DirectReactionRateHelper", "ChainFissionHelper", + "ConstantFissionYieldHelper", "FissionYieldCutoffHelper", + "AveragedFissionYieldHelper") # ------------------------------------- # Helpers for generating reaction rates @@ -142,3 +154,465 @@ class ChainFissionHelper(EnergyHelper): isotopes in all materials have the same Q value. """ self._energy += dot(fission_rates, self._fission_q_vector) + + +# ------------------------------------ +# Helper for collapsing fission yields +# ------------------------------------ + + +class ConstantFissionYieldHelper(FissionYieldHelper): + """Class that uses a single set of fission yields on each isotope + + Parameters + ---------- + chain_nuclides : iterable of openmc.deplete.Nuclide + Nuclides tracked in the depletion chain. All nuclides are + not required to have fission yield data. + energy : float, optional + Key in :attr:`openmc.deplete.Nuclide.yield_data` corresponding + to the desired set of fission yield data. Typically one of + ``{0.0253, 500000, 14000000}`` corresponding to 0.0253 eV, + 500 keV, and 14 MeV yield libraries. If the specific key is not + found, will fall back to closest energy present. + Default: 0.0253 eV for thermal yields + + Attributes + ---------- + constant_yields : dict of str to :class:`openmc.deplete.FissionYield` + Fission yields for all nuclides that only have one set of + fission yield data. Can be accessed as ``{parent: {product: yield}}`` + energy : float + Energy of fission yield libraries. + """ + + def __init__(self, chain_nuclides, energy=0.0253): + check_type("energy", energy, Real) + check_greater_than("energy", energy, 0.0, equality=True) + self._energy = energy + super().__init__(chain_nuclides) + # Iterate over all nuclides with > 1 set of yields + for name, nuc in self._chain_nuclides.items(): + yield_data = nuc.yield_data.get(energy) + if yield_data is not None: + self._constant_yields[name] = yield_data + continue + # Specific energy not found, use closest energy + distances = [abs(energy - ene) for ene in nuc.yield_energies] + min_E = min(nuc.yield_energies, key=lambda e: abs(e - energy)) + self._constant_yields[name] = nuc.yield_data[min_E] + + @classmethod + def from_operator(cls, operator, **kwargs): + """Return a new ConstantFissionYieldHelper using operator data + + All keyword arguments should be identical to their counterpart + in the main ``__init__`` method + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + operator with a depletion chain + kwargs: + Additional keyword arguments to be used in construction + + Returns + ------- + ConstantFissionYieldHelper + """ + return cls(operator.chain.nuclides, **kwargs) + + @property + def energy(self): + return self._energy + + def weighted_yields(self, _local_mat_index=None): + """Return fission yields for all nuclides requested + + Parameters + ---------- + _local_mat_index : int, optional + Current material index. Not used since all yields are + constant + + Returns + ------- + library : dict + Dictionary of ``{parent: {product: fyield}}`` + """ + return self.constant_yields + + +class FissionYieldCutoffHelper(TalliedFissionYieldHelper): + """Helper that computes fission yields based on a cutoff energy + + Tally fission rates above and below the cutoff energy. + Assume that all fissions below cutoff energy have use thermal fission + product yield distributions, while all fissions above use a faster + set of yield distributions. + + Uses a limit of 20 MeV for tallying fission. + + Parameters + ---------- + chain_nuclides : iterable of openmc.deplete.Nuclide + Nuclides tracked in the depletion chain. All nuclides are + not required to have fission yield data. + n_bmats : int, optional + Number of burnable materials tracked in the problem + cutoff : float, optional + Cutoff energy in [eV] below which all fissions will be + use thermal yields. All other fissions will use a + faster set of yields. Default: 112 [eV] + thermal_energy : float, optional + Energy of yield data corresponding to thermal yields. + Default: 0.0253 [eV] + fast_energy : float, optional + Energy of yield data corresponding to fast yields. + Default: 500 [kev] + + Attributes + ---------- + n_bmats : int + Number of burnable materials tracked in the problem. + Must be set prior to generating tallies + thermal_yields : dict + Dictionary of the form ``{parent: {product: yield}}`` + with thermal yields + fast_yields : dict + Dictionary of the form ``{parent: {product: yield}}`` + with fast yields + results : numpy.ndarray + Array of fission rate fractions with shape + ``(n_mats, 2, n_nucs)``. ``results[:, 0]`` + corresponds to the fraction of all fissions + that occured below ``cutoff``. The number + of materials in the first axis corresponds + to the number of materials burned by the + :class:`openmc.deplete.Operator` + """ + + def __init__(self, chain_nuclides, n_bmats, cutoff=112.0, + thermal_energy=0.0253, fast_energy=500.0e3): + check_type("cutoff", cutoff, Real) + check_type("thermal_energy", thermal_energy, Real) + check_type("fast_energy", fast_energy, Real) + check_greater_than("thermal_energy", thermal_energy, 0.0, equality=True) + check_greater_than("cutoff", cutoff, thermal_energy, equality=False) + check_greater_than("fast_energy", fast_energy, cutoff, equality=False) + self.n_bmats = n_bmats + super().__init__(chain_nuclides) + self._cutoff = cutoff + self._thermal_yields = {} + self._fast_yields = {} + convert_to_constant = set() + for name, nuc in self._chain_nuclides.items(): + yields = nuc.yield_data + energies = nuc.yield_energies + thermal = yields.get(thermal_energy) + fast = yields.get(fast_energy) + if thermal is None or fast is None: + if cutoff <= energies[0]: + # use lowest energy yields as constant + self._constant_yields[name] = yields[energies[0]] + convert_to_constant.add(name) + continue + if cutoff >= energies[-1]: + # use highest energy yields as constant + self._constant_yields[name] = yields[energies[-1]] + convert_to_constant.add(name) + continue + cutoff_ix = bisect.bisect_left(energies, cutoff) + # find closest energy to requested thermal, fast energies + if thermal is None: + min_E = min(energies[:cutoff_ix], + key=lambda e: abs(e - thermal_energy)) + thermal = yields[min_E] + if fast is None: + min_E = min(energies[cutoff_ix:], + key=lambda e: abs(e - fast_energy)) + fast = yields[min_E] + self._thermal_yields[name] = thermal + self._fast_yields[name] = fast + for name in convert_to_constant: + self._chain_nuclides.pop(name) + + @classmethod + def from_operator(cls, operator, **kwargs): + """Construct a helper from an operator + + All keyword arguments should be identical to their counterpart + in the main ``__init__`` method + + Parameters + ---------- + operator : openmc.deplete.Operator + Operator with a chain and burnable materials + kwargs: + Additional keyword arguments to be used in construction + + Returns + ------- + FissionYieldCutoffHelper + + """ + return cls(operator.chain.nuclides, len(operator.burnable_mats), + **kwargs) + + def generate_tallies(self, materials, mat_indexes): + """Use C API to produce a fission rate tally in burnable materials + + Include a :class:`openmc.capi.EnergyFilter` to tally fission rates + above and below cutoff energy. + + Parameters + ---------- + materials : iterable of :class:`openmc.capi.Material` + Materials to be used in :class:`openmc.capi.MaterialFilter` + mat_indexes : iterable of int + Indices of tallied materials that will have their fission + yields computed by this helper. Necessary as the + :class:`openmc.deplete.Operator` that uses this helper + may only burn a subset of all materials when running + in parallel mode. + """ + super().generate_tallies(materials, mat_indexes) + energy_filter = EnergyFilter([0.0, self._cutoff, self._upper_energy]) + self._fission_rate_tally.filters = ( + self._fission_rate_tally.filters + [energy_filter]) + + def unpack(self): + """Obtain fast and thermal fission fractions from tally""" + if not self._tally_nucs: + self.results = None + return + fission_rates = self._fission_rate_tally.results[..., 1].reshape( + self.n_bmats, 2, len(self._tally_nucs)) + self.results = fission_rates[self._local_indexes] + total_fission = self.results.sum(axis=1) + nz_mat, nz_nuc = total_fission.nonzero() + self.results[nz_mat, :, nz_nuc] /= total_fission[nz_mat, newaxis, nz_nuc] + + def weighted_yields(self, local_mat_index): + """Return fission yields for a specific material + + For nuclides with both yield data above and below + the cutoff energy, the effective yield for nuclide ``A`` + will be a weighted sum of fast and thermal yields. The + weights will be the fraction of ``A`` fission events + in the above and below the cutoff energy. + + If ``A`` has fission product distribution ``F`` + for fast fissions and ``T`` for thermal fissions, and + 70% of ``A`` fissions are considered thermal, then + the effective fission product yield distributions + for ``A`` is ``0.7 * T + 0.3 * F`` + + Parameters + ---------- + local_mat_index : int + Index for specific burnable material. Effective + yields will be produced using + ``self.results[local_mat_index]`` + + Returns + ------- + library : dict + Dictionary of ``{parent: {product: fyield}}`` + """ + yields = self.constant_yields + if not self._tally_nucs: + return yields + rates = self.results[local_mat_index] + # iterate over thermal then fast yields, prefer __mul__ to __rmul__ + for therm_frac, fast_frac, nuc in zip(rates[0], rates[1], self._tally_nucs): + yields[nuc.name] = (self._thermal_yields[nuc.name] * therm_frac + + self._fast_yields[nuc.name] * fast_frac) + return yields + + @property + def thermal_yields(self): + return deepcopy(self._thermal_yields) + + @property + def fast_yields(self): + return deepcopy(self._fast_yields) + + +class AveragedFissionYieldHelper(TalliedFissionYieldHelper): + r"""Class that computes fission yields based on average fission energy + + Computes average energy at which fission events occured with + + .. math:: + + \bar{E} = \frac{ + \int_0^\infty E\sigma_f(E)\phi(E)dE + }{ + \int_0^\infty\sigma_f(E)\phi(E)dE + } + + If the average energy for a nuclide is below the lowest energy + with yield data, that set of fission yields is taken. + Conversely, if the average energy is above the highest energy + with yield data, that set of fission yields is used. + For the case where the average energy is between two sets + of yields, the effective fission yield computed by + linearly interpolating between yields provided at the + nearest energies above and below the average. + + Parameters + ---------- + chain_nuclides : iterable of openmc.deplete.Nuclide + Nuclides tracked in the depletion chain. All nuclides are + not required to have fission yield data. + + Attributes + ---------- + constant_yields : dict of str to :class:`openmc.deplete.FissionYield` + Fission yields for all nuclides that only have one set of + fission yield data. Can be accessed as ``{parent: {product: yield}}`` + results : None or numpy.ndarray + If tallies have been generated and unpacked, then the array will + have shape ``(n_mats, n_tnucs)``, where ``n_mats`` is the number + of materials where fission reactions were tallied and ``n_tnucs`` + is the number of nuclides with multiple sets of fission yields. + Data in the array are the average energy of fission events for + tallied nuclides across burnable materials. + """ + + def __init__(self, chain_nuclides): + super().__init__(chain_nuclides) + self._weighted_tally = None + + def generate_tallies(self, materials, mat_indexes): + """Construct tallies to determine average energy of fissions + + Parameters + ---------- + materials : iterable of :class:`openmc.capi.Material` + Materials to be used in :class:`openmc.capi.MaterialFilter` + mat_indexes : iterable of int + Indices of tallied materials that will have their fission + yields computed by this helper. Necessary as the + :class:`openmc.deplete.Operator` that uses this helper + may only burn a subset of all materials when running + in parallel mode. + """ + super().generate_tallies(materials, mat_indexes) + fission_tally = self._fission_rate_tally + filters = fission_tally.filters + + ene_filter = EnergyFilter([0, self._upper_energy]) + fission_tally.filters = filters + [ene_filter] + + func_filter = EnergyFunctionFilter() + func_filter.set_data((0, self._upper_energy), (0, self._upper_energy)) + weighted_tally = Tally() + weighted_tally.scores = ['fission'] + weighted_tally.filters = filters + [func_filter] + self._weighted_tally = weighted_tally + + def update_tally_nuclides(self, nuclides): + """Tally nuclides with non-zero density and multiple yields + + Must be run after :meth:`generate_tallies`. + + Parameters + ---------- + nuclides : iterable of str + Potential nuclides to be tallied, such as those with + non-zero density at this stage. + + Returns + ------- + nuclides : tuple of str + Union of input nuclides and those that have multiple sets + of yield data. Sorted by nuclide name + + Raises + ------ + AttributeError + If tallies not generated + """ + tally_nucs = super().update_tally_nuclides(nuclides) + self._weighted_tally.nuclides = tally_nucs + return tally_nucs + + def unpack(self): + """Unpack tallies and populate :attr:`results` with average energies""" + if not self._tally_nucs: + self.results = None + return + fission_results = ( + self._fission_rate_tally.results[self._local_indexes, :, 1]) + self.results = ( + self._weighted_tally.results[self._local_indexes, :, 1]).copy() + nz_mat, nz_nuc = fission_results.nonzero() + self.results[nz_mat, nz_nuc] /= fission_results[nz_mat, nz_nuc] + + def weighted_yields(self, local_mat_index): + """Return fission yields for a specific material + + Use the computed average energy of fission + events to determine fission yields. If average + energy is between two sets of yields, linearly + interpolate bewteen the two. + Otherwise take the closet set of yields. + + Parameters + ---------- + local_mat_index : int + Index for specific burnable material. Effective + yields will be produced using + ``self.results[local_mat_index]`` + + Returns + ------- + library : dict + Dictionary of ``{parent: {product: fyield}}`` + """ + if not self._tally_nucs: + return self.constant_yields + mat_yields = {} + average_energies = self.results[local_mat_index] + for avg_e, nuc in zip(average_energies, self._tally_nucs): + nuc_energies = nuc.yield_energies + if avg_e <= nuc_energies[0]: + mat_yields[nuc.name] = nuc.yield_data[nuc_energies[0]] + continue + if avg_e >= nuc_energies[-1]: + mat_yields[nuc.name] = nuc.yield_data[nuc_energies[-1]] + continue + # in-between two energies + # linear search since there are usually ~3 energies + for ix, ene in enumerate(nuc_energies[:-1]): + if nuc_energies[ix + 1] > avg_e: + break + lower, upper = nuc_energies[ix:ix + 2] + fast_frac = (avg_e - lower) / (upper - lower) + mat_yields[nuc.name] = ( + nuc.yield_data[lower] * (1 - fast_frac) + + nuc.yield_data[upper] * fast_frac) + mat_yields.update(self.constant_yields) + return mat_yields + + @classmethod + def from_operator(cls, operator, **kwargs): + """Return a new helper with data from an operator + + All keyword arguments should be identical to their counterpart + in the main ``__init__`` method + + Parameters + ---------- + operator : openmc.deplete.TransportOperator + Operator with a depletion chain + kwargs : + Additional keyword arguments to be used in construction + + Returns + ------- + AveragedFissionYieldHelper + """ + return cls(operator.chain.nuclides) diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 538b697904..8896180a88 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -3,13 +3,20 @@ Contains the per-nuclide components of a depletion chain. """ +import bisect +from collections.abc import Mapping from collections import namedtuple, defaultdict from warnings import warn +from numbers import Real try: import lxml.etree as ET except ImportError: import xml.etree.ElementTree as ET +from numpy import empty + +from openmc.checkvalue import check_type + DecayTuple = namedtuple('DecayTuple', 'type target branching_ratio') DecayTuple.__doc__ = """\ @@ -62,11 +69,16 @@ except AttributeError: class Nuclide(object): """Decay modes, reactions, and fission yields for a single nuclide. + Parameters + ---------- + name : str, optional + GND name of this nuclide, e.g. ``"He4"``, ``"Am242_m1"`` + Attributes ---------- - name : str + name : str or None Name of nuclide. - half_life : float + half_life : float or None Half life of nuclide in [s]. decay_energy : float Energy deposited from decay in [eV]. @@ -80,17 +92,16 @@ class Nuclide(object): reactions : list of openmc.deplete.ReactionTuple Reaction information. Each element of the list is a named tuple with attribute 'type', 'target', 'Q', and 'branching_ratio'. - yield_data : dict of float to list - Maps tabulated energy to list of (product, yield) for all - neutron-induced fission products. - yield_energies : list of float - Energies at which fission product yiels exist - + yield_data : FissionYieldDistribution or None + Fission product yields at tabulated energies for this nuclide. Can be + treated as a nested dictionary ``{energy: {product: yield}}`` + yield_energies : tuple of float or None + Energies at which fission product yields exist """ - def __init__(self): + def __init__(self, name=None): # Information about the nuclide - self.name = None + self.name = name self.half_life = None self.decay_energy = 0.0 @@ -101,8 +112,7 @@ class Nuclide(object): self.reactions = [] # Neutron fission yields, if present - self.yield_data = {} - self.yield_energies = [] + self._yield_data = None @property def n_decay_modes(self): @@ -112,6 +122,29 @@ class Nuclide(object): def n_reaction_paths(self): return len(self.reactions) + @property + def yield_data(self): + if self._yield_data is None: + return None + return self._yield_data + + @yield_data.setter + def yield_data(self, fission_yields): + if fission_yields is None: + self._yield_data = None + else: + check_type("fission_yields", fission_yields, Mapping) + if isinstance(fission_yields, FissionYieldDistribution): + self._yield_data = fission_yields + else: + self._yield_data = FissionYieldDistribution(fission_yields) + + @property + def yield_energies(self): + if self._yield_data is None: + return None + return self.yield_data.energies + @classmethod def from_xml(cls, element, fission_q=None): """Read nuclide from an XML element. @@ -166,13 +199,7 @@ class Nuclide(object): fpy_elem = element.find('neutron_fission_yields') if fpy_elem is not None: - for yields_elem in fpy_elem.iter('fission_yields'): - E = float(yields_elem.get('energy')) - products = yields_elem.find('products').text.split() - yields = [float(y) for y in - yields_elem.find('data').text.split()] - nuc.yield_data[E] = list(zip(products, yields)) - nuc.yield_energies = list(sorted(nuc.yield_data.keys())) + nuc.yield_data = FissionYieldDistribution.from_xml_element(fpy_elem) return nuc @@ -212,15 +239,7 @@ class Nuclide(object): fpy_elem = ET.SubElement(elem, 'neutron_fission_yields') energy_elem = ET.SubElement(fpy_elem, 'energies') energy_elem.text = ' '.join(str(E) for E in self.yield_energies) - - for E in self.yield_energies: - yields_elem = ET.SubElement(fpy_elem, 'fission_yields') - yields_elem.set('energy', str(E)) - - products_elem = ET.SubElement(yields_elem, 'products') - products_elem.text = ' '.join(x[0] for x in self.yield_data[E]) - data_elem = ET.SubElement(yields_elem, 'data') - data_elem.text = ' '.join(str(x[1]) for x in self.yield_data[E]) + self.yield_data.to_xml_element(fpy_elem) return elem @@ -304,8 +323,8 @@ class Nuclide(object): valid = False if self.yield_data: - for energy, yield_list in self.yield_data.items(): - sum_yield = sum(y[1] for y in yield_list) + for energy, fission_yield in self.yield_data.items(): + sum_yield = fission_yield.yields.sum() stat = 2.0 - tolerance <= sum_yield <= 2.0 + tolerance if stat: continue @@ -321,3 +340,227 @@ class Nuclide(object): valid = False return valid + + +class FissionYieldDistribution(Mapping): + """Energy-dependent fission product yields for a single nuclide + + Can be used as a dictionary mapping energies and products to fission + yields:: + + >>> fydist = FissionYieldDistribution{ + ... {0.0253: {"Xe135": 0.021}}) + >>> fydist[0.0253]["Xe135"] + 0.021 + + Parameters + ---------- + fission_yields : dict + Dictionary of energies and fission product yields for that energy. + Expected to be of the form ``{float: {str: float}}``. The first + float is the energy, typically in eV, that represents this + distribution. The underlying dictionary maps fission products + to their respective yields. + + Attributes + ---------- + energies : tuple + Energies for which fission yields exist. Sorted by + increasing energy + products : tuple + Fission products produced at all energies. Sorted by name. + yield_matrix : numpy.ndarray + Array ``(n_energy, n_products)`` where + ``yield_matrix[g, j]`` is the fission yield of product + ``j`` for energy group ``g``. + + See Also + -------- + * :meth:`from_xml_element` - Construction methods + * :class:`FissionYield` - Class used for storing yields at a given energy + """ + + def __init__(self, fission_yields): + # mapping {energy: {product: value}} + energies = sorted(fission_yields) + + # Get a consistent set of products to produce a matrix of yields + shared_prod = set.union(*(set(x) for x in fission_yields.values())) + ordered_prod = sorted(shared_prod) + + yield_matrix = empty((len(energies), len(shared_prod))) + + for g_index, energy in enumerate(energies): + prod_map = fission_yields[energy] + for prod_ix, product in enumerate(ordered_prod): + yield_val = prod_map.get(product) + yield_matrix[g_index, prod_ix] = ( + 0.0 if yield_val is None else yield_val) + self.energies = tuple(energies) + self.products = tuple(ordered_prod) + self.yield_matrix = yield_matrix + + def __len__(self): + return len(self.energies) + + def __getitem__(self, energy): + if energy not in self.energies: + raise KeyError(energy) + return FissionYield( + self.products, self.yield_matrix[self.energies.index(energy)]) + + def __iter__(self): + return iter(self.energies) + + def __repr__(self): + return "<{} with {} products at {} energies>".format( + self.__class__.__name__, self.yield_matrix.shape[1], + len(self.energies)) + + @classmethod + def from_xml_element(cls, element): + """Construct a distribution from a depletion chain xml file + + Parameters + ---------- + element : xml.etree.ElementTree.Element + XML element to pull fission yield data from + + Returns + ------- + FissionYieldDistribution + """ + all_yields = {} + for elem_index, yield_elem in enumerate(element.iter("fission_yields")): + energy = float(yield_elem.get("energy")) + products = yield_elem.find("products").text.split() + yields = map(float, yield_elem.find("data").text.split()) + # Get a map of products to their corresponding yield + all_yields[energy] = dict(zip(products, yields)) + + return cls(all_yields) + + def to_xml_element(self, root): + """Write fission yield data to an xml element + + Parameters + ---------- + root : xml.etree.ElementTree.Element + Element to write distribution data to + """ + for energy, yield_obj in self.items(): + yield_element = ET.SubElement(root, "fission_yields") + yield_element.set("energy", str(energy)) + product_elem = ET.SubElement(yield_element, "products") + product_elem.text = " ".join(map(str, yield_obj.products)) + data_elem = ET.SubElement(yield_element, "data") + data_elem.text = " ".join(map(str, yield_obj.yields)) + + +class FissionYield(Mapping): + """Mapping for fission yields of a parent at a specific energy + + Separated to support nested dictionary-like behavior for + :class:`FissionYieldDistribution`, and allowing math operations + on a single vector of yields. Can in turn be used like a + dictionary to fetch fission yields. + Supports multiplication of a scalar to scale the fission + yields and addition of another set of yields. + + Does not support resizing / inserting new products that do + not exist. + + Parameters + ---------- + products : tuple of str + Products for this specific distribution + yields : numpy.ndarray + Fission product yields for each product in ``products`` + + Attributes + ---------- + products : tuple of str + Products for this specific distribution + yields : numpy.ndarray + Fission product yields for each product in ``products`` + + Examples + -------- + >>> import numpy + >>> fy_vector = FissionYield( + ... ("Xe135", "I129", "Sm149"), + ... numpy.array((0.002, 0.001, 0.0003))) + >>> fy_vector["Xe135"] + 0.002 + >>> new = fy_vector.copy() + >>> fy_vector *= 2 + >>> fy_vector["Xe135"] + 0.004 + >>> new["Xe135"] + 0.002 + >>> (new + fy_vector)["Sm149"] + 0.0009 + >>> dict(new) + {"Xe135": 0.002, "I129": 0.001, "Sm149": 0.0003} + """ + + def __init__(self, products, yields): + self.products = products + self.yields = yields + + def __contains__(self, product): + ix = bisect.bisect_left(self.products, product) + return ix != len(self.products) and self.products[ix] == product + + def __getitem__(self, product): + ix = bisect.bisect_left(self.products, product) + if ix == len(self.products) or self.products[ix] != product: + raise KeyError(product) + return self.yields[ix] + + def __len__(self): + return len(self.products) + + def __iter__(self): + return iter(self.products) + + def items(self): + """Return pairs of product, yield""" + return zip(self.products, self.yields) + + def __add__(self, other): + if not isinstance(other, FissionYield): + return NotImplemented + new = FissionYield(self.products, self.yields.copy()) + new += other + return new + + def __iadd__(self, other): + """Increment value from other fission yield""" + if not isinstance(other, FissionYield): + return NotImplemented + self.yields += other.yields + return self + + def __radd__(self, other): + return self + other + + def __imul__(self, scalar): + if not isinstance(scalar, Real): + return NotImplemented + self.yields *= scalar + return self + + def __mul__(self, scalar): + if not isinstance(scalar, Real): + return NotImplemented + new = FissionYield(self.products, self.yields.copy()) + new *= scalar + return new + + def __rmul__(self, scalar): + return self * scalar + + def __repr__(self): + return "<{} containing {} products and yields>".format( + self.__class__.__name__, len(self)) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index f8288841f7..01bf92dcc9 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -25,7 +25,9 @@ from .abc import TransportOperator, OperatorResult from .atom_number import AtomNumber from .reaction_rates import ReactionRates from .results_list import ResultsList -from .helpers import DirectReactionRateHelper, ChainFissionHelper +from .helpers import ( + DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, + FissionYieldCutoffHelper, AveragedFissionYieldHelper) def _distribute(items): @@ -84,6 +86,21 @@ class Operator(TransportOperator): in initial condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. Defaults to 1.0e3. + fission_yield_mode : {"constant", "cutoff", "average"} + Key indicating what fission product yield scheme to use. The + key determines what fission energy helper is used: + + * "constant": :class:`~openmc.deplete.helpers.ConstantFissionYieldHelper` + * "cutoff": :class:`~openmc.deplete.helpers.FissionYieldCutoffHelper` + * "average": :class:`~openmc.deplete.helpers.AveragedFissionYieldHelper` + + The documentation on these classes describe their methodology + and differences. Default: ``"constant"`` + fission_yield_opts : dict of str to option, optional + Optional arguments to pass to the helper determined by + ``fission_yield_mode``. Will be passed directly on to the + helper. Passing a value of None will use the defaults for + the associated helper. Attributes ---------- @@ -120,9 +137,20 @@ class Operator(TransportOperator): diff_burnable_mats : bool Whether to differentiate burnable materials with multiple instances """ + _fission_helpers = { + "average": AveragedFissionYieldHelper, + "constant": ConstantFissionYieldHelper, + "cutoff": FissionYieldCutoffHelper, + } + def __init__(self, geometry, settings, chain_file=None, prev_results=None, diff_burnable_mats=False, fission_q=None, - dilute_initial=1.0e3): + dilute_initial=1.0e3, fission_yield_mode="constant", + fission_yield_opts=None): + if fission_yield_mode not in self._fission_helpers: + raise KeyError( + "fission_yield_mode must be one of {}, not {}".format( + ", ".join(self._fission_helpers), fission_yield_mode)) super().__init__(chain_file, fission_q, dilute_initial, prev_results) self.round_number = False self.prev_res = None @@ -178,6 +206,13 @@ class Operator(TransportOperator): self.reaction_rates.n_nuc, self.reaction_rates.n_react) self._energy_helper = ChainFissionHelper() + # Select and create fission yield helper + fission_helper = self._fission_helpers[fission_yield_mode] + fission_yield_opts = ( + {} if fission_yield_opts is None else fission_yield_opts) + self._yield_helper = fission_helper.from_operator( + self, **fission_yield_opts) + def __call__(self, vec, power): """Runs a simulation. @@ -204,8 +239,10 @@ class Operator(TransportOperator): # Update material compositions and tally nuclides self._update_materials() - self._rate_helper.nuclides = self._get_tally_nuclides() - self._energy_helper.nuclides = self._rate_helper.nuclides + nuclides = self._get_tally_nuclides() + self._rate_helper.nuclides = nuclides + self._energy_helper.nuclides = nuclides + self._yield_helper.update_tally_nuclides(nuclides) # Run OpenMC openmc.capi.reset() @@ -409,6 +446,10 @@ class Operator(TransportOperator): self._rate_helper.generate_tallies(materials, self.chain.reactions) self._energy_helper.prepare( self.chain.nuclides, self.reaction_rates.index_nuc, materials) + # Tell fission yield helper what materials this process is + # responsible for + self._yield_helper.generate_tallies( + materials, tuple(sorted(self._mat_index_map.values()))) # Return number density vector return list(self.number.get_mat_slice(np.s_[:])) @@ -554,6 +595,10 @@ class Operator(TransportOperator): # Keep track of energy produced from all reactions in eV per source # particle self._energy_helper.reset() + self._yield_helper.unpack() + + # Store fission yield dictionaries + fission_yields = [] # Create arrays to store fission Q values, reaction rates, and nuclide # numbers, zeroed out in material iteration @@ -576,6 +621,9 @@ class Operator(TransportOperator): tally_rates = self._rate_helper.get_material_rates( mat_index, nuc_ind, react_ind) + # Compute fission yields for this material + fission_yields.append(self._yield_helper.weighted_yields(i)) + # Accumulate energy from fission self._energy_helper.update(tally_rates[:, fission_ind], mat_index) @@ -589,6 +637,9 @@ class Operator(TransportOperator): # Scale reaction rates to obtain units of reactions/sec rates *= power / energy + # Store new fission yields on the chain + self.chain.fission_yields = fission_yields + return OperatorResult(k_combined, rates) def _get_nuclides_with_data(self): diff --git a/tests/dummy_operator.py b/tests/dummy_operator.py index 37269b7bd6..85c689743a 100644 --- a/tests/dummy_operator.py +++ b/tests/dummy_operator.py @@ -73,6 +73,49 @@ SCHEMES = { } +class TestChain(object): + """Empty chain to assist with unit testing depletion routines + + Only really provides the form_matrix function, but acts like + a real Chain + """ + + fission_yields = [None] + + @staticmethod + def get_thermal_fission_yields(): + return None + + def form_matrix(self, rates, _fission_yields=None): + """Forms the f(y) matrix in y' = f(y)y. + + Nominally a depletion matrix, this is abstracted on the off chance + that the function f has nothing to do with depletion at all. + + Parameters + ---------- + rates : numpy.ndarray + Slice of reaction rates for a single material + _fission_yields : optional + Not used + + Returns + ------- + scipy.sparse.csr_matrix + Sparse matrix representing f(y). + """ + + y_1 = rates[0, 0] + y_2 = rates[1, 0] + + a11 = np.sin(y_2) + a12 = np.cos(y_1) + a21 = -np.cos(y_2) + a22 = np.sin(y_1) + + return sp.csr_matrix(np.array([[a11, a12], [a21, a22]])) + + class DummyOperator(TransportOperator): """This is a dummy operator class with no statistical uncertainty. @@ -88,6 +131,7 @@ class DummyOperator(TransportOperator): """ def __init__(self, previous_results=None): self.prev_res = previous_results + self.chain = TestChain() self.output_dir = "." def __call__(self, vec, power, print_out=False): @@ -120,37 +164,6 @@ class DummyOperator(TransportOperator): # Create a fake rates object return OperatorResult(ufloat(0.0, 0.0), reaction_rates) - @property - def chain(self): - return self - - def form_matrix(self, rates): - """Forms the f(y) matrix in y' = f(y)y. - - Nominally a depletion matrix, this is abstracted on the off chance - that the function f has nothing to do with depletion at all. - - Parameters - ---------- - rates : numpy.ndarray - Slice of reaction rates for a single material - - Returns - ------- - scipy.sparse.csr_matrix - Sparse matrix representing f(y). - """ - - y_1 = rates[0, 0] - y_2 = rates[1, 0] - - a11 = np.sin(y_2) - a12 = np.cos(y_1) - a21 = -np.cos(y_2) - a22 = np.sin(y_1) - - return sp.csr_matrix(np.array([[a11, a12], [a21, a22]])) - @property def volume(self): """ diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 1fd974e054..3dcd3ad750 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -3,10 +3,11 @@ from collections.abc import Mapping import os from pathlib import Path +from itertools import product import numpy as np from openmc.data import zam, ATOMIC_SYMBOL -from openmc.deplete import comm, Chain, reaction_rates, nuclide +from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram import pytest from tests import cdtemp @@ -128,9 +129,10 @@ def test_from_xml(simple_chain): assert [r.branching_ratio for r in nuc.reactions] == [1.0, 0.7, 0.3] # Yield tests - assert nuc.yield_energies == [0.0253] + assert nuc.yield_energies == (0.0253,) assert list(nuc.yield_data) == [0.0253] - assert nuc.yield_data[0.0253] == [("A", 0.0292737), ("B", 0.002566345)] + assert nuc.yield_data[0.0253].products == ("A", "B") + assert (nuc.yield_data[0.0253].yields == [0.0292737, 0.002566345]).all() def test_export_to_xml(run_in_tmpdir): @@ -139,8 +141,7 @@ def test_export_to_xml(run_in_tmpdir): # Prevent different MPI ranks from conflicting filename = 'test{}.xml'.format(comm.rank) - A = nuclide.Nuclide() - A.name = "A" + A = nuclide.Nuclide("A") A.half_life = 2.36520e4 A.decay_modes = [ nuclide.DecayTuple("beta1", "B", 0.6), @@ -148,21 +149,19 @@ def test_export_to_xml(run_in_tmpdir): ] A.reactions = [nuclide.ReactionTuple("(n,gamma)", "C", 0.0, 1.0)] - B = nuclide.Nuclide() - B.name = "B" + B = nuclide.Nuclide("B") B.half_life = 3.29040e4 B.decay_modes = [nuclide.DecayTuple("beta", "A", 1.0)] B.reactions = [nuclide.ReactionTuple("(n,gamma)", "C", 0.0, 1.0)] - C = nuclide.Nuclide() - C.name = "C" + C = nuclide.Nuclide("C") C.reactions = [ nuclide.ReactionTuple("fission", None, 2.0e8, 1.0), nuclide.ReactionTuple("(n,gamma)", "A", 0.0, 0.7), nuclide.ReactionTuple("(n,gamma)", "B", 0.0, 0.3) ] - C.yield_energies = [0.0253] - C.yield_data = {0.0253: [("A", 0.0292737), ("B", 0.002566345)]} + C.yield_data = nuclide.FissionYieldDistribution({ + 0.0253: {"A": 0.0292737, "B": 0.002566345}}) chain = Chain() chain.nuclides = [A, B, C] @@ -220,6 +219,13 @@ def test_form_matrix(simple_chain): assert mat[1, 2] == mat12 assert mat[2, 2] == mat22 + # Pass equivalent fission yields directly + # Ensure identical matrix is formed + f_yields = {"C": {"A": 0.0292737, "B": 0.002566345}} + new_mat = chain.form_matrix(react[0], f_yields) + for r, c in product(range(3), range(3)): + assert new_mat[r, c] == mat[r, c] + def test_getitem(): """ Test nuc_by_ind converter function. """ @@ -283,8 +289,7 @@ def test_capture_branch_infer_ground(): chain = Chain.from_xml(chain_file) # Create nuclide to be added into the chain - xe136m = nuclide.Nuclide() - xe136m.name = "Xe136_m1" + xe136m = nuclide.Nuclide("Xe136_m1") chain.nuclides.append(xe136m) chain.nuclide_dict[xe136m.name] = len(chain.nuclides) - 1 @@ -301,8 +306,7 @@ def test_capture_branch_no_rxn(): chain_file = Path(__file__).parents[1] / "chain_simple.xml" chain = Chain.from_xml(chain_file) - u5m = nuclide.Nuclide() - u5m.name = "U235_m1" + u5m = nuclide.Nuclide("U235_m1") chain.nuclides.append(u5m) chain.nuclide_dict[u5m.name] = len(chain.nuclides) - 1 @@ -378,6 +382,32 @@ def test_set_alpha_branches(): raise ValueError("Helium has been removed and should not have been") +def test_simple_fission_yields(simple_chain): + """Check the default fission yields that can be used to form the matrix + """ + fission_yields = simple_chain.get_thermal_fission_yields() + assert fission_yields == {"C": {"A": 0.0292737, "B": 0.002566345}} + + +def test_fission_yield_attribute(simple_chain): + """Test the fission_yields property""" + thermal_yields = simple_chain.get_thermal_fission_yields() + # generate default with property + assert simple_chain.fission_yields[0] == thermal_yields + empty_chain = Chain() + empty_chain.fission_yields = thermal_yields + assert empty_chain.fission_yields[0] == thermal_yields + empty_chain.fission_yields = [thermal_yields] * 2 + assert empty_chain.fission_yields[0] == thermal_yields + assert empty_chain.fission_yields[1] == thermal_yields + + # test failure with deplete function + # number fission yields != number of materials + 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) + def test_validate(simple_chain): """Test the validate method""" @@ -394,7 +424,7 @@ def test_validate(simple_chain): # Fix fission yields but keep to restore later old_yields = simple_chain["C"].yield_data - simple_chain["C"].yield_data = {0.0253: [("A", 1.4), ("B", 0.6)]} + simple_chain["C"].yield_data = {0.0253: {"A": 1.4, "B": 0.6}} assert simple_chain.validate(strict=True, tolerance=0.0) with pytest.warns(None) as record: diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py new file mode 100644 index 0000000000..bedd702a90 --- /dev/null +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -0,0 +1,293 @@ +"""Test the FissionYieldHelpers""" + +import os +from collections import namedtuple +from unittest.mock import Mock +import bisect + +import pytest +import numpy as np +import openmc +from openmc import capi +from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution +from openmc.deplete.helpers import ( + FissionYieldCutoffHelper, ConstantFissionYieldHelper, + AveragedFissionYieldHelper) + + +@pytest.fixture(scope="module") +def materials(tmpdir_factory): + """Use C API to construct realistic materials for testing tallies""" + tmpdir = tmpdir_factory.mktemp("capi") + orig = tmpdir.chdir() + # Create proxy problem to please openmc + mfuel = openmc.Material(name="test_fuel") + mfuel.volume = 1.0 + for nuclide in ["U235", "U238", "Xe135", "Pu239"]: + mfuel.add_nuclide(nuclide, 1.0) + openmc.Materials([mfuel]).export_to_xml() + # Geometry + box = openmc.rectangular_prism(1.0, 1.0, boundary_type="reflective") + cell = openmc.Cell(fill=mfuel, region=box) + root = openmc.Universe(cells=[cell]) + openmc.Geometry(root).export_to_xml() + # settings + settings = openmc.Settings() + settings.particles = 100 + settings.inactive = 0 + settings.batches = 10 + settings.verbosity = 1 + settings.export_to_xml() + + try: + with capi.run_in_memory(): + yield [capi.Material(), capi.Material()] + finally: + # Convert to strings as os.remove in py 3.5 doesn't support Paths + for file_path in ("settings.xml", "geometry.xml", "materials.xml", + "summary.h5"): + os.remove(str(tmpdir / file_path)) + orig.chdir() + os.rmdir(str(tmpdir)) + + +def proxy_tally_data(tally, fill=None): + """Construct an empty matrix built from a C tally + + The shape of tally.results will be + ``(n_bins, n_nuc * n_scores, 3)`` + """ + n_nucs = max(len(tally.nuclides), 1) + n_scores = max(len(tally.scores), 1) + n_bins = 1 + for tfilter in tally.filters: + if not hasattr(tfilter, "bins"): + continue + this_bins = len(tfilter.bins) + if isinstance(tfilter, capi.EnergyFilter): + this_bins -= 1 + n_bins *= max(this_bins, 1) + data = np.empty((n_bins, n_nucs * n_scores, 3)) + if fill is not None: + data.fill(fill) + return data + + +@pytest.fixture(scope="module") +def nuclide_bundle(): + u5yield_dict = { + 0.0253: {"Xe135": 7.85e-4, "Gd155": 4.08e-12, "Sm149": 1.71e-12}, + 5.0e5: {"Xe135": 7.85e-4, "Sm149": 1.71e-12}, + 1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8}} + u235 = Nuclide("U235") + u235.yield_data = FissionYieldDistribution(u5yield_dict) + + u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}} + u238 = Nuclide("U238") + u238.yield_data = FissionYieldDistribution(u8yield_dict) + + xe135 = Nuclide("Xe135") + + pu239 = Nuclide("Pu239") + pu239.yield_data = FissionYieldDistribution({ + 5.0e5: {"Xe135": 6.14e-3, "Sm149": 9.429e-10, "Gd155": 5.24e-9}, + 2e6: {"Xe135": 6.15e-3, "Sm149": 9.42e-10, "Gd155": 5.29e-9}}) + + NuclideBundle = namedtuple("NuclideBundle", "u235 u238 xe135 pu239") + return NuclideBundle(u235, u238, xe135, pu239) + + +@pytest.mark.parametrize( + "input_energy, yield_energy", + ((0.0253, 0.0253), (0.01, 0.0253), (4e5, 5e5))) +def test_constant_helper(nuclide_bundle, input_energy, yield_energy): + helper = ConstantFissionYieldHelper(nuclide_bundle, energy=input_energy) + assert helper.energy == input_energy + assert helper.constant_yields == { + "U235": nuclide_bundle.u235.yield_data[yield_energy], + "U238": nuclide_bundle.u238.yield_data[5.00e5], + "Pu239": nuclide_bundle.pu239.yield_data[5e5]} + assert helper.constant_yields == helper.weighted_yields(1) + + +def test_cutoff_construction(nuclide_bundle): + u235 = nuclide_bundle.u235 + u238 = nuclide_bundle.u238 + pu239 = nuclide_bundle.pu239 + + # defaults + helper = FissionYieldCutoffHelper(nuclide_bundle, 1) + assert helper.constant_yields == { + "U238": u238.yield_data[5.0e5], + "Pu239": pu239.yield_data[5e5]} + assert helper.thermal_yields == {"U235": u235.yield_data[0.0253]} + assert helper.fast_yields == {"U235": u235.yield_data[5e5]} + + # use 14 MeV yields + helper = FissionYieldCutoffHelper(nuclide_bundle, 1, fast_energy=14e6) + assert helper.constant_yields == { + "U238": u238.yield_data[5.0e5], + "Pu239": pu239.yield_data[5e5]} + assert helper.thermal_yields == {"U235": u235.yield_data[0.0253]} + assert helper.fast_yields == {"U235": u235.yield_data[14e6]} + + # specify missing thermal yields -> use 0.0253 + helper = FissionYieldCutoffHelper(nuclide_bundle, 1, thermal_energy=1) + assert helper.thermal_yields == {"U235": u235.yield_data[0.0253]} + assert helper.fast_yields == {"U235": u235.yield_data[5e5]} + + # request missing fast yields -> use epithermal + helper = FissionYieldCutoffHelper(nuclide_bundle, 1, fast_energy=1e4) + assert helper.thermal_yields == {"U235": u235.yield_data[0.0253]} + assert helper.fast_yields == {"U235": u235.yield_data[5e5]} + + # higher cutoff energy -> obtain fast and "faster" yields + helper = FissionYieldCutoffHelper(nuclide_bundle, 1, cutoff=1e6, + thermal_energy=5e5, fast_energy=14e6) + assert helper.constant_yields == {"U238": u238.yield_data[5e5]} + assert helper.thermal_yields == { + "U235": u235.yield_data[5e5], "Pu239": pu239.yield_data[5e5]} + assert helper.fast_yields == { + "U235": u235.yield_data[14e6], "Pu239": pu239.yield_data[2e6]} + + # test super low and super high cutoff energies + helper = FissionYieldCutoffHelper( + nuclide_bundle, 1, thermal_energy=0.001, cutoff=0.002) + assert helper.fast_yields == {} + assert helper.thermal_yields == {} + assert helper.constant_yields == { + "U235": u235.yield_data[0.0253], "U238": u238.yield_data[5e5], + "Pu239": pu239.yield_data[5e5]} + + helper = FissionYieldCutoffHelper( + nuclide_bundle, 1, cutoff=15e6, fast_energy=17e6) + assert helper.thermal_yields == {} + assert helper.fast_yields == {} + assert helper.constant_yields == { + "U235": u235.yield_data[14e6], "U238": u238.yield_data[5e5], + "Pu239": pu239.yield_data[2e6]} + + +@pytest.mark.parametrize("key", ("cutoff", "thermal_energy", "fast_energy")) +def test_cutoff_failure(key): + with pytest.raises(TypeError, match=key): + FissionYieldCutoffHelper(None, None, **{key: None}) + with pytest.raises(ValueError, match=key): + FissionYieldCutoffHelper(None, None, **{key: -1}) + + +# emulate some split between fast and thermal U235 fissions +@pytest.mark.parametrize("therm_frac", (0.5, 0.2, 0.8)) +def test_cutoff_helper(materials, nuclide_bundle, therm_frac): + helper = FissionYieldCutoffHelper(nuclide_bundle, len(materials), + cutoff=1e6, fast_energy=14e6) + helper.generate_tallies(materials, [0]) + + non_zero_nucs = [n.name for n in nuclide_bundle] + tally_nucs = helper.update_tally_nuclides(non_zero_nucs) + assert tally_nucs == ["Pu239", "U235"] + + # Check tallies + fission_tally = helper._fission_rate_tally + assert fission_tally is not None + filters = fission_tally.filters + assert len(filters) == 2 + assert isinstance(filters[0], capi.MaterialFilter) + assert len(filters[0].bins) == len(materials) + assert isinstance(filters[1], capi.EnergyFilter) + # lower, cutoff, and upper energy + assert len(filters[1].bins) == 3 + + # Emulate building tallies + # material x energy, tallied_nuclides, 3 + tally_data = proxy_tally_data(fission_tally) + helper._fission_rate_tally = Mock() + helper_flux = 1e6 + tally_data[0, :, 1] = therm_frac * helper_flux + tally_data[1, :, 1] = (1 - therm_frac) * helper_flux + helper._fission_rate_tally.results = tally_data + + helper.unpack() + # expected results of shape (n_mats, 2, n_tnucs) + expected_results = np.empty((1, 2, len(tally_nucs))) + expected_results[:, 0] = therm_frac + expected_results[:, 1] = 1 - therm_frac + assert helper.results == pytest.approx(expected_results) + + actual_yields = helper.weighted_yields(0) + assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5] + for nuc in tally_nucs: + assert actual_yields[nuc] == ( + helper.thermal_yields[nuc] * therm_frac + + helper.fast_yields[nuc] * (1 - therm_frac)) + + +@pytest.mark.parametrize("avg_energy", (0.01, 6e5, 15e6)) +def test_averaged_helper(materials, nuclide_bundle, avg_energy): + helper = AveragedFissionYieldHelper(nuclide_bundle) + helper.generate_tallies(materials, [0]) + tallied_nucs = helper.update_tally_nuclides( + [n.name for n in nuclide_bundle]) + assert tallied_nucs == ["Pu239", "U235"] + + # check generated tallies + fission_tally = helper._fission_rate_tally + assert fission_tally is not None + fission_filters = fission_tally.filters + assert len(fission_filters) == 2 + assert isinstance(fission_filters[0], capi.MaterialFilter) + assert len(fission_filters[0].bins) == len(materials) + assert isinstance(fission_filters[1], capi.EnergyFilter) + assert len(fission_filters[1].bins) == 2 + assert fission_tally.scores == ["fission"] + assert fission_tally.nuclides == list(tallied_nucs) + + weighted_tally = helper._weighted_tally + assert weighted_tally is not None + weighted_filters = weighted_tally.filters + assert len(weighted_filters) == 2 + assert isinstance(weighted_filters[0], capi.MaterialFilter) + assert len(weighted_filters[0].bins) == len(materials) + assert isinstance(weighted_filters[1], capi.EnergyFunctionFilter) + assert len(weighted_filters[1].energy) == 2 + assert len(weighted_filters[1].y) == 2 + assert weighted_tally.scores == ["fission"] + assert weighted_tally.nuclides == list(tallied_nucs) + + helper_flux = 1e16 + fission_results = proxy_tally_data(fission_tally, helper_flux) + weighted_results = proxy_tally_data( + weighted_tally, helper_flux * avg_energy) + + helper._fission_rate_tally = Mock() + helper._weighted_tally = Mock() + helper._fission_rate_tally.results = fission_results + helper._weighted_tally.results = weighted_results + + helper.unpack() + expected_results = np.ones((1, len(tallied_nucs))) * avg_energy + assert helper.results == pytest.approx(expected_results) + + actual_yields = helper.weighted_yields(0) + # constant U238 => no interpolation + assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5] + # construct expected yields + exp_u235_yields = interp_average_yields(nuclide_bundle.u235, avg_energy) + assert actual_yields["U235"] == exp_u235_yields + exp_pu239_yields = interp_average_yields(nuclide_bundle.pu239, avg_energy) + assert actual_yields["Pu239"] == exp_pu239_yields + + +def interp_average_yields(nuc, avg_energy): + """Construct a set of yields by interpolation between neighbors""" + energies = nuc.yield_energies + yields = nuc.yield_data + if avg_energy < energies[0]: + return yields[energies[0]] + if avg_energy > energies[-1]: + return yields[energies[-1]] + thermal_ix = bisect.bisect_left(energies, avg_energy) + thermal_E, fast_E = energies[thermal_ix - 1:thermal_ix + 1] + assert thermal_E < avg_energy < fast_E + split = (avg_energy - thermal_E)/(fast_E - thermal_E) + return yields[thermal_E]*(1 - split) + yields[fast_E]*split diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index 38848e9c19..7d68b3a3d8 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -2,8 +2,8 @@ import xml.etree.ElementTree as ET +import numpy import pytest - from openmc.deplete import nuclide @@ -71,18 +71,20 @@ def test_from_xml(): nuclide.ReactionTuple('(n,gamma)', 'U236', 6545200.0, 1.0), nuclide.ReactionTuple('fission', None, 193405400.0, 1.0), ] - assert u235.yield_energies == [0.0253] - assert u235.yield_data == { - 0.0253: [('Te134', 0.062155), ('Zr100', 0.0497641), - ('Xe138', 0.0481413)] - } + expected_yield_data = nuclide.FissionYieldDistribution({ + 0.0253: {"Xe138": 0.0481413, "Zr100": 0.0497641, "Te134": 0.062155}}) + assert u235.yield_data == expected_yield_data + # test accessing the yield energies through the FissionYieldDistribution + assert u235.yield_energies == (0.0253,) + assert u235.yield_energies is u235.yield_data.energies + with pytest.raises(AttributeError): # not settable + u235.yield_energies = [0.0253, 5e5] def test_to_xml_element(): """Test writing nuclide data to an XML element.""" - C = nuclide.Nuclide() - C.name = "C" + C = nuclide.Nuclide("C") C.half_life = 0.123 C.decay_modes = [ nuclide.DecayTuple('beta-', 'B', 0.99), @@ -92,8 +94,8 @@ def test_to_xml_element(): nuclide.ReactionTuple('fission', None, 2.0e8, 1.0), nuclide.ReactionTuple('(n,gamma)', 'A', 0.0, 1.0) ] - C.yield_energies = [0.0253] - C.yield_data = {0.0253: [("A", 0.0292737), ("B", 0.002566345)]} + C.yield_data = nuclide.FissionYieldDistribution( + {0.0253: {"A": 0.0292737, "B": 0.002566345}}) element = C.to_xml_element() assert element.get("half_life") == "0.123" @@ -118,6 +120,57 @@ def test_to_xml_element(): assert element.find('neutron_fission_yields') is not None +def test_fission_yield_distribution(): + """Test an energy-dependent yield distribution""" + yield_dict = { + 0.0253: {"Xe135": 7.85e-4, "Gd155": 4.08e-12, "Sm149": 1.71e-12}, + 1.40e7: {"Xe135": 4.54e-3, "Gd155": 5.83e-8, "Sm149": 2.69e-8}, + 5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}, # drop Sm149 + } + yield_dist = nuclide.FissionYieldDistribution(yield_dict) + assert len(yield_dist) == len(yield_dict) + assert yield_dist.energies == tuple(sorted(yield_dict.keys())) + for exp_ene, exp_dist in yield_dict.items(): + act_dist = yield_dict[exp_ene] + for exp_prod, exp_yield in exp_dist.items(): + assert act_dist[exp_prod] == exp_yield + exp_yield = numpy.array([ + [4.08e-12, 1.71e-12, 7.85e-4], + [1.32e-12, 0.0, 1.12e-3], + [5.83e-8, 2.69e-8, 4.54e-3]]) + assert numpy.array_equal(yield_dist.yield_matrix, exp_yield) + + # Test the operations / special methods for fission yield + orig_yields = yield_dist[0.0253] + assert len(orig_yields) == len(yield_dict[0.0253]) + for key, value in yield_dict[0.0253].items(): + assert key in orig_yields + assert orig_yields[key] == value + # __getitem__ return yields as a view into yield matrix + assert orig_yields.yields.base is yield_dist.yield_matrix + + # Fission yield feature uses scaled and incremented + mod_yields = orig_yields * 2 + assert numpy.array_equal(orig_yields.yields * 2, mod_yields.yields) + mod_yields += orig_yields + assert numpy.array_equal(orig_yields.yields * 3, mod_yields.yields) + + # Failure modes for adding, multiplying yields + similar = numpy.empty_like(orig_yields.yields) + with pytest.raises(TypeError): + orig_yields + similar + with pytest.raises(TypeError): + similar + orig_yields + with pytest.raises(TypeError): + orig_yields += similar + with pytest.raises(TypeError): + orig_yields * similar + with pytest.raises(TypeError): + similar * orig_yields + with pytest.raises(TypeError): + orig_yields *= similar + + def test_validate(): nuc = nuclide.Nuclide() @@ -141,8 +194,8 @@ def test_validate(): # fission yields nuc.yield_data = { - 0.0253: [("0", 1.5), ("1", 0.5)], - 1e6: [("0", 1.5), ("1", 0.5)], + 0.0253: {"0": 1.5, "1": 0.5}, + 1e6: {"0": 1.5, "1": 0.5}, } # nuclide is good and should have no warnings raise @@ -176,7 +229,7 @@ def test_validate(): # restore reactions, invalidate fission yields nuc.reactions.append(reaction) - nuc.yield_data[1e6].pop() + nuc.yield_data[1e6].yields *= 2 with pytest.raises(ValueError, match=r"fission yields.*1\.0*e"): nuc.validate(strict=True, quiet=False, tolerance=0.0)