Add k-infinity estimator for IndependentOperator

docs/source/usersguide/depletion.rst

@@ -270,6 +270,22 @@ transport-depletion calculation and follow the same steps from there.
    the depletion chain with at least one reaction, that reaction will not be
    simulated.
 
+If the microscopic cross section data includes 'fission' and 'nu-fission'
+cross sections, :class:`~openmc.deplete.IndependentOperator` can also estimate
+the infinite multiplication factor at each depletion step by passing
+``calculate_kinf=True``::
+
+    op = openmc.deplete.IndependentOperator(materials, fluxes, micros,
+                                            calculate_kinf=True)
+
+The estimate is computed as the ratio of the neutron production rate to the
+neutron loss rate based on the one-group reaction rates and is reported as the
+eigenvalue in the depletion results, which can be retrieved with
+:meth:`~openmc.deplete.Results.get_keff`. Consistent with the definition of
+the multiplication factor used elsewhere in OpenMC, neutrons produced in
+(n,xn) reactions are not counted as production; instead, each (n,xn) reaction
+reduces the loss term by :math:`x - 1`.
+
 .. _micros:
 
 Loading and Generating Microscopic Cross Sections

openmc/deplete/independent_operator.py

@@ -8,6 +8,7 @@
 from __future__ import annotations
 from collections.abc import Iterable
 import copy
+import re
 
 import numpy as np
 from uncertainties import ufloat
@@ -22,6 +23,18 @@
 from .results import Results
 from .helpers import ChainFissionHelper, ConstantFissionYieldHelper, SourceRateHelper
 
+# Regular expression matching reactions that emit one or more neutrons, e.g.,
+# (n,2n) or (n,np), with the number of emitted neutrons captured
+_XN_REACTION = re.compile(r'\(n,(\d*)n')
+
+
+def _neutrons_emitted(reaction: str) -> int:
+    """Number of neutrons in the exit channel of a transmutation reaction."""
+    match = _XN_REACTION.match(reaction)
+    if match is None:
+        return 0
+    return int(match.group(1)) if match.group(1) else 1
+
 
 class IndependentOperator(OpenMCOperator):
     """Transport-independent transport operator based on multigroup data.
@@ -55,6 +68,14 @@ class IndependentOperator(OpenMCOperator):
         Defaults to ``openmc.config['chain_file']``.
     keff : 2-tuple of float, optional
        keff eigenvalue and uncertainty from transport calculation.
+    calculate_kinf : bool, optional
+        If True, the infinite multiplication factor is estimated from the
+        material compositions and one-group cross sections at each depletion
+        step and reported as the eigenvalue in the depletion results. Requires
+        that each :class:`~openmc.deplete.MicroXS` instance contains 'fission'
+        and 'nu-fission' cross sections. Mutually exclusive with ``keff``.
+
+        .. versionadded:: 0.15.4
     prev_results : Results, optional
         Results from a previous depletion calculation.
     normalization_mode : {"fission-q", "source-rate"}
@@ -116,7 +137,8 @@ def __init__(self,
                  fission_q=None,
                  prev_results=None,
                  reduce_chain_level=None,
-                 fission_yield_opts=None):
+                 fission_yield_opts=None,
+                 calculate_kinf=False):
         # Validate micro-xs parameters
         check_type('materials', materials, Iterable, openmc.Material)
         check_type('micros', micros, Iterable, MicroXS)
@@ -134,6 +156,20 @@ def __init__(self,
 
         self._keff = keff
 
+        check_type('calculate_kinf', calculate_kinf, bool)
+        if calculate_kinf:
+            if keff is not None:
+                raise ValueError("The 'keff' and 'calculate_kinf' arguments "
+                                 "are mutually exclusive.")
+            for micro in micros:
+                missing = {'fission', 'nu-fission'} - set(micro.reactions)
+                if missing:
+                    raise ValueError(
+                        "Estimating k-infinity requires 'fission' and "
+                        "'nu-fission' cross sections in each MicroXS "
+                        f"instance (missing {sorted(missing)}).")
+        self._calculate_kinf = calculate_kinf
+
         if fission_yield_opts is None:
             fission_yield_opts = {}
         helper_kwargs = {'normalization_mode': normalization_mode,
@@ -165,7 +201,8 @@ def from_nuclides(cls, volume, nuclides,
                       fission_q=None,
                       prev_results=None,
                       reduce_chain_level=None,
-                      fission_yield_opts=None):
+                      fission_yield_opts=None,
+                      calculate_kinf=False):
         """
         Alternate constructor from a dictionary of nuclide concentrations
 
@@ -187,6 +224,13 @@ def from_nuclides(cls, volume, nuclides,
         keff : 2-tuple of float, optional
            keff eigenvalue and uncertainty from transport calculation.
            Default is None.
+        calculate_kinf : bool, optional
+            If True, the infinite multiplication factor is estimated from the
+            material compositions and one-group cross sections at each
+            depletion step. Requires that ``micro_xs`` contains 'fission' and
+            'nu-fission' cross sections. Mutually exclusive with ``keff``.
+
+            .. versionadded:: 0.15.4
         normalization_mode : {"fission-q", "source-rate"}
             Indicate how reaction rates should be calculated.
             ``"fission-q"`` uses the fission Q values from the depletion
@@ -222,7 +266,8 @@ def from_nuclides(cls, volume, nuclides,
                    fission_q=fission_q,
                    prev_results=prev_results,
                    reduce_chain_level=reduce_chain_level,
-                   fission_yield_opts=fission_yield_opts)
+                   fission_yield_opts=fission_yield_opts,
+                   calculate_kinf=calculate_kinf)
 
     @staticmethod
     def _consolidate_nuclides_to_material(nuclides, nuc_units, volume):
@@ -407,14 +452,73 @@ def __call__(self, vec, source_rate) -> OperatorResult:
         if source_rate == 0.0:
             rates = self.reaction_rates.copy()
             rates.fill(0.0)
-            return OperatorResult(ufloat(0.0, 0.0), rates)
+            if self._calculate_kinf:
+                keff = self._estimate_k_inf()
+            else:
+                keff = ufloat(0.0, 0.0)
+            return OperatorResult(keff, rates)
 
         rates = self._calculate_reaction_rates(source_rate)
-        keff = self._keff
+        if self._calculate_kinf:
+            keff = self._estimate_k_inf()
+        else:
+            keff = self._keff
 
         op_result = OperatorResult(keff, rates)
         return copy.deepcopy(op_result)
 
+    def _estimate_k_inf(self):
+        r"""Estimate the infinite multiplication factor.
+
+        The estimate is computed as the ratio of the neutron production rate
+        to the neutron loss rate:
+
+        .. math::
+            k_\infty = \frac{\sum_i N_i (\nu\sigma_f)_i}
+                            {\sum_i N_i \sum_j (1 - x_j) \sigma_{i,j}}
+
+        where :math:`N_i` is the number of atoms of nuclide :math:`i`,
+        :math:`(\nu\sigma_f)_i` is its one-group fission neutron production
+        cross section, :math:`\sigma_{i,j}` is the one-group cross section of
+        transmutation reaction :math:`j`, and :math:`x_j` is the number of
+        neutrons emitted by reaction :math:`j`. This is consistent with the
+        definition of the multiplication factor used elsewhere in OpenMC:
+        neutrons produced in (n,xn) reactions are not counted as production;
+        instead, each (n,xn) reaction reduces the loss term by :math:`x - 1`.
+
+        Returns
+        -------
+        uncertainties.UFloat
+            Estimated k-infinity with zero uncertainty
+
+        """
+        production = 0.0
+        loss = 0.0
+        for mat in self.local_mats:
+            i_mat = self._mat_index_map[mat]
+            flux = self.fluxes[i_mat]
+            micro_xs = self.cross_sections[i_mat]
+            for nuc in micro_xs.nuclides:
+                if nuc not in self.number.index_nuc:
+                    continue
+                atoms = self.number[mat, nuc]
+                if atoms <= 0.0:
+                    continue
+                for rxn in micro_xs.reactions:
+                    rate = atoms * (micro_xs[nuc, rxn] * flux).sum()
+                    if rxn == 'nu-fission':
+                        production += rate
+                    elif rxn != 'damage-energy':
+                        loss += (1 - _neutrons_emitted(rxn)) * rate
+
+        # Sum contributions over all MPI processes
+        production = comm.allreduce(production)
+        loss = comm.allreduce(loss)
+
+        if loss <= 0.0:
+            return ufloat(0.0, 0.0)
+        return ufloat(production / loss, 0.0)
+
     def _update_materials(self):
         """Updates material compositions in OpenMC on all processes."""
 

openmc/deplete/microxs.py

@@ -17,7 +17,7 @@
 from openmc.checkvalue import check_type, check_value, check_iterable_type, PathLike
 from openmc import StatePoint
 from openmc.mgxs import GROUP_STRUCTURES
-from openmc.data import REACTION_MT
+from openmc.data import DataLibrary, REACTION_MT, Reaction
 import openmc
 from .chain import Chain, REACTIONS, _get_chain
 from .coupled_operator import _find_cross_sections, _get_nuclides_with_data
@@ -28,6 +28,7 @@
 _valid_rxns = list(REACTIONS)
 _valid_rxns.append('fission')
 _valid_rxns.append('damage-energy')
+_valid_rxns.append('nu-fission')
 
 
 # TODO: Replace with type statement when support is Python 3.12+
@@ -81,7 +82,10 @@ def get_microxs_and_flux(
         nuclides from the depletion chain file are used.
     reactions : list of str
         Reactions to get cross sections for. If not specified, all neutron
-        reactions listed in the depletion chain file are used.
+        reactions listed in the depletion chain file are used. In addition to
+        transmutation reactions, 'nu-fission' may be specified to obtain the
+        fission neutron production cross section, which is needed to estimate
+        k-infinity with :class:`~openmc.deplete.IndependentOperator`.
     energies : iterable of float or str
         Energy group boundaries in [eV] or the name of the group structure.
         If left as None energies will default to [0.0, 100e6]
@@ -303,6 +307,82 @@ def get_microxs_and_flux(
     return fluxes, micros
 
 
+def _collapse_nu_fission(
+    path: PathLike,
+    nuclide: str,
+    temperature: float,
+    energies: Sequence[float],
+    flux: np.ndarray
+) -> float:
+    r"""Compute a one-group fission neutron production cross section.
+
+    The fission neutron production cross section,
+    :math:`\nu(E)\sigma_f(E)`, is integrated against a flux that is assumed
+    to be constant in energy within each group, matching the treatment used
+    for other reactions in :meth:`openmc.lib.Nuclide.collapse_rate`.
+
+    Parameters
+    ----------
+    path : PathLike
+        Path to the HDF5 data file containing the nuclide.
+    nuclide : str
+        Name of the nuclide, e.g., 'U235'.
+    temperature : float
+        Temperature in [K]. The closest available temperature is used for the
+        fission cross section.
+    energies : iterable of float
+        Energy group boundaries in [eV] in ascending order.
+    flux : numpy.ndarray
+        Flux in each energy group, normalized to sum to unity.
+
+    Returns
+    -------
+    float
+        Flux-averaged fission neutron production cross section in [b]. Zero if
+        the nuclide has no fission data.
+
+    """
+    with h5py.File(path, 'r') as h5:
+        group = h5[nuclide]
+        if 'reactions/reaction_018' not in group:
+            return 0.0
+
+        # Select the available temperature closest to the requested one
+        temp_keys = list(group['energy'])
+        temps = np.array([float(t[:-1]) for t in temp_keys])
+        temp_key = temp_keys[np.argmin(np.abs(temps - temperature))]
+
+        energy_grid = {temp_key: group['energy'][temp_key][()]}
+        rx = Reaction.from_hdf5(group['reactions/reaction_018'], energy_grid)
+
+    xs = rx.xs[temp_key]
+
+    # Total nu(E) is the sum of the yields of all neutron products. If a
+    # product with emission mode 'total' is present, use it alone to avoid
+    # double counting prompt and delayed neutrons.
+    neutron_products = [p for p in rx.products if p.particle == 'neutron']
+    total_products = [p for p in neutron_products if p.emission_mode == 'total']
+    if total_products:
+        neutron_products = total_products
+    if not neutron_products:
+        return 0.0
+
+    def nu(e):
+        return sum(p.yield_(e) for p in neutron_products)
+
+    # Integrate nu(E)*sigma_f(E) against a histogram flux
+    nu_fission = 0.0
+    for g, flux_g in enumerate(flux):
+        if flux_g == 0.0:
+            continue
+        e_low, e_high = energies[g], energies[g + 1]
+        inside = xs.x[(xs.x > e_low) & (xs.x < e_high)]
+        e = np.concatenate([[e_low], inside, [e_high]])
+        nu_fission += np.trapezoid(nu(e) * xs(e), e) * flux_g / (e_high - e_low)
+
+    return nu_fission
+
+
 class MicroXS:
     """Microscopic cross section data for use in transport-independent depletion.
 
@@ -385,7 +465,14 @@ def from_multigroup_flux(
             nuclides from the depletion chain file are used.
         reactions : list of str, optional
             Reactions to get cross sections for. If not specified, all neutron
-            reactions listed in the depletion chain file are used.
+            reactions listed in the depletion chain file are used. In addition
+            to transmutation reactions, 'nu-fission' may be specified to
+            obtain the fission neutron production cross section, which is
+            needed to estimate k-infinity with
+            :class:`~openmc.deplete.IndependentOperator`.
+
+            .. versionchanged:: 0.15.4
+                Added support for 'nu-fission'.
         **init_kwargs : dict
             Keyword arguments passed to :func:`openmc.lib.init`
 
@@ -418,10 +505,12 @@ def from_multigroup_flux(
             nuclides = [nuc.name for nuc in nuclides]
 
         # Get reaction MT values. If no reactions specified, default to the
-        # reactions available in the chain file
+        # reactions available in the chain file. The 'nu-fission' reaction is
+        # handled separately since it does not correspond to a single MT value.
         if reactions is None:
             reactions = chain.reactions
-        mts = [REACTION_MT[name] for name in reactions]
+        mts = [REACTION_MT[name] if name != 'nu-fission' else None
+               for name in reactions]
 
         # Create 3D array for microscopic cross sections
         microxs_arr = np.zeros((len(nuclides), len(mts), 1))
@@ -434,6 +523,10 @@ def from_multigroup_flux(
         # Normalize multigroup flux
         multigroup_flux /= flux_sum
 
+        # If nu-fission was requested, get paths to pointwise data files
+        if 'nu-fission' in reactions:
+            data_library = DataLibrary.from_xml(cross_sections)
+
         # Compute microscopic cross sections within a temporary session
         with openmc.lib.TemporarySession(**init_kwargs):
             # For each nuclide and reaction, compute the flux-averaged xs
@@ -442,9 +535,15 @@ def from_multigroup_flux(
                     continue
                 lib_nuc = openmc.lib.load_nuclide(nuc)
                 for mt_index, mt in enumerate(mts):
-                    microxs_arr[nuc_index, mt_index, 0] = lib_nuc.collapse_rate(
-                        mt, temperature, energies, multigroup_flux
-                    )
+                    if mt is None:
+                        path = data_library.get_by_material(nuc)['path']
+                        microxs_arr[nuc_index, mt_index, 0] = \
+                            _collapse_nu_fission(path, nuc, temperature,
+                                                 energies, multigroup_flux)
+                    else:
+                        microxs_arr[nuc_index, mt_index, 0] = \
+                            lib_nuc.collapse_rate(
+                                mt, temperature, energies, multigroup_flux)
 
         return cls(microxs_arr, nuclides, reactions)