From d8b9f5db54d28c6dc16330cd056e47eed90bedec Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 13 Aug 2019 15:39:13 -0500 Subject: [PATCH] Add ConstantFissionYieldHelper concrete class Given a requested energy, will take fission yields from that energy on all nuclides and hold them as constant throughout the simulation. If the requested energy is not found for a specific nuclide, then the closest energy is used. The default is to take the thermal 0.0253 eV yields. This is now the helper attached to the Operator. The user can select what energy of yields they would like to use. --- docs/source/pythonapi/deplete.rst | 1 + openmc/deplete/helpers.py | 173 +++++------------- openmc/deplete/operator.py | 12 +- .../unit_tests/test_deplete_fission_yields.py | 42 +++++ tests/unit_tests/test_deplete_helpers.py | 96 ---------- 5 files changed, 99 insertions(+), 225 deletions(-) create mode 100644 tests/unit_tests/test_deplete_fission_yields.py delete mode 100644 tests/unit_tests/test_deplete_helpers.py 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)