diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 44f8b70e0..5699f3f2a 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -77,6 +77,7 @@ data, such as number densities and reaction rates for each material. :template: myclass.rst AtomNumber + AveragedFissionYieldHelper ChainFissionHelper ConstantFissionYieldHelper DirectReactionRateHelper diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 5efda8702..d1dc2cdc4 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -504,6 +504,8 @@ class TalliedFissionYieldHelper(FissionYieldHelper): 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 @@ -513,7 +515,8 @@ class TalliedFissionYieldHelper(FissionYieldHelper): self.n_bmats = n_bmats self._local_indexes = None self._fission_rate_tally = None - self._tally_index = {} + self._tally_nucs = [] + self.results = None def generate_tallies(self, materials, mat_indexes): """Construct the fission rate tally @@ -554,10 +557,10 @@ class TalliedFissionYieldHelper(FissionYieldHelper): if len(overlap) == 0: # tally no nuclides, but keep the Tally alive self._fission_rate_tally.nuclides = None - self._tally_index = [] + self._tally_nucs = [] return tuple() nuclides = tuple(sorted(overlap)) - self._tally_index = [self._chain_nuclides[n] for n in nuclides] + self._tally_nucs = [self._chain_nuclides[n] for n in nuclides] self._fission_rate_tally.nuclides = nuclides return nuclides @@ -587,6 +590,7 @@ class TalliedFissionYieldHelper(FissionYieldHelper): return cls(operator.chain.nuclides, len(operator.burnable_mats), **kwargs) + class Integrator(ABC): """Abstract class for solving the time-integration for depletion diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index b66c9ba5e..19e51628d 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -9,14 +9,15 @@ from numpy import dot, zeros, newaxis from openmc.checkvalue import check_type, check_greater_than from openmc.capi import ( - Tally, MaterialFilter, EnergyFilter) + Tally, MaterialFilter, EnergyFilter, EnergyFunctionFilter) from .abc import ( ReactionRateHelper, EnergyHelper, FissionYieldHelper, TalliedFissionYieldHelper) __all__ = ( "DirectReactionRateHelper", "ChainFissionHelper", - "ConstantFissionYieldHelper", "FissionYieldCutoffHelper") + "ConstantFissionYieldHelper", "FissionYieldCutoffHelper", + "AveragedFissionYieldHelper") # ------------------------------------- # Helpers for generating reaction rates @@ -386,7 +387,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): def unpack(self): """Obtain fast and thermal fission fractions from tally""" fission_rates = self._fission_rate_tally.results[..., 1].reshape( - self.n_bmats, 2, len(self._tally_index)) + 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() @@ -422,10 +423,10 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): rates = self.results[local_mat_index] yields = self.constant_yields # iterate over thermal then fast yields, prefer __mul__ to __rmul__ - for therm_frac, nuc in zip(rates[0], self._tally_index): + for therm_frac, nuc in zip(rates[0], self._tally_nucs): yields[nuc.name] = self._thermal_yields[nuc.name] * therm_frac - for fast_frac, nuc in zip(rates[1], self._tally_index): + for fast_frac, nuc in zip(rates[1], self._tally_nucs): yields[nuc.name] += self._fast_yields[nuc.name] * fast_frac return yields @@ -442,3 +443,149 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): for key, sub in self._fast_yields.items(): out[key] = sub.copy() return out + + +class AveragedFissionYieldHelper(TalliedFissionYieldHelper): + r"""Class that computes fission yields based on average fission energy + + Computes average energy at which fission events occured + reactions for all nuclides with multiple sets of fission yields + by + + .. 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, a geometric mean of the yield distributions is + used. + + 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 + ---------- + n_bmats : int + Number of burnable materials tracked in the problem. + 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. + """ + + def __init__(self, chain_nuclides, n_bmats): + super().__init__(chain_nuclides, n_bmats) + 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 + Indexes for materials in ``materials`` tracked on this + process + """ + super().generate_tallies(materials, mat_indexes) + fission_tally = self._fission_rate_tally + + weighted_tally = Tally() + weighted_tally.filters = fission_tally.filters.copy() + weighted_tally.nuclides = fission_tally.nuclides + weighted_tally.scores = ['fission'] + + ene_bin = EnergyFilter() + ene_bin.bins = (0, self._upper_energy) + fission_tally.filters.append(ene_bin) + + ene_filter = EnergyFunctionFilter() + ene_filter.set_data((0, self._upper_energy), (0, self._upper_energy)) + weighted_tally.filters.append(ene_filter) + self._weighted_tally = weighted_tally + + def unpack(self): + """Unpack tallies and populate :attr:`results` with average energies""" + 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}}`` + """ + 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): + """Return a new helper with data from an operator + + Parameters + ---------- + operator : openmc.deplete.Operator + Operator with a depletion chain + + Returns + ------- + AveragedFissionYieldHelper + """ + return cls(operator.chain.nuclides, len(operator.burnable_mats)) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index e967dc85d..b7df9b8fc 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -27,7 +27,7 @@ from .reaction_rates import ReactionRates from .results_list import ResultsList from .helpers import ( DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, - FissionYieldCutoffHelper) + FissionYieldCutoffHelper, AveragedFissionYieldHelper) def _distribute(items): @@ -92,6 +92,7 @@ class Operator(TransportOperator): * "constant": :class:`openmc.deplete.ConstantFissionYieldHelper` * "cutoff": :class:`openmc.deplete.FissionYieldCutoffHelper` + * "average": :class:`openmc.deplete.AveragedFissionYieldHelper` The documentation on these classes describe their methodology and differences. ``"constant"`` will treat fission yields as @@ -140,6 +141,7 @@ class Operator(TransportOperator): Whether to differentiate burnable materials with multiple instances """ _fission_helpers_ = { + "average": AveragedFissionYieldHelper, "constant": ConstantFissionYieldHelper, "cutoff": FissionYieldCutoffHelper, } diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 233074236..0a87c4fb6 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -8,7 +8,11 @@ import numpy from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution from openmc.deplete.helpers import ( - FissionYieldCutoffHelper, ConstantFissionYieldHelper) + FissionYieldCutoffHelper, ConstantFissionYieldHelper, + AveragedFissionYieldHelper) + + +MATERIALS = ["1", "2"] @pytest.fixture(scope="module") @@ -45,11 +49,6 @@ def test_constant_helper(nuclide_bundle, input_energy, u5_yield_energy): assert helper.constant_yields == helper.weighted_yields(1) -# --------------------------------- -# Test the FissionYieldCutoffHelper -# --------------------------------- - - def test_cutoff_construction(nuclide_bundle): # defaults helper = FissionYieldCutoffHelper(nuclide_bundle, 1) @@ -92,22 +91,23 @@ def test_cutoff_failure(key): FissionYieldCutoffHelper(None, None, **{key: -1}) -class CutoffProxy(FissionYieldCutoffHelper): - """Proxy that supplies a set of tallies""" - +class ProxyMixin(object): + """Mixing that overloads the tally generation""" def generate_tallies(self, materials, mat_indexes): self._fission_rate_tally = Mock() self._local_indexes = numpy.asarray(mat_indexes) +class CutoffProxy(ProxyMixin, FissionYieldCutoffHelper): + """Proxy that supplies a set of tallies""" + + # emulate some split between fast and thermal U235 fissions @pytest.mark.parametrize("therm_frac", (0.5, 0.2, 0.8)) def test_cutoff_helper(nuclide_bundle, therm_frac): - materials = ["1", "2"] # TODO Use real C API materials? - n_bmats = len(materials) - local_mats = [0] + n_bmats = len(MATERIALS) proxy = CutoffProxy(nuclide_bundle, n_bmats) - proxy.generate_tallies(materials, local_mats) + proxy.generate_tallies(MATERIALS, [0]) non_zero_nucs = [n.name for n in nuclide_bundle] tally_nucs = proxy.update_nuclides_from_operator(non_zero_nucs) assert tally_nucs == ("U235",) @@ -130,5 +130,44 @@ def test_cutoff_helper(nuclide_bundle, therm_frac): assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5] assert actual_yields["U235"] == ( proxy.thermal_yields["U235"] * therm_frac - + proxy.fast_yields["U235"] * (1 - therm_frac) - ) + + proxy.fast_yields["U235"] * (1 - therm_frac)) + + +class AverageProxy(ProxyMixin, AveragedFissionYieldHelper): + """Proxy for generating mock set of tallies""" + def generate_tallies(self, materials, mat_indexes): + super().generate_tallies(materials, mat_indexes) + self._weighted_tally = Mock() + + +@pytest.mark.parametrize("avg_energy", (0.01, 100, 15e6)) +def test_averaged_helper(nuclide_bundle, avg_energy): + proxy = AverageProxy(nuclide_bundle, len(MATERIALS)) + proxy.generate_tallies(MATERIALS, [0]) + tallied_nucs = proxy.update_nuclides_from_operator( + [n.name for n in nuclide_bundle]) + assert tallied_nucs == ("U235", ) + # enforce some average energy + proxy_flux = 1e16 + fission_results = numpy.ones((len(MATERIALS), 1, 3)) * proxy_flux + weighted_results = fission_results * avg_energy + proxy._fission_rate_tally.results = fission_results + proxy._weighted_tally.results = weighted_results + proxy.unpack() + expected_results = numpy.ones((1, 1)) * avg_energy + assert proxy.results == pytest.approx(expected_results) + + actual_yields = proxy.weighted_yields(0) + # constant U238 => no interpolation + assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5] + # construct expected yields + if avg_energy < 0.0253: # take thermal U235 yields + exp_u235_yields = nuclide_bundle.u235.yield_data[0.0253] + elif avg_energy > 14e6: # take fastest U235 yields + exp_u235_yields = nuclide_bundle.u235.yield_data[14e6] + else: # reconstruct between thermal and epithermal + thermal = nuclide_bundle.u235.yield_data[0.0253] + epithermal = nuclide_bundle.u235.yield_data[5e5] + split = (avg_energy - 0.0253) / (5e5 - 0.0253) + exp_u235_yields = thermal * (1 - split) + epithermal * split + assert actual_yields["U235"] == exp_u235_yields