Implement hybrid depletion tallies

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
 # ------------------------------------------

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.")

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