From 599d473888be6e443cac4465987078a58eec9c19 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 2 Aug 2019 17:51:33 -0500 Subject: [PATCH] Add FissionYieldHelper for energy-dependent fission yields Operator now has a new helper for tallying fission rates and computing fission yields using energy dependency. The FissionYieldHelper creates fission rate tallies across burnable materials with an energy filter. The energy filter corresponds to each energy group for which there are fission yields. During the operator's unpacking of tallies, this helper computes a relative fission contribution by normalizing fission rates, such that the sum of all group fission rates in each material and each nuclide is unity. A single fission yield for each parent nuclide is computed by weighting the energy-dependent fission yields by these relative contributions. A nested dictionary {parent: {product: yield}} is added to a library, and then attached to the depletion chain at the end of the unpacking. These libraries are used during the depletion process instead of the only-thermal yields from before. A test was added that creates a proxy helper that doesn't use the C-API, but unpacks the "fission tally" and computes the fission yield library in an identical manner as the actual class. Since the proxy helper is a subclass of the one used in real transport, the functionality is consistent. --- openmc/deplete/chain.py | 3 - openmc/deplete/helpers.py | 193 ++++++++++++++++++++++- openmc/deplete/operator.py | 22 ++- tests/unit_tests/test_deplete_helpers.py | 99 ++++++++++++ 4 files changed, 309 insertions(+), 8 deletions(-) create mode 100644 tests/unit_tests/test_deplete_helpers.py diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 62158255a7..31b074bddb 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -452,9 +452,6 @@ 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.. for product, y in fission_yields[nuc.name].items(): yield_val = y * path_rate if yield_val != 0.0: diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 66722c00a2..07c8e20242 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -1,11 +1,13 @@ """ Class for normalizing fission energy deposition """ +from collections import defaultdict from itertools import product -from numpy import dot, zeros +from numpy import dot, zeros, newaxis, divide, asarray -from openmc.capi import Tally, MaterialFilter +from openmc.capi import ( + Tally, MaterialFilter, EnergyFilter) from .abc import ReactionRateHelper, EnergyHelper # ------------------------------------- @@ -142,3 +144,190 @@ 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 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. + + 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 + + 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)`` + libraries : list of dict + List of fission yield dictionaries of the form + ``{parent: {product: yield}}``. Populated in + :meth:`compute_yields` and reset during + :meth:`unpack`. + """ + + def __init__(self, chain_nuclides, n_bmats): + self.libraries = [] + 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: + 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 + + @property + def energy_bounds(self): + return self._energy_bounds + + def generate_tallies(self, materials, mat_indexes): + """Construct the fission rate tally + + 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 + + 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 + """ + # 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` + + Resets :attr:`libraries` under the assumption this is called + during the :class:`openmc.deplete.Operator` unpackign process + """ + # clear old libraries + self.libraries = [] + + # 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 + self.results = result_view[self.local_indexes, ...] + + # scale fission yields proportional to total fission rate + fsn_rate = self.results.sum(axis=1) + # TODO Guard against divide by zero + self.results /= fsn_rate[:, newaxis, :] + + def compute_yields(self, local_mat_index): + """Compute single fission yields using :attr:`results` + + Produces a new library in :attr:`self.libraries` + + 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` + """ + tally_results = self.results[local_mat_index] + + # Dictionary {parent_nuclide : [product, yield_vector]} + initial_library = {} + for i_energy, energy in enumerate(self._energy_bounds[1:]): + for i_nuc, fiss_frac in enumerate(tally_results[i_energy]): + parent = self._tally_index[i_nuc] + yield_data = parent.yield_data.get(energy) + if yield_data is None: + continue + if parent not in initial_library: + initial_library[parent] = yield_data * fiss_frac + continue + initial_library[parent] += yield_data * fiss_frac + + # convert to dictionary that can be passed to Chain.form_matrix + # {parent: {product: yield}} + library = {} + for k, yield_obj in initial_library.items(): + library[k.name] = dict(zip(yield_obj.products, yield_obj.yields)) + + self.libraries.append(library) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index f8288841f7..0c8c72987b 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -25,7 +25,8 @@ 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, FissionYieldHelper) def _distribute(items): @@ -177,6 +178,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)) def __call__(self, vec, power): """Runs a simulation. @@ -204,8 +207,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._fsn_yield_helper.set_fissionable_nuclides(nuclides) # Run OpenMC openmc.capi.reset() @@ -409,6 +414,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._fsn_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 +563,7 @@ class Operator(TransportOperator): # Keep track of energy produced from all reactions in eV per source # particle self._energy_helper.reset() + self._fsn_yield_helper.unpack() # Create arrays to store fission Q values, reaction rates, and nuclide # numbers, zeroed out in material iteration @@ -576,6 +586,9 @@ class Operator(TransportOperator): tally_rates = self._rate_helper.get_material_rates( mat_index, nuc_ind, react_ind) + # Compute fission yields for this material + self._fsn_yield_helper.compute_yields(i) + # Accumulate energy from fission self._energy_helper.update(tally_rates[:, fission_ind], mat_index) @@ -589,6 +602,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 = self._fsn_yield_helper.libraries + return OperatorResult(k_combined, rates) def _get_nuclides_with_data(self): diff --git a/tests/unit_tests/test_deplete_helpers.py b/tests/unit_tests/test_deplete_helpers.py new file mode 100644 index 0000000000..5db979098b --- /dev/null +++ b/tests/unit_tests/test_deplete_helpers.py @@ -0,0 +1,99 @@ +"""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], + } + } + + assert len(helper.libraries) == 0 + helper.compute_yields(0) + assert len(helper.libraries) == 1 + act_library = helper.libraries[0] + for parent, sub_yields in exp_lib.items(): + assert act_library[parent] == pytest.approx(sub_yields)