From 454e6cc6b6739cd6f358951757354f654171ce80 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 14 Aug 2019 17:39:57 -0500 Subject: [PATCH] Add AveragedFissionYieldHelper operator helper Computes the effective fission yields for nuclides with multiple sets of yields by 1) Tallying fission rate and energy-weighted fission rate 2) Determining average energy at which fission events occur 3) Interpolating [lin-lin] between adjacent fission yields based on this average energy Uses a second tally with an EnergyFunctionFilter to tally E * sigma_f * phi. Added unit tests with the proxy-style class that doesn't use the C API to generate tallies but tests all other aspects of the implementation. Tests cover three possible average energies (below minimum supplied yield energy, above max supplied yield energy, and between two) The Operator uses fission_yield_mode="average" to select this helper. The helper also has a from_operator class method. --- docs/source/pythonapi/deplete.rst | 1 + openmc/deplete/abc.py | 10 +- openmc/deplete/helpers.py | 157 +++++++++++++++++- openmc/deplete/operator.py | 4 +- .../unit_tests/test_deplete_fission_yields.py | 69 ++++++-- 5 files changed, 217 insertions(+), 24 deletions(-) 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