diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 260838439..8e8cf755b 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -78,6 +78,7 @@ data, such as number densities and reaction rates for each material. AtomNumber ChainFissionHelper + ConstantFissionYieldHelper DirectReactionRateHelper OperatorResult ReactionRates diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index ac81d0ac7..44ed9293d 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -2,12 +2,19 @@ Class for normalizing fission energy deposition """ from itertools import product +from numbers import Real from numpy import dot, zeros, newaxis, asarray, empty_like, where +from openmc.checkvalue import check_type, check_greater_than from openmc.capi import ( Tally, MaterialFilter, EnergyFilter) -from .abc import ReactionRateHelper, EnergyHelper +from .abc import ( + ReactionRateHelper, EnergyHelper, FissionYieldHelper) + +__all__ = ( + "DirectReactionRateHelper", "ChainFissionHelper", + "ConstantFissionYieldHelper") # ------------------------------------- # Helpers for generating reaction rates @@ -150,152 +157,68 @@ class ChainFissionHelper(EnergyHelper): # ------------------------------------ -class FissionYieldHelper(object): - """Class for using energy-dependent fission yields in depletion chain - - Creates a tally across all burnable materials to score the fission - rate in nuclides with yield data. An energy filter is used to - compute this rates in a group structure corresponding to the - fission yield data. This tally data is used to compute the - relative number of fission events in each energy region, - which serve as the weights for each energy-dependent fission - yield distribution. +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. Not necessary that all have yield data. - n_bmats : int - Number of burnable materials tracked in the problem + 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 ---------- - energy_bounds : tuple of float - Sorted energy bounds from the tally filter - results : numpy.ndarray - Array of tally results for this process with shape - ``(n_local_mat, n_energy, n_nucs)`` + 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, n_bmats): - self._chain_nuclides = {} - self._chain_set = set() - self._tally_map = {} - # TODO Support user-requested minimum energy? - yield_energies = {0.0} - - # Get nuclides with fission yield data, names - # and all energy points - # Names are provided from operator tally nuclides - for nuc in chain_nuclides: - if len(nuc.yield_data) == 0: + 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 - self._chain_nuclides[nuc.name] = nuc - self._chain_set.add(nuc.name) - yield_energies.update(nuc.yield_energies) - - # Create energy grid - self._energy_bounds = tuple(sorted(yield_energies)) - self.n_bmats = n_bmats - - self._reaction_tally = None - self.results = None - self.local_indexes = None + # Specific energy not found, use closest energy + distances = [abs(energy - ene) for ene in nuc.yield_energies] + min_index = min( + range(len(nuc.yield_energies)), key=distances.__getitem__) + self._constant_yields[name] = ( + nuc.yield_data[nuc.yield_energies[min_index]]) @property - def energy_bounds(self): - return self._energy_bounds + def energy(self): + return self._energy - def generate_tallies(self, materials, mat_indexes): - """Construct the fission rate tally + def weighted_yields(self, _local_mat_index=None): + """Return fission yields for all nuclides requested Parameters ---------- - materials : iterable of C-API materials - Materials to be used in :class:`openmc.capi.MaterialFilter` - mat_indexes : iterable of int - Indexes for materials in ``materials`` tracked on this - process - """ - # Tally group-wise fission reaction rates - self._reaction_tally = Tally() - self._reaction_tally.scores = ['fission'] - - # Tally energy-weighted group-wise fission reaction rate - # Used to evaluated linear interpolation between fission yield points - self._weighted_reaction_tally = Tally() - self._weighted_reaction_tally.scores = ['fission'] - - filters = [ - MaterialFilter(materials), EnergyFilter(self._energy_bounds)] - - self._reaction_tally.filters = filters - self.local_indexes = asarray(mat_indexes) - - def set_fissionable_nuclides(self, nuclides): - """List of string of nuclides with data to be tallied - - Parameters - ---------- - nuclides : iterable of str - Nuclides with non-zero densities that are candidates - for the fission tally. Not necessary that all are nuclides - with fission yields, but at least one must be + _local_mat_index : int, optional + Current material index. Not used since all yields are + constant Returns ------- - nuclides : tuple of str - Nuclides ordered as they appear in the tally and in - the nuclide column of :attr:`results` - - Raises - ------ - ValueError - If no nuclides in ``nuclides`` are tracked on this - object + library : dict + Dictionary of ``{parent: {product: fyield}}`` """ - # Set of all nuclides with positive density - # and fission yield data - nuc_set = self._chain_set & set(nuclides) - if len(nuc_set) == 0: - raise ValueError( - "No overlap between chain nuclides with fission yields and " - "requested tally nuclides") - nuclides = tuple(sorted(nuc_set)) - self._tally_index = [self._chain_nuclides[n] for n in nuclides] - self._reaction_tally.nuclides = nuclides - return nuclides - - def unpack(self): - """Unpack fission rate tallies to produce :attr:`results` - """ - # if this process is not responsible for depleting anything - # [more processes than burnable materials] - # don't do anything - if self.local_indexes.size == 0: - return - - # get view into tally results - # new shape: [material, energy, parent nuclide] - result_view = self._reaction_tally.results[..., 1].reshape( - self.n_bmats, len(self._energy_bounds) - 1, - len(self._reaction_tally.nuclides)) - - # Get results specific to this process - fission_rates = result_view[self.local_indexes] - self.results = empty_like(fission_rates) - - # scale group fission rates proportional to total fission rate - fission_total = fission_rates.sum(axis=1) - nz_mat, nz_nuc = fission_total.nonzero() - self.results[nz_mat, :, nz_nuc] = ( - fission_rates[nz_mat, :, nz_nuc] - / fission_total[nz_mat, newaxis, nz_nuc]) - - # directly set values to zero where total fission rate is zero - z_mat, z_nuc = where(fission_total == 0.0) - self.results[z_mat, :, z_nuc] = 0.0 + return self.constant_yields def compute_yields(self, local_mat_index): """Compute single fission yields using :attr:`results` diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 15f85aa4b..660539af4 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -26,7 +26,7 @@ from .atom_number import AtomNumber from .reaction_rates import ReactionRates from .results_list import ResultsList from .helpers import ( - DirectReactionRateHelper, ChainFissionHelper, FissionYieldHelper) + DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper) def _distribute(items): @@ -85,6 +85,10 @@ 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_energy : float, optional + Energy [eV] to pull fission product yields from. Passed to the + :class:`openmc.deplete.ConstantFissionYieldHelper`. Default: + 0.0253 eV Attributes ---------- @@ -123,7 +127,7 @@ class Operator(TransportOperator): """ 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_energy=0.0253): super().__init__(chain_file, fission_q, dilute_initial, prev_results) self.round_number = False self.prev_res = None @@ -178,8 +182,8 @@ class Operator(TransportOperator): self._rate_helper = DirectReactionRateHelper( self.reaction_rates.n_nuc, self.reaction_rates.n_react) self._energy_helper = ChainFissionHelper() - self._fsn_yield_helper = FissionYieldHelper( - self.chain.nuclides, len(self.burnable_mats)) + self._fsn_yield_helper = ConstantFissionYieldHelper( + self.chain.nuclides, energy=fission_yield_energy) def __call__(self, vec, power): """Runs a simulation. 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 000000000..b9eaf9d2f --- /dev/null +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -0,0 +1,42 @@ +"""Test the FissionYieldHelpers""" + +from collections import namedtuple + +import pytest + +from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution +from openmc.deplete.helpers import ConstantFissionYieldHelper + + +@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.name = "U235" + u235.yield_energies = (0.0253, 5.0e5, 1.40e7) + u235.yield_data = FissionYieldDistribution.from_dict(u5yield_dict) + + u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}} + u238 = Nuclide() + u238.name = "U238" + u238.yield_energies = (5.00e5, ) + u238.yield_data = FissionYieldDistribution.from_dict(u8yield_dict) + + xe135 = Nuclide() + xe135.name = "Xe135" + + NuclideBundle = namedtuple("NuclideBundle", "u235 u238 xe135") + return NuclideBundle(u235, u238, xe135) +@pytest.mark.parametrize("input_energy, u5_yield_energy", ( + (0.0253, 0.0253), (0.01, 0.0253), (4e5, 5e5))) +def test_constant_helper(nuclide_bundle, input_energy, u5_yield_energy): + helper = ConstantFissionYieldHelper( + nuclide_bundle, energy=input_energy) + assert helper.energy == input_energy + assert helper.constant_yields == { + "U235": nuclide_bundle.u235.yield_data[u5_yield_energy], + "U238": nuclide_bundle.u238.yield_data[5.00e5]} # only epithermal + assert helper.constant_yields == helper.weighted_yields(1) diff --git a/tests/unit_tests/test_deplete_helpers.py b/tests/unit_tests/test_deplete_helpers.py deleted file mode 100644 index 15ae1550e..000000000 --- a/tests/unit_tests/test_deplete_helpers.py +++ /dev/null @@ -1,96 +0,0 @@ -"""Test the Operator helpers""" - -from unittest.mock import Mock - -import pytest -import numpy -from numpy.testing import assert_array_equal - -from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution -from openmc.deplete.helpers import FissionYieldHelper - - -class FissionYieldHelperProxy(FissionYieldHelper): - - def __init__(self, chain_nuclides, n_bmats): - super().__init__(chain_nuclides, n_bmats) - self._reaction_tally = Mock() - self.local_indexes = numpy.array([0]) - - def generate_tallies(self, *args, **kwargs): - # Avoid calls to the C-API - pass - - -def test_fission_yield_helper(): - """Test the collection of fission yield data using approximated tallies - """ - u5yield_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}} - u235 = Nuclide() - u235.name = "U235" - u235.yield_energies = (0.0253, 1.40e7) - u235.yield_data = FissionYieldDistribution.from_dict(u5yield_dict) - - u8yield_dict = {5.00e5: {"Xe135": 1.12e-3, "Gd155": 1.32e-12}} - u238 = Nuclide() - u238.name = "U238" - u238.yield_energies = (5.00e5, ) - u238.yield_data = FissionYieldDistribution.from_dict(u8yield_dict) - - xe135 = Nuclide() - xe135.name = "Xe135" - - n_bmats = 2 - helper = FissionYieldHelperProxy([u235, u238, xe135], n_bmats) - - assert helper.energy_bounds == (0, 0.0253, 5.00e5, 1.40e7) - - # test that tally must be created with nuclides with yields - with pytest.raises(ValueError, match="No overlap"): - helper.set_fissionable_nuclides(["Xe135", ]) - - act_nucs = helper.set_fissionable_nuclides(["U235", "U238", "Xe135"]) - assert act_nucs == ("U235", "U238") - - # Emulate getting tally data from transport run - # Test as if this Helper is responsible for one of two materials - # Tally results ordered [n_mat * n_ene, n_fiss_nuc, 3] - - tally_res = numpy.zeros((3 * n_bmats, 2, 3)) - u5_fiss_rates = numpy.array([1.0, 1.5, 2.0]) - u8_fiss_rates = numpy.array([0.0, 1.0, 1.0]) - tally_res[:3, 0, 1] = u5_fiss_rates - tally_res[:3, 1, 1] = u8_fiss_rates - - helper._reaction_tally.results = tally_res - helper.unpack() # compute yield fractions - - # Compare fraction fission rate from helper.reset - exp_results = numpy.empty((1, 3, 2)) - # Fraction of fission events in each energy range - u5_vec = u5_fiss_rates / u5_fiss_rates.sum() - u8_vec = u8_fiss_rates / u8_fiss_rates.sum() - exp_results[..., 0] = u5_vec - exp_results[..., 1] = u8_vec - assert_array_equal(helper.results, exp_results) - - # Compute and compare the library of fission yields - exp_lib = { - "U235": { - "Xe135": (u5yield_dict[0.0253]["Xe135"] * u5_vec[0] - + u5yield_dict[1.40e7]["Xe135"] * u5_vec[2]), - "Gd155": (u5yield_dict[0.0253]["Gd155"] * u5_vec[0] - + u5yield_dict[1.40e7]["Gd155"] * u5_vec[2]), - "Sm149": u5yield_dict[0.0253]["Sm149"] * u5_vec[0], - }, - "U238": { - "Xe135": u8yield_dict[5.00e5]["Xe135"] * u8_vec[1], - "Gd155": u8yield_dict[5.00e5]["Gd155"] * u8_vec[1], - } - } - - act_library = helper.compute_yields(0) - for parent, sub_yields in exp_lib.items(): - assert act_library[parent] == pytest.approx(sub_yields)