diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 7951f18adf..c2ddb0e9ec 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -201,6 +201,145 @@ class FluxCollapseHelper(ReactionRateHelper): return self._results_cache + +class HybridReactionHelper(ReactionRateHelper): + """Class that generates tallies for one-group rates + + .. versionadded:: 0.12.1 + + Parameters + ---------- + n_nucs : int + Number of burnable nuclides tracked by :class:`openmc.deplete.Operator` + n_react : int + Number of reactions tracked by :class:`openmc.deplete.Operator` + energies : iterable of float + Energy group boundaries for flux spectrum in [eV] + reactions : iterable of str + Reactions for which rates should be directly tallied + nuclides : iterable of str + Nuclides for which some reaction rates should be directly tallied. If + None, then all ``reactions`` will be used for all nuclides. + + Attributes + ---------- + nuclides : list of str + All nuclides with desired reaction rates. + + """ + def __init__(self, n_nucs, n_reacts, energies, reactions, nuclides=None): + super().__init__(n_nucs, n_reacts) + self._energies = asarray(energies) + self._reactions_direct = list(reactions) + self._nuclides_direct = list(nuclides) if nuclides is not None else None + + @ReactionRateHelper.nuclides.setter + def nuclides(self, nuclides): + ReactionRateHelper.nuclides.fset(self, nuclides) + if self._nuclides_direct is None: + self._rate_tally.nuclides = nuclides + + def generate_tallies(self, materials, scores): + """Produce multigroup flux spectrum tally + + Uses the :mod:`openmc.lib` module to generate a multigroup flux tally + for each burnable material. + + Parameters + ---------- + materials : iterable of :class:`openmc.Material` + Burnable materials in the problem. Used to construct a + :class:`openmc.MaterialFilter` + scores : iterable of str + Reaction identifiers, e.g. ``"(n, fission)"``, ``"(n, gamma)"``, + needed for the reaction rate tally. + """ + self._materials = materials + + # Convert reactions to MT values (needed when collapsing) + mt_values = {v: k for k, v in REACTION_NAME.items()} + mt_values['fission'] = 18 + self._mts = [mt_values[x] for x in scores] + self._scores = scores + + # Create flux tally with material and energy filters + self._flux_tally = Tally() + self._flux_tally.writable = False + self._flux_tally.filters = [ + MaterialFilter(materials), + EnergyFilter(self._energies) + ] + self._flux_tally.scores = ['flux'] + + # Create reaction rate tally + self._rate_tally = Tally() + self._rate_tally.writable = False + self._rate_tally.scores = self._reactions_direct + self._rate_tally.filters = [MaterialFilter(materials)] + if self._nuclides_direct is not None: + self._rate_tally.nuclides = self._nuclides_direct + + def get_material_rates(self, mat_index, nuc_index, react_index): + """Return an array of reaction rates for a material + + Parameters + ---------- + mat_index : int + Index for material + nuc_index : iterable of int + Index for each nuclide in :attr:`nuclides` in the + desired reaction rate matrix + react_index : iterable of int + Index for each reaction scored in the tally + + Returns + ------- + rates : numpy.ndarray + Array with shape ``(n_nuclides, n_rxns)`` with the reaction rates in + this material + + """ + self._results_cache.fill(0.0) + + # Get flux for specified material + shape = (len(self._materials), len(self._energies) - 1) + mean_value = self._flux_tally.mean.reshape(shape) + flux = mean_value[mat_index] + + # Get direct reaction rates + nuclides_direct = self._rate_tally.nuclides + shape = (len(nuclides_direct), len(self._reactions_direct)) + rx_rates = self._rate_tally.mean[mat_index].reshape(shape) + + mat = self._materials[mat_index] + + # Build nucname: density mapping to enable O(1) lookup in loop below + densities = dict(zip(mat.nuclides, mat.densities)) + + for name, i_nuc in zip(self.nuclides, nuc_index): + # Determine density of nuclide + density = densities[name] + + for mt, score, i_rx in zip(self._mts, self._scores, react_index): + if score in self._reactions_direct and name in nuclides_direct: + # Determine index in rx_rates + i_rx_direct = self._reactions_direct.index(score) + i_nuc_direct = nuclides_direct.index(name) + + # Get reaction rate from tally + self._results_cache[i_nuc, i_rx] = rx_rates[i_nuc_direct, i_rx_direct] + else: + # Use flux to collapse reaction rate (per N) + nuc = openmc.lib.nuclides[name] + rate_per_nuc = nuc.collapse_rate( + mt, mat.temperature, self._energies, flux) + + # Multiply by density to get absolute reaction rate + self._results_cache[i_nuc, i_rx] = rate_per_nuc * density + + return self._results_cache + + # ------------------------------------------ # Helpers for obtaining normalization factor # ------------------------------------------ diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index f1a2912e9a..bc106278b1 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -27,7 +27,7 @@ from .results_list import ResultsList from .helpers import ( DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper, - SourceRateHelper, FluxCollapseHelper) + SourceRateHelper, FluxCollapseHelper, HybridReactionHelper) __all__ = ["Operator", "OperatorResult"] @@ -120,7 +120,7 @@ class Operator(TransportOperator): rates after a transport solve. .. versionadded:: 0.12.1 - reaction_rate_energies : iterable of float + reaction_rate_opts : iterable of float Energy group boundaries that are to be used for calculating a multigroup flux spectrum when the "flux" based ``reaction_rate_mode`` is being used. @@ -182,7 +182,7 @@ class Operator(TransportOperator): diff_burnable_mats=False, normalization_mode="fission-q", fission_q=None, dilute_initial=1.0e3, fission_yield_mode="constant", fission_yield_opts=None, - reaction_rate_mode="direct", reaction_rate_energies=None, + reaction_rate_mode="direct", reaction_rate_opts=None, reduce_chain=False, reduce_chain_level=None): check_value('fission yield mode', fission_yield_mode, self._fission_helpers.keys()) @@ -255,18 +255,27 @@ class Operator(TransportOperator): if reaction_rate_mode == "direct": self._rate_helper = DirectReactionRateHelper( self.reaction_rates.n_nuc, self.reaction_rates.n_react) - elif reaction_rate_mode == "flux": + elif reaction_rate_mode in ("flux", "hybrid"): + if reaction_rate_opts is None: + reaction_rate_opts = {} + # Ensure energy group boundaries were specified - if reaction_rate_energies is None: + if 'energies' not in reaction_rate_opts: raise ValueError( "Energy group boundaries must be specified in the " - "reaction_rate_energies argument when reaction_rate_mode is" - "set to 'flux'.") + "reaction_rate_opts argument when reaction_rate_mode is" + "set to 'flux' or 'hybrid'.") - self._rate_helper = FluxCollapseHelper( + if reaction_rate_mode == "flux": + cls = FluxCollapseHelper + else: + cls = HybridReactionHelper + + + self._rate_helper = cls( self.reaction_rates.n_nuc, self.reaction_rates.n_react, - reaction_rate_energies + **reaction_rate_opts ) else: raise ValueError("Invalid reaction rate mode.") diff --git a/tests/unit_tests/test_deplete_activation.py b/tests/unit_tests/test_deplete_activation.py index 18111fea74..0c4d0cc89f 100644 --- a/tests/unit_tests/test_deplete_activation.py +++ b/tests/unit_tests/test_deplete_activation.py @@ -61,7 +61,7 @@ def test_activation(run_in_tmpdir, model, reaction_rate_mode): model.geometry, model.settings, 'test_chain.xml', normalization_mode="source-rate", reaction_rate_mode=reaction_rate_mode, - reaction_rate_energies=energies, + reaction_rate_opts={'energies': energies}, ) # To determine the source rate necessary to reduce W186 density in half, we