Repository: terrapower/armi
Pull request #1887 (2024-09-16 17:18:53)

Improving efficiency of reaction rate calcs
## What is the change?

Improve the efficiency of reaction rate calculations.

## Why is the change being made?

A common post-processing step for a `GlobalFluxInterface` is to compute reaction rates from the flux solution and the microscopic cross sections that exist on the reactor model. The implementation of this calculation is done in a straightforward way, but it is naive with-respect-to the typical data structure in an ARMI simulation. The efficiency for typical problems can be greatly improved by refactoring the loops so that the outermost loop is over cross section type.

Often, there are thousands of blocks in a reactor model but only a handful of unique microscopic cross section types; each set of microscopic cross sections might be shared across several hundred blocks. An outer iteration over cross section types with an inner iteration over blocks can significantly reduce the number of times a microscopic cross section set needs to be retrieved.

---

## Checklist

- [x] This PR has only [one purpose or idea](https://terrapower.github.io/armi/developer/tooling.html#one-idea-one-pr).
- [x] [Tests](https://terrapower.github.io/armi/developer/tooling.html#test-it) have been added/updated to verify any new/changed code.

<!-- Check the code quality -->

- [x] The code style follows [good practices](https://terrapower.github.io/armi/developer/standards_and_practices.html).
- [x] The commit message(s) follow [good practices](https://terrapower.github.io/armi/developer/tooling.html).

<!-- Check the project-level cruft -->

- [x] The [release notes](https://terrapower.github.io/armi/developer/tooling.html#add-release-notes) have been updated if necessary.
- [x] The [documentation](https://terrapower.github.io/armi/developer/tooling.html#document-it) is still up-to-date in the `doc` folder.
- [x] The dependencies are still up-to-date in `pyproject.toml`.

john-science: Just a head's up, I know you're not done. But this seems like an SCR situation.

Patch:
diff --git a/armi/nuclearDataIO/xsCollections.py b/armi/nuclearDataIO/xsCollections.py
index faf7978b..91914fcf 100644
--- a/armi/nuclearDataIO/xsCollections.py
+++ b/armi/nuclearDataIO/xsCollections.py
@@ -50,7 +50,7 @@ from armi.utils import units
 # they have some origin in the ISOTXS file format card 04 definition
 # fmt: off
 NGAMMA = "nGamma"      # radiative capture
-NAPLHA = "nalph"       # (n, alpha)
+NALPHA = "nalph"       # (n, alpha)
 NP = "np"              # (n, proton)
 ND = "nd"              # (n, deuteron)
 NT = "nt"              # (n, triton)
@@ -59,7 +59,7 @@ N2N_XS = "n2n"         # (n,2n)
 NUSIGF = "nuSigF"
 NU = "neutronsPerFission"
 # fmt: on
-CAPTURE_XS = [NGAMMA, NAPLHA, NP, ND, NT]
+CAPTURE_XS = [NGAMMA, NALPHA, NP, ND, NT]
 
 # Cross section types that are represented by 2-D matrices in the multigroup approximation
 BASIC_SCAT_MATRIX = ["elasticScatter", "inelasticScatter", "n2nScatter"]
diff --git a/armi/physics/neutronics/globalFlux/__init__.py b/armi/physics/neutronics/globalFlux/__init__.py
index 4d4496f5..98c44bab 100644
--- a/armi/physics/neutronics/globalFlux/__init__.py
+++ b/armi/physics/neutronics/globalFlux/__init__.py
@@ -13,3 +13,6 @@
 # limitations under the License.
 
 """Global flux solvers."""
+
+RX_ABS_MICRO_LABELS = ["nGamma", "fission", "nalph", "np", "nd", "nt"]
+RX_PARAM_NAMES = ["rateCap", "rateFis", "rateProdN2n", "rateProdFis", "rateAbs"]
diff --git a/armi/physics/neutronics/globalFlux/globalFluxInterface.py b/armi/physics/neutronics/globalFlux/globalFluxInterface.py
index ca97f299..9e51adab 100644
--- a/armi/physics/neutronics/globalFlux/globalFluxInterface.py
+++ b/armi/physics/neutronics/globalFlux/globalFluxInterface.py
@@ -27,6 +27,7 @@ from armi.physics import executers
 from armi.physics import neutronics
 from armi.reactor import geometry
 from armi.reactor import reactors
+from armi.physics.neutronics.globalFlux import RX_PARAM_NAMES, RX_ABS_MICRO_LABELS
 from armi.reactor.blocks import Block
 from armi.reactor.converters import geometryConverters
 from armi.reactor.converters import uniformMesh
@@ -36,9 +37,6 @@ from armi.utils import units, codeTiming, getMaxBurnSteps, getBurnSteps
 
 ORDER = interfaces.STACK_ORDER.FLUX
 
-RX_ABS_MICRO_LABELS = ["nGamma", "fission", "nalph", "np", "nd", "nt"]
-RX_PARAM_NAMES = ["rateCap", "rateFis", "rateProdN2n", "rateProdFis", "rateAbs"]
-
 
 class GlobalFluxInterface(interfaces.Interface):
     """
diff --git a/armi/reactor/converters/uniformMesh.py b/armi/reactor/converters/uniformMesh.py
index 47d67092..47f4917d 100644
--- a/armi/reactor/converters/uniformMesh.py
+++ b/armi/reactor/converters/uniformMesh.py
@@ -72,6 +72,8 @@ from armi.reactor.converters.geometryConverters import GeometryConverter
 from armi.reactor import parameters
 from armi.reactor.reactors import Reactor
 from armi.settings.fwSettings.globalSettings import CONF_UNIFORM_MESH_MINIMUM_SIZE
+from armi.physics.neutronics.globalFlux import RX_PARAM_NAMES, RX_ABS_MICRO_LABELS
+
 
 HEAVY_METAL_PARAMS = ["molesHmBOL", "massHmBOL"]
 
@@ -1100,13 +1102,42 @@ class UniformMeshGeometryConverter(GeometryConverter):
                     aDest,
                     self.paramMapper,
                     mapNumberDensities,
-                    calcReactionRates=self.calcReactionRates,
+                    calcReactionRates=False,
+                )
+
+            # If requested, the reaction rates will be calculated based on the
+            # mapped neutron flux and the XS library.
+            if self.calcReactionRates:
+                self._calculateReactionRatesEfficient(
+                    destReactor.core, sourceReactor.core.p.keff
                 )
 
         # Clear the cached data after it has been mapped to prevent issues with
         # holding on to block data long-term.
         self._cachedReactorCoreParamData = {}
 
+    @staticmethod
+    def _calculateReactionRatesEfficient(core, keff):
+        """
+        First, sort blocks into groups by XS type. Then, we just need to grab micros for each XS type once.
+
+        Iterate over list of blocks with the given XS type; calculate reaction rates for these blocks
+        """
+        xsTypeGroups = collections.defaultdict(list)
+        for b in core.getBlocks():
+            xsTypeGroups[b.getMicroSuffix()].append(b)
+
+        for xsID, blockList in xsTypeGroups.items():
+            nucSet = set()
+            for b in blockList:
+                nucSet.update(
+                    nuc for nuc, ndens in b.getNumberDensities().items() if ndens > 0.0
+                )
+            xsNucDict = {nuc: core.lib.getNuclide(nuc, xsID) for nuc in nucSet}
+            UniformMeshGeometryConverter._calcReactionRatesBlockList(
+                blockList, keff, xsNucDict
+            )
+
     @staticmethod
     def _calculateReactionRates(lib, keff, assem):
         """
@@ -1130,6 +1161,87 @@ class UniformMeshGeometryConverter(GeometryConverter):
                 continue
             globalFluxInterface.calcReactionRates(b, keff, lib)
 
+    @staticmethod
+    def _calcReactionRatesBlockList(objList, keff, xsNucDict):
+        r"""
+        Compute 1-group reaction rates for the objects in objList (usually a block).
+
+        :meta public:
+
+        .. impl:: Return the reaction rates for a given ArmiObject
+            :id: I_ARMI_FLUX_RX_RATES_BY_XS_ID
+            :implements: R_ARMI_FLUX_RX_RATES
+
+            This is an alternative implementation of :need:`I_ARMI_FLUX_RX_RATES` that
+            is more efficient when computing reaction rates for a large set of blocks
+            that share a common set of microscopic cross sections.
+
+            For more detail on the reation rate calculations, see :need:`I_ARMI_FLUX_RX_RATES`.
+
+        Parameters
+        ----------
+        objList : List[Block]
+            The list of objects to compute reaction rates on. Notionally this could be upgraded to be
+            any kind of ArmiObject but with params defined as they are it currently is only
+            implemented for a block.
+
+        keff : float
+            The keff of the core. This is required to get the neutron production rate correct
+            via the neutron balance statement (since nuSigF has a 1/keff term).
+
+        xsNucDict: Dict[str, XSNuclide]
+            Microscopic cross sections to use in computing the reaction rates. Keys are
+            nuclide names (e.g., "U235") and values are the associated XSNuclide objects
+            from the cross section library, which contain the microscopic cross section
+            data for a given nuclide in the current cross section group.
+        """
+        for obj in objList:
+            rate = collections.defaultdict(float)
+
+            numberDensities = obj.getNumberDensities()
+            try:
+                mgFlux = np.array(obj.getMgFlux())
+            except TypeError:
+                continue
+
+            for nucName, numberDensity in numberDensities.items():
+                if numberDensity == 0.0:
+                    continue
+                nucRate = collections.defaultdict(float)
+
+                micros = xsNucDict[nucName].micros
+
+                # absorption is fission + capture (no n2n here)
+                for name in RX_ABS_MICRO_LABELS:
+                    volumetricRR = numberDensity * mgFlux.dot(micros[name])
+                    nucRate["rateAbs"] += volumetricRR
+                    if name != "fission":
+                        nucRate["rateCap"] += volumetricRR
+                    else:
+                        nucRate["rateFis"] += volumetricRR
+                        # scale nu by keff.
+                        nusigmaF = micros["fission"] * micros.neutronsPerFission
+                        nucRate["rateProdFis"] += (
+                            numberDensity * mgFlux.dot(nusigmaF) / keff
+                        )
+
+                nucRate["rateProdN2n"] += 2.0 * numberDensity * mgFlux.dot(micros.n2n)
+
+                for rx in RX_PARAM_NAMES:
+                    if nucRate[rx]:
+                        rate[rx] += nucRate[rx]
+
+            for paramName, val in rate.items():
+                obj.p[paramName] = val  # put in #/cm^3/s
+
+            if rate["rateFis"] > 0.0:
+                fuelVolFrac = obj.getComponentAreaFrac(Flags.FUEL)
+                obj.p.fisDens = rate["rateFis"] / fuelVolFrac
+                obj.p.fisDensHom = rate["rateFis"]
+            else:
+                obj.p.fisDens = 0.0
+                obj.p.fisDensHom = 0.0
+
     def updateReactionRates(self):
         """
         Update reaction rates on converted assemblies.
@@ -1147,10 +1259,9 @@ class UniformMeshGeometryConverter(GeometryConverter):
                     self.convReactor.core.lib, self.convReactor.core.p.keff, assem
                 )
         else:
-            for assem in self.convReactor.core.getAssemblies():
-                self._calculateReactionRates(
-                    self.convReactor.core.lib, self.convReactor.core.p.keff, assem
-                )
+            self._calculateReactionRatesEfficient(
+                self.convReactor.core, self.convReactor.core.p.keff
+            )
 
 
 class NeutronicsUniformMeshConverter(UniformMeshGeometryConverter):
diff --git a/doc/release/0.4.rst b/doc/release/0.4.rst
index 80ef1f4a..c2f54a32 100644
--- a/doc/release/0.4.rst
+++ b/doc/release/0.4.rst
@@ -19,11 +19,11 @@ New Features
 #. Provide ``Parameter.hasCategory`` for quickly checking if a parameter is defined with a given category. (`PR#1899 <https://github.com/terrapower/armi/pull/1899>`_)
 #. Provide ``ParameterCollection.where`` for efficient iteration over parameters who's definition matches a given condition. (`PR#1899 <https://github.com/terrapower/armi/pull/1899>`_)
 #. Flags can now be defined with letters and numbers. (`PR#1966 <https://github.com/terrapower/armi/pull/1966>`_)
-#. Plugins can provide the ``getAxialExpansionChanger`` hook to customize axial expansion. (`PR#1870 <https://github.com/terrapower/armi/pull/1870`_)
 #. Plugins can provide the ``getAxialExpansionChanger`` hook to customize axial expansion. (`PR#1870 <https://github.com/terrapower/armi/pull/1870>`_)
 #. New plugin hook ``beforeReactorConstruction`` added to enable plugins to process case settings before reactor init. (`PR#1945 <https://github.com/terrapower/armi/pull/1945>`_)
 #. Provide ``Block.getInputHeight`` for determining the height of a block from blueprints. (`PR#1927 <https://github.com/terrapower/armi/pull/1927`_)
 #. Improve performance by changing the lattice physics interface so that cross sections are not updated on ``everyNode`` calls during coupled calculations (`PR#1963 <https://github.com/terrapower/armi/pull/1963>`_)
+#. Improve efficiency of reaction rate calculations. (`PR#1887 <https://github.com/terrapower/armi/pull/1887>`_)
 #. TBD
 
 API Changes


Test patch:
diff --git a/armi/reactor/converters/tests/test_uniformMesh.py b/armi/reactor/converters/tests/test_uniformMesh.py
index 11f78153..2edaf6a4 100644
--- a/armi/reactor/converters/tests/test_uniformMesh.py
+++ b/armi/reactor/converters/tests/test_uniformMesh.py
@@ -16,6 +16,7 @@ import collections
 import os
 import random
 import unittest
+import copy
 
 from armi.nuclearDataIO.cccc import isotxs
 from armi.physics.neutronics.settings import CONF_XS_KERNEL
@@ -23,6 +24,7 @@ from armi.settings.fwSettings.globalSettings import CONF_UNIFORM_MESH_MINIMUM_SI
 from armi.reactor.converters import uniformMesh
 from armi.reactor.flags import Flags
 from armi.reactor.tests import test_assemblies
+from armi.reactor.tests import test_blocks
 from armi.reactor.tests.test_reactors import loadTestReactor, reduceTestReactorRings
 from armi.tests import TEST_ROOT, ISOAA_PATH
 
@@ -498,6 +500,31 @@ class TestUniformMesh(unittest.TestCase):
                 self.assertTrue(b.p.rateCap)
 
 
+class TestCalcReationRates(unittest.TestCase):
+    def test_calcReactionRatesBlockList(self):
+        """
+        Test that the efficient reaction rate code executes and sets a param > 0.0.
+
+        .. test:: Return the reaction rates for a given list of ArmiObjects.
+            :id: T_ARMI_FLUX_RX_RATES_BY_XS_ID
+            :tests: R_ARMI_FLUX_RX_RATES
+        """
+        b = test_blocks.loadTestBlock()
+        test_blocks.applyDummyData(b)
+        self.assertAlmostEqual(b.p.rateAbs, 0.0)
+        blockList = [copy.deepcopy(b) for _i in range(3)]
+        xsID = b.getMicroSuffix()
+        xsNucDict = {nuc: b.core.lib.getNuclide(nuc, xsID) for nuc in b.getNuclides()}
+        uniformMesh.UniformMeshGeometryConverter._calcReactionRatesBlockList(
+            blockList, 1.01, xsNucDict
+        )
+        for b in blockList:
+            self.assertGreater(b.p.rateAbs, 0.0)
+            vfrac = b.getComponentAreaFrac(Flags.FUEL)
+            self.assertEqual(b.p.fisDens, b.p.rateFis / vfrac)
+            self.assertEqual(b.p.fisDensHom, b.p.rateFis)
+
+
 class TestGammaUniformMesh(unittest.TestCase):
     """
     Tests gamma uniform mesh converter.