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)