Refactor MicroXS class and usage in IndependentOperator

docs/source/pythonapi/deplete.rst

@@ -57,6 +57,16 @@ The :class:`CoupledOperator` and :class:`IndependentOperator` classes must also
 have some knowledge of how nuclides transmute and decay. This is handled by the
 :class:`Chain` class.
 
+The :class:`IndependentOperator` class requires a set of fluxes and microscopic
+cross sections. The following function can be used to generate this information:
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: myfunction.rst
+
+   get_microxs_and_flux
+
 Minimal Example
 ---------------
 

docs/source/usersguide/depletion.rst

@@ -197,43 +197,54 @@ across all material instances.
 Transport-independent depletion
 ===============================
 
-.. warning::
-
-   This feature is still under heavy development and has yet to be rigorously
-   verified. API changes and feature additions are possible and likely in
-   the near future.
-
-This category of operator uses one-group microscopic cross sections to obtain
-transmutation reaction rates. The cross sections are pre-calculated, so there is
-no need for direct coupling between a transport-independent operator and a
-transport solver. The :mod:`openmc.deplete` module offers a single
-transport-independent operator, :class:`~openmc.deplete.IndependentOperator`,
-and only one operator is needed since, in theory, any transport code could
-calcuate the one-group microscopic cross sections.
+This category of operator uses multigroup microscopic cross sections along with
+multigroup flux spectra to obtain transmutation reaction rates. The cross
+sections are pre-calculated, so there is no need for direct coupling between a
+transport-independent operator and a transport solver. The :mod:`openmc.deplete`
+module offers a single transport-independent operator,
+:class:`~openmc.deplete.IndependentOperator`, and only one operator is needed
+since, in theory, any transport code could calculate the multigroup microscopic
+cross sections. The :class:`~openmc.deplete.IndependentOperator` class has two
+constructors. The default constructor requires a :class:`openmc.Materials`
+instance, a list of multigroup flux arrays, and a list of
+:class:`~openmc.deplete.MicroXS` instances containing multigroup microscopic
+cross sections in units of barns. This might look like the following::
+
+    materials = openmc.Materials([m1, m2, m3])
+    ...
 
-The :class:`~openmc.deplete.IndependentOperator` class has two constructors.
-The default constructor requires a :class:`openmc.Materials` instance, a
-:class:`~openmc.deplete.MicroXS` instance containing one-group microscoic cross
-sections in units of barns, and a path to a depletion chain file::
+    # Assign fluxes (generated from any code)
+    flux_m1 = numpy.array([...])
+    flux_m2 = numpy.array([...])
+    flux_m3 = numpy.array([...])
+    fluxes = [flux_m1, flux_m2, flux_m3]
 
-    materials = openmc.Materials()
-    ...
+    # Assign microscopic cross sections
+    micro_m1 = openmc.deplete.MicroXS.from_csv('xs_m1.csv')
+    micro_m2 = openmc.deplete.MicroXS.from_csv('xs_m2.csv')
+    micro_m3 = openmc.deplete.MicroXS.from_csv('xs_m3.csv')
+    micros = [micro_m1, micro_m2, micro_m3]
 
-    # load in the microscopic cross sections
-    micro_xs = openmc.deplete.MicroXS.from_csv(micro_xs_path)
+    # Create operator
+    op = openmc.deplete.IndependentOperator(materials, fluxes, micros)
 
-    op = openmc.deplete.IndependentOperator(materials, micro_xs, chain_file)
+For more details on the :class:`~openmc.deplete.MicroXS` class, including how to
+use OpenMC's transport solver to generate microscopic cross sections and fluxes
+for use with :class:`~openmc.deplete.IndependentOperator`, see :ref:`micros`.
 
 .. note::
 
-   The same statements from :ref:`coupled-depletion` about which
-   materials are depleted and the requirement for depletable materials to have
-   a specified volume also apply here.
+   The same statements from :ref:`coupled-depletion` about which materials are
+   depleted and the requirement for depletable materials to have a specified
+   volume also apply here.
 
 An alternate constructor,
 :meth:`~openmc.deplete.IndependentOperator.from_nuclides`, accepts a volume and
 dictionary of nuclide concentrations in place of the :class:`openmc.Materials`
-instance::
+instance. Note that while the normal constructor allows multiple materials to be
+depleted with a single operator, the
+:meth:`~openmc.deplete.IndependentOperator.from_nuclides` classmethod only works
+for a single material::
 
     nuclides = {'U234': 8.92e18,
                 'U235': 9.98e20,
@@ -244,6 +255,7 @@ instance::
     volume = 0.5
     op = openmc.deplete.IndependentOperator.from_nuclides(volume,
                                                           nuclides,
+                                                          flux,
                                                           micro_xs,
                                                           chain_file,
                                                           nuc_units='atom/cm3')
@@ -253,18 +265,21 @@ transport-depletion calculation and follow the same steps from there.
 
 .. note::
 
-   Ideally, one-group cross section data should be available for every
-   reaction in the depletion chain. If cross section data is not present for
-   a nuclide in the depletion chain with at least one reaction, that reaction
-   will not be simulated.
+   Ideally, multigroup cross section data should be available for every reaction
+   in the depletion chain. If cross section data is not present for a nuclide in
+   the depletion chain with at least one reaction, that reaction will not be
+   simulated.
+
+.. _micros:
 
 Loading and Generating Microscopic Cross Sections
 -------------------------------------------------
 
-As mentioned earlier, any transport code could be used to calculate one-group
-microscopic cross sections. The :mod:`openmc.deplete` module provides the
-:class:`~openmc.deplete.MicroXS` class, which contains methods to read in
-pre-calculated cross sections from a ``.csv`` file or from data arrays::
+As mentioned above, any transport code could be used to calculate multigroup
+microscopic cross sections and fluxes. The :mod:`openmc.deplete` module provides
+the :class:`~openmc.deplete.MicroXS` class, which can either be instantiated
+from pre-calculated cross sections in a ``.csv`` file or from data arrays
+directly::
 
     micro_xs = MicroXS.from_csv(micro_xs_path)
 
@@ -273,37 +288,31 @@ pre-calculated cross sections from a ``.csv`` file or from data arrays::
     data = np.array([[0.1, 0.2],
                      [0.3, 0.4],
                      [0.01, 0.5]])
-    micro_xs = MicroXS.from_array(nuclides, reactions, data)
+    micro_xs = MicroXS(data, nuclides, reactions)
 
 .. important::
 
-   Both :meth:`~openmc.deplete.MicroXS.from_csv()` and
-   :meth:`~openmc.deplete.MicroXS.from_array()` assume the cross section values
-   provided are in barns by defualt, but have no way of verifying this. Make
-   sure your cross sections are in the correct units before passing to a
+   The cross section values are assumed to be in units of barns. Make sure your
+   cross sections are in the correct units before passing to a
    :class:`~openmc.deplete.IndependentOperator` object.
 
-The :class:`~openmc.deplete.MicroXS` class also contains a method to generate one-group microscopic cross sections using OpenMC's transport solver. The
-:meth:`~openmc.deplete.MicroXS.from_model()` method will produce a
-:class:`~openmc.deplete.MicroXS` instance with microscopic cross section data in
-units of barns::
-
-    import openmc
+Additionally, a convenience function,
+:func:`~openmc.deplete.get_microxs_and_flux`, can provide the needed fluxes and
+cross sections using OpenMC's transport solver::
 
-    model = openmc.Model.from_xml()
+    model = openmc.Model()
+    ...
 
-    micro_xs = openmc.deplete.MicroXS.from_model(model,
-                                                 model.materials[0],
-                                                 chain_file)
+    fluxes, micros = openmc.deplete.get_microxs_and_flux(model, materials)
 
-If you are running :meth:`~openmc.deplete.MicroXS.from_model()` on a cluster
+If you are running :func:`~openmc.deplete.get_microxs_and_flux` on a cluster
 where temporary files are created on a local filesystem that is not shared
 across nodes, you'll need to set an environment variable pointing to a local
 directoy so that each MPI process knows where to store output files used to
 calculate the microscopic cross sections. In order of priority, they are
-:envvar:`TMPDIR`. :envvar:`TEMP`, and :envvar:`TMP`. Users interested in
-further details can read the documentation for the `tempfile <https://docs.python.org/3/library/tempfile.html#tempfile.gettempdir>`_ module.
-
+:envvar:`TMPDIR`. :envvar:`TEMP`, and :envvar:`TMP`. Users interested in further
+details can read the documentation for the `tempfile
+<https://docs.python.org/3/library/tempfile.html#tempfile.gettempdir>`_ module.
 
 Caveats
 -------

openmc/deplete/helpers.py

@@ -196,13 +196,13 @@ def reset_tally_means(self):
         """
         self._rate_tally_means_cache = None
 
-    def get_material_rates(self, mat_id, nuc_index, rx_index):
+    def get_material_rates(self, mat_index, nuc_index, rx_index):
         """Return an array of reaction rates for a material
 
         Parameters
         ----------
-        mat_id : int
-            Unique ID for the requested material
+        mat_index : int
+            Index for the material
         nuc_index : iterable of int
             Index for each nuclide in :attr:`nuclides` in the
             desired reaction rate matrix
@@ -216,7 +216,7 @@ def get_material_rates(self, mat_id, nuc_index, rx_index):
             reaction rates in this material
         """
         self._results_cache.fill(0.0)
-        full_tally_res = self.rate_tally_means[mat_id]
+        full_tally_res = self.rate_tally_means[mat_index]
         for i_tally, (i_nuc, i_rx) in enumerate(product(nuc_index, rx_index)):
             self._results_cache[i_nuc, i_rx] = full_tally_res[i_tally]
 

openmc/deplete/independent_operator.py

@@ -5,8 +5,10 @@
 
 """
 
+from __future__ import annotations
+from collections.abc import Iterable
 import copy
-from itertools import product
+from typing import List, Set
 
 import numpy as np
 from uncertainties import ufloat
@@ -37,14 +39,20 @@ class IndependentOperator(OpenMCOperator):
 
     .. versionadded:: 0.13.1
 
+    .. versionchanged:: 0.13.4
+        Arguments updated to include list of fluxes and microscopic cross
+        sections.
+
     Parameters
     ----------
     materials : openmc.Materials
         Materials to deplete.
-    micro_xs : MicroXS
-        One-group microscopic cross sections in [b]. If the
-        :class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation will
-        be run.
+    fluxes : list of numpy.ndarray
+        Flux in each group in [n-cm/src] for each domain
+    micros : list of MicroXS
+        Cross sections in [b] for each domain. If the
+        :class:`~openmc.deplete.MicroXS` object is empty, a decay-only
+        calculation will be run.
     chain_file : str
         Path to the depletion chain XML file. Defaults to
         ``openmc.config['chain_file']``.
@@ -109,7 +117,8 @@ class IndependentOperator(OpenMCOperator):
 
     def __init__(self,
                  materials,
-                 micro_xs,
+                 fluxes,
+                 micros,
                  chain_file=None,
                  keff=None,
                  normalization_mode='fission-q',
@@ -120,7 +129,7 @@ def __init__(self,
                  fission_yield_opts=None):
         # Validate micro-xs parameters
         check_type('materials', materials, openmc.Materials)
-        check_type('micro_xs', micro_xs, MicroXS)
+        check_type('micros', micros, Iterable, MicroXS)
         if keff is not None:
             check_type('keff', keff, tuple, float)
             keff = ufloat(*keff)
@@ -132,10 +141,15 @@ def __init__(self,
         helper_kwargs = {'normalization_mode': normalization_mode,
                          'fission_yield_opts': fission_yield_opts}
 
-        cross_sections = micro_xs * 1e-24
+        # Sort fluxes and micros in same order that materials get sorted
+        index_sort = np.argsort([mat.id for mat in materials])
+        fluxes = [fluxes[i] for i in index_sort]
+        micros = [micros[i] for i in index_sort]
+
+        self.fluxes = fluxes
         super().__init__(
             materials,
-            cross_sections,
+            micros,
             chain_file,
             prev_results,
             fission_q=fission_q,
@@ -145,6 +159,7 @@ def __init__(self,
 
     @classmethod
     def from_nuclides(cls, volume, nuclides,
+                      flux,
                       micro_xs,
                       chain_file=None,
                       nuc_units='atom/b-cm',
@@ -163,10 +178,11 @@ def from_nuclides(cls, volume, nuclides,
         nuclides : dict of str to float
             Dictionary with nuclide names as keys and nuclide concentrations as
             values.
+        flux : numpy.ndarray
+            Flux in each group in [n-cm/src]
         micro_xs : MicroXS
-            One-group microscopic cross sections in [b]. If the
-            :class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation
-            will be run.
+            Cross sections in [b]. If the :class:`~openmc.deplete.MicroXS`
+            object is empty, a decay-only calculation will be run.
         chain_file : str, optional
             Path to the depletion chain XML file. Defaults to
             ``openmc.config['chain_file']``.
@@ -203,8 +219,11 @@ def from_nuclides(cls, volume, nuclides,
         """
         check_type('nuclides', nuclides, dict, str)
         materials = cls._consolidate_nuclides_to_material(nuclides, nuc_units, volume)
+        fluxes = [flux]
+        micros = [micro_xs]
         return cls(materials,
-                   micro_xs,
+                   fluxes,
+                   micros,
                    chain_file,
                    keff=keff,
                    normalization_mode=normalization_mode,
@@ -256,9 +275,9 @@ def _load_previous_results(self):
                 self.prev_res.append(new_res)
 
 
-    def _get_nuclides_with_data(self, cross_sections):
+    def _get_nuclides_with_data(self, cross_sections: List[MicroXS]) -> Set[str]:
         """Finds nuclides with cross section data"""
-        return set(cross_sections.index)
+        return set(cross_sections[0].nuclides)
 
     class _IndependentRateHelper(ReactionRateHelper):
         """Class for generating one-group reaction rates with flux and
@@ -285,7 +304,7 @@ class _IndependentRateHelper(ReactionRateHelper):
 
         """
 
-        def __init__(self, op):
+        def __init__(self, op: IndependentOperator):
             rates = op.reaction_rates
             super().__init__(rates.n_nuc, rates.n_react)
 
@@ -301,34 +320,32 @@ def reset_tally_means(self):
             """Unused in this case"""
             pass
 
-        def get_material_rates(self, mat_id, nuc_index, react_index):
+        def get_material_rates(self, mat_index, nuc_index, react_index):
             """Return 2D array of [nuclide, reaction] reaction rates
 
             Parameters
             ----------
-            mat_id : int
-                Unique ID for the requested material
+            mat_index : int
+                Index for the material
             nuc_index : list of str
                 Ordering of desired nuclides
             react_index : list of str
                 Ordering of reactions
             """
             self._results_cache.fill(0.0)
 
-            # Get volume in units of [b-cm]
-            volume_b_cm = 1e24 * self._op.number.get_mat_volume(mat_id)
+            # Get flux and microscopic cross sections from operator
+            flux = self._op.fluxes[mat_index]
+            xs = self._op.cross_sections[mat_index]
 
-            for i_nuc, i_react in product(nuc_index, react_index):
+            for i_nuc in nuc_index:
                 nuc = self.nuc_ind_map[i_nuc]
-                rx = self.rx_ind_map[i_react]
-
-                # OK, this is kind of weird, but we multiply by volume in [b-cm]
-                # only because OpenMCOperator._calculate_reaction_rates has to
-                # divide it out later. It might make more sense to account for
-                # the source rate (flux) here rather than in the normalization
-                # helper.
-                self._results_cache[i_nuc, i_react] = \
-                    self._op.cross_sections[rx][nuc] * volume_b_cm
+                for i_rx in react_index:
+                    rx = self.rx_ind_map[i_rx]
+
+                    # Determine reaction rate by multiplying xs in [b] by flux
+                    # in [n-cm/src] to give [(reactions/src)*b-cm/atom]
+                    self._results_cache[i_nuc, i_rx] = (xs[nuc, rx] * flux).sum()
 
             return self._results_cache