Repository: terrapower/armi
Pull request #1954 (2024-10-15 21:02:18)

Enabling axial expansion with detailed depletion
## What is the change?

Detailed depletion calculations currently do not allow for axial expansion as the untruncated isotope vector is saved at the block level. Since axial expansion needs the fuel component isotopic composition, some changes are made in the framework to support calculating and storing detailed depletion and the isotope vector at the component level.

- Add component level parameters `detailedNDens`, `powerDecay`, and `pdensDecay` 
- Add functionality and testing to update component `detailedNDens` to `changeNDensByFactor` 

## Why is the change being made?
To support downstream depletion upgrades.
 
<!-- MANDATORY: Explain why the change is necessary -->
<!-- Optional: Link to any related GitHub Issues -->

---

## 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`.

onufer: this looks good to me, but i think we need @albeanth  to confirm that he thinks it will hit all his axial expansion corner cases
john-science: bump

Patch:
diff --git a/armi/physics/neutronics/parameters.py b/armi/physics/neutronics/parameters.py
index 78268b8b..fc84a5df 100644
--- a/armi/physics/neutronics/parameters.py
+++ b/armi/physics/neutronics/parameters.py
@@ -655,8 +655,6 @@ def _getNeutronicsBlockParams():
             categories=[parameters.Category.neutronics],
         )
 
-        pb.defParam("powerDecay", units=units.WATTS, description="Total decay power")
-
         pb.defParam(
             "powerGamma",
             units=units.WATTS,
diff --git a/armi/reactor/components/componentParameters.py b/armi/reactor/components/componentParameters.py
index 14292b12..92bde7dc 100644
--- a/armi/reactor/components/componentParameters.py
+++ b/armi/reactor/components/componentParameters.py
@@ -15,6 +15,7 @@
 """Component parameter definitions."""
 from armi.reactor import parameters
 from armi.reactor.parameters import ParamLocation
+from armi.reactor.parameters.parameterDefinitions import isNumpyArray
 from armi.utils import units
 
 
@@ -35,7 +36,7 @@ def getComponentParameterDefinitions():
         pb.defParam(
             "mult",
             units=units.UNITLESS,
-            description="The multiplicity of this component, i.e. how many of them there are. ",
+            description="The multiplicity of this component, i.e. how many of them there are.",
             default=1,
         )
 
@@ -63,6 +64,20 @@ def getComponentParameterDefinitions():
             description="Number densities of each nuclide.",
         )
 
+        pb.defParam(
+            "detailedNDens",
+            setter=isNumpyArray("detailedNDens"),
+            units=f"atoms/(bn*{units.CM})",
+            description=(
+                "High-fidelity number density vector with up to thousands of nuclides. "
+                "Used in high-fi depletion runs where low-fi depletion may also be occurring. "
+                "This param keeps the hi-fi and low-fi depletion values from interfering. "
+                "See core.p.detailedNucKeys for keys."
+            ),
+            saveToDB=True,
+            default=None,
+        )
+
         pb.defParam(
             "percentBu",
             units=f"{units.PERCENT_FIMA}",
@@ -98,7 +113,7 @@ def getComponentParameterDefinitions():
         pb.defParam(
             "customIsotopicsName",
             units=units.UNITLESS,
-            description="Label of isotopics applied to this component. ",
+            description="Label of isotopics applied to this component.",
         )
 
         pb.defParam(
@@ -111,7 +126,7 @@ def getComponentParameterDefinitions():
         pb.defParam(
             "zrFrac",
             units=units.UNITLESS,
-            description="Original Zr frac of this, used for material properties. ",
+            description="Original Zr frac of this, used for material properties.",
         )
 
         pb.defParam(
diff --git a/armi/reactor/composites.py b/armi/reactor/composites.py
index 3ae44b78..7f09d865 100644
--- a/armi/reactor/composites.py
+++ b/armi/reactor/composites.py
@@ -1575,6 +1575,9 @@ class ArmiObject(metaclass=CompositeModelType):
             nuc: val * factor for nuc, val in self.getNumberDensities().items()
         }
         self.setNumberDensities(densitiesScaled)
+        # Update detailedNDens
+        if self.p.detailedNDens is not None:
+            self.p.detailedNDens *= factor
 
     def clearNumberDensities(self):
         """
diff --git a/armi/reactor/converters/axialExpansionChanger/axialExpansionChanger.py b/armi/reactor/converters/axialExpansionChanger/axialExpansionChanger.py
index 86f88187..b38a7722 100644
--- a/armi/reactor/converters/axialExpansionChanger/axialExpansionChanger.py
+++ b/armi/reactor/converters/axialExpansionChanger/axialExpansionChanger.py
@@ -348,11 +348,7 @@ class AxialExpansionChanger:
                             c.zbottom = self.linked.linkedBlocks[b].lower.p.ztop
                     c.ztop = c.zbottom + c.height
                     # update component number densities
-                    newNumberDensities = {
-                        nuc: c.getNumberDensity(nuc) / growFrac
-                        for nuc in c.getNuclides()
-                    }
-                    c.setNumberDensities(newNumberDensities)
+                    c.changeNDensByFactor(1.0 / growFrac)
                     # redistribute block boundaries if on the target component
                     if self.expansionData.isTargetComponent(c):
                         b.p.ztop = c.ztop
@@ -393,7 +389,7 @@ class AxialExpansionChanger:
         if not self._detailedAxialExpansion:
             # loop through again now that the reference is adjusted and adjust the non-fuel assemblies.
             for a in r.core.getAssemblies():
-                a.setBlockMesh(r.core.refAssem.getAxialMesh())
+                a.setBlockMesh(r.core.refAssem.getAxialMesh(), conserveMassFlag="auto")
 
         oldMesh = r.core.p.axialMesh
         r.core.updateAxialMesh()
diff --git a/doc/release/0.4.rst b/doc/release/0.4.rst
index edda7276..056fc906 100644
--- a/doc/release/0.4.rst
+++ b/doc/release/0.4.rst
@@ -27,6 +27,7 @@ New Features
 #. Improve efficiency of reaction rate calculations. (`PR#1887 <https://github.com/terrapower/armi/pull/1887>`_)
 #. Adding new options for simplifying 1D cross section modeling. (`PR#1949 <https://github.com/terrapower/armi/pull/1949>`_)
 #. Updating ``copyOrWarn`` and ``getFileSHA1Hash`` to support directories. (`PR#1984 <https://github.com/terrapower/armi/pull/1984>`_)
+#. Exposing ``detailedNDens`` to components. (`PR#1954 <https://github.com/terrapower/armi/pull/1954>`_)
 #. TBD
 
 API Changes


Test patch:
diff --git a/.github/workflows/unittests.yaml b/.github/workflows/unittests.yaml
index 3e0b66db..a7271814 100644
--- a/.github/workflows/unittests.yaml
+++ b/.github/workflows/unittests.yaml
@@ -38,6 +38,6 @@ jobs:
         run: |
           pip install -e .[memprof,mpi,test]
           pytest -n 4 armi
-          mpiexec -n 2 --use-hwthread-cpus coverage run --rcfile=pyproject.toml -m pytest --cov=armi --cov-config=pyproject.toml --ignore=venv armi/tests/test_mpiFeatures.py || true
-          mpiexec -n 2 --use-hwthread-cpus coverage run --rcfile=pyproject.toml -m pytest --cov=armi --cov-config=pyproject.toml --ignore=venv armi/tests/test_mpiParameters.py || true
-          mpiexec -n 2 --use-hwthread-cpus coverage run --rcfile=pyproject.toml -m pytest --cov=armi --cov-config=pyproject.toml --ignore=venv armi/utils/tests/test_directoryChangersMpi.py || true
+          mpiexec -n 2 --use-hwthread-cpus pytest armi/tests/test_mpiFeatures.py
+          mpiexec -n 2 --use-hwthread-cpus pytest armi/tests/test_mpiParameters.py
+          mpiexec -n 2 --use-hwthread-cpus pytest armi/utils/tests/test_directoryChangersMpi.py
diff --git a/armi/bookkeeping/db/tests/test_comparedb3.py b/armi/bookkeeping/db/tests/test_comparedb3.py
index a19b9495..2b07a3dc 100644
--- a/armi/bookkeeping/db/tests/test_comparedb3.py
+++ b/armi/bookkeeping/db/tests/test_comparedb3.py
@@ -228,15 +228,17 @@ class TestCompareDB3(unittest.TestCase):
             refData4 = f4.create_dataset("numberDensities", data=a2)
             refData4.attrs["shapes"] = "2"
             refData4.attrs["numDens"] = a2
+            refData4.attrs["specialFormatting"] = True
             f5 = h5py.File("test_diffSpecialData5.hdf5", "w")
             srcData5 = f5.create_dataset("numberDensities", data=a2)
             srcData5.attrs["shapes"] = "2"
             srcData5.attrs["numDens"] = a2
+            srcData5.attrs["specialFormatting"] = True
 
-            # there should an exception
-            with self.assertRaises(Exception) as e:
+            # there should a log message
+            with mockRunLogs.BufferLog() as mock:
                 _diffSpecialData(refData4, srcData5, out, dr)
-                self.assertIn("Unable to unpack special data for paramName", e)
+                self.assertIn("Unable to unpack special data for", mock.getStdout())
 
             # make an H5 datasets that will add a np.inf diff because keys don't match
             f6 = h5py.File("test_diffSpecialData6.hdf5", "w")
diff --git a/armi/reactor/converters/tests/test_axialExpansionChanger.py b/armi/reactor/converters/tests/test_axialExpansionChanger.py
index 98cc6ea0..a8a4708b 100644
--- a/armi/reactor/converters/tests/test_axialExpansionChanger.py
+++ b/armi/reactor/converters/tests/test_axialExpansionChanger.py
@@ -14,6 +14,7 @@
 
 """Test axialExpansionChanger."""
 import collections
+import copy
 import os
 import unittest
 from statistics import mean
@@ -28,15 +29,15 @@ from armi.reactor.components import Component, DerivedShape, UnshapedComponent
 from armi.reactor.components.basicShapes import Circle, Hexagon, Rectangle
 from armi.reactor.components.complexShapes import Helix
 from armi.reactor.converters.axialExpansionChanger import (
-    AxialExpansionChanger,
     AssemblyAxialLinkage,
+    AxialExpansionChanger,
     ExpansionData,
     getSolidComponents,
     iterSolidComponents,
 )
 from armi.reactor.converters.axialExpansionChanger.assemblyAxialLinkage import (
-    areAxiallyLinked,
     AxialLink,
+    areAxiallyLinked,
 )
 from armi.reactor.flags import Flags
 from armi.reactor.tests.test_reactors import loadTestReactor, reduceTestReactorRings
@@ -270,6 +271,7 @@ class TestConservation(AxialExpansionTestBase, unittest.TestCase):
         a = buildTestAssemblyWithFakeMaterial(name="HT9")
         origMesh = a.getAxialMesh()[:-1]
         origMasses, origNDens = self._getComponentMassAndNDens(a)
+        origDetailedNDens = self._setComponentDetailedNDens(a, origNDens)
         axialExpChngr = AxialExpansionChanger(detailedAxialExpansion=True)
 
         tempGrid = linspace(0.0, a.getHeight())
@@ -284,16 +286,20 @@ class TestConservation(AxialExpansionTestBase, unittest.TestCase):
             # Set new isothermal temp and expand
             tempField = array([temp] * len(tempGrid))
             oldMasses, oldNDens = self._getComponentMassAndNDens(a)
+            oldDetailedNDens = self._getComponentDetailedNDens(a)
             axialExpChngr.performThermalAxialExpansion(a, tempGrid, tempField)
             newMasses, newNDens = self._getComponentMassAndNDens(a)
+            newDetailedNDens = self._getComponentDetailedNDens(a)
             self._checkMass(oldMasses, newMasses)
             self._checkNDens(oldNDens, newNDens, totGrowthFrac)
+            self._checkDetailedNDens(oldDetailedNDens, newDetailedNDens, totGrowthFrac)
 
         # make sure that the assembly returned to the original state
         for orig, new in zip(origMesh, a.getAxialMesh()):
             self.assertAlmostEqual(orig, new, places=12)
         self._checkMass(origMasses, newMasses)
         self._checkNDens(origNDens, newNDens, 1.0)
+        self._checkDetailedNDens(origDetailedNDens, newDetailedNDens, 1.0)
 
     def test_thermalExpansionContractionConservation_complex(self):
         """Thermally expand and then contract to ensure original state is recovered.
@@ -416,6 +422,17 @@ class TestConservation(AxialExpansionTestBase, unittest.TestCase):
                 if prev:
                     self.assertAlmostEqual(prev / new, ratio, msg=f"{prev} / {new}")
 
+    def _checkDetailedNDens(self, prevDetailedNDen, newDetailedNDens, ratio):
+        """Check whether the detailedNDens of two input dictionaries containing the
+        detailedNDens arrays for all components of an assembly are conserved.
+        """
+        for prevComp, newComp in zip(
+            prevDetailedNDen.values(), newDetailedNDens.values()
+        ):
+            for prev, new in zip(prevComp, newComp):
+                if prev:
+                    self.assertAlmostEqual(prev / new, ratio, msg=f"{prev} / {new}")
+
     @staticmethod
     def _getComponentMassAndNDens(a):
         masses = {}
@@ -426,6 +443,30 @@ class TestConservation(AxialExpansionTestBase, unittest.TestCase):
                 nDens[c] = c.getNumberDensities()
         return masses, nDens
 
+    @staticmethod
+    def _setComponentDetailedNDens(a, nDens):
+        """Returns a dictionary that contains detailedNDens for all components in an
+        assembly object input which are set to the corresponding component number densities
+        from a number density dictionary input.
+        """
+        detailedNDens = {}
+        for b in a:
+            for c in getSolidComponents(b):
+                c.p.detailedNDens = copy.deepcopy([val for val in nDens[c].values()])
+                detailedNDens[c] = c.p.detailedNDens
+        return detailedNDens
+
+    @staticmethod
+    def _getComponentDetailedNDens(a):
+        """Returns a dictionary containing all solid components and their corresponding
+        detailedNDens from an assembly object input.
+        """
+        detailedNDens = {}
+        for b in a:
+            for c in getSolidComponents(b):
+                detailedNDens[c] = copy.deepcopy(c.p.detailedNDens)
+        return detailedNDens
+
     def test_targetComponentMassConservation(self):
         """Tests mass conservation for target components."""
         self.expandAssemForMassConservationTest()
@@ -571,20 +612,65 @@ class TestManageCoreMesh(unittest.TestCase):
         reduceTestReactorRings(self.r, o.cs, 3)
 
         self.oldAxialMesh = self.r.core.p.axialMesh
+        self.componentLst = []
+        for b in self.r.core.refAssem:
+            if b.hasFlags([Flags.FUEL, Flags.PLENUM]):
+                self.componentLst.extend(getSolidComponents(b))
         # expand refAssem by 1.01 L1/L0
-        componentLst = [c for b in self.r.core.refAssem for c in b]
-        expansionGrowthFracs = 1.01 + zeros(len(componentLst))
+        expansionGrowthFracs = 1.01 + zeros(len(self.componentLst))
+        (
+            self.origDetailedNDens,
+            self.origVolumes,
+        ) = self._getComponentDetailedNDensAndVol(self.componentLst)
         self.axialExpChngr.performPrescribedAxialExpansion(
-            self.r.core.refAssem, componentLst, expansionGrowthFracs, setFuel=True
+            self.r.core.refAssem, self.componentLst, expansionGrowthFracs, setFuel=True
         )
 
     def test_manageCoreMesh(self):
         self.axialExpChngr.manageCoreMesh(self.r)
         newAxialMesh = self.r.core.p.axialMesh
-        # skip first and last entries as they do not change
-        for old, new in zip(self.oldAxialMesh[1:-1], newAxialMesh[1:-1]):
+        # the top and bottom and top of the grid plate block are not expected to change
+        for old, new in zip(self.oldAxialMesh[2:-1], newAxialMesh[2:-1]):
             self.assertLess(old, new)
 
+    def test_componentConservation(self):
+        self.axialExpChngr.manageCoreMesh(self.r)
+        newDetailedNDens, newVolumes = self._getComponentDetailedNDensAndVol(
+            self.componentLst
+        )
+        for c in newVolumes.keys():
+            self._checkMass(
+                self.origDetailedNDens[c],
+                self.origVolumes[c],
+                newDetailedNDens[c],
+                newVolumes[c],
+                c,
+            )
+
+    def _getComponentDetailedNDensAndVol(self, componentLst):
+        """Returns a tuple containing dictionaries of detailedNDens and volumes of
+        all components from a component list input.
+        """
+        detailedNDens = {}
+        volumes = {}
+        for c in componentLst:
+            c.p.detailedNDens = [val for val in c.getNumberDensities().values()]
+            detailedNDens[c] = copy.deepcopy(c.p.detailedNDens)
+            volumes[c] = c.getVolume()
+        return (detailedNDens, volumes)
+
+    def _checkMass(self, origDetailedNDens, origVolume, newDetailedNDens, newVolume, c):
+        for prevMass, newMass in zip(
+            origDetailedNDens * origVolume, newDetailedNDens * newVolume
+        ):
+            if c.parent.hasFlags(Flags.FUEL):
+                self.assertAlmostEqual(
+                    prevMass, newMass, delta=1e-12, msg=f"{c}, {c.parent}"
+                )
+            else:
+                # should not conserve mass here as it is structural material above active fuel
+                self.assertAlmostEqual(newMass / prevMass, 0.99, msg=f"{c}, {c.parent}")
+
 
 class TestExceptions(AxialExpansionTestBase, unittest.TestCase):
     """Verify exceptions are caught."""
diff --git a/armi/reactor/tests/test_components.py b/armi/reactor/tests/test_components.py
index 46b93609..7cad72c9 100644
--- a/armi/reactor/tests/test_components.py
+++ b/armi/reactor/tests/test_components.py
@@ -697,9 +697,11 @@ class TestCircle(TestShapedComponent):
     def test_changeNumberDensities(self):
         """Test that demonstates that the number densities on a component can be modified."""
         self.component.p.numberDensities = {"NA23": 1.0}
+        self.component.p.detailedNDens = [1.0]
         self.assertEqual(self.component.getNumberDensity("NA23"), 1.0)
         self.component.changeNDensByFactor(3.0)
         self.assertEqual(self.component.getNumberDensity("NA23"), 3.0)
+        self.assertEqual(self.component.p.detailedNDens[0], 3.0)
 
     def test_fuelMass(self):
         nominalMass = self.component.getMass()
diff --git a/armi/utils/tests/test_pathTools.py b/armi/utils/tests/test_pathTools.py
index 2b6f046d..8ac3affe 100644
--- a/armi/utils/tests/test_pathTools.py
+++ b/armi/utils/tests/test_pathTools.py
@@ -19,6 +19,7 @@ import types
 import unittest
 
 from armi import context
+from armi.tests import mockRunLogs
 from armi.utils import pathTools
 from armi.utils.directoryChangers import TemporaryDirectoryChanger
 
@@ -36,11 +37,10 @@ class PathToolsTests(unittest.TestCase):
             pathTools.copyOrWarn("Test File", path, pathCopy)
             self.assertTrue(os.path.exists(pathCopy))
 
-            # Test an exception
-            with self.assertRaises(Exception) as e:
+            # Test a non-existant file
+            with mockRunLogs.BufferLog() as mock:
                 pathTools.copyOrWarn("Test File", "FileDoesntExist.txt", pathCopy)
-                self.assertIn("Could not copy", e)
-                self.assertIn("No such file or directory", e)
+                self.assertIn("Could not copy", mock.getStdout())
 
     def test_copyOrWarnDir(self):
         with TemporaryDirectoryChanger():
@@ -55,11 +55,10 @@ class PathToolsTests(unittest.TestCase):
             self.assertTrue(os.path.exists(pathDirCopy))
             self.assertTrue(os.path.exists(os.path.join(pathDirCopy, "test.txt")))
 
-            # Test an exception
-            with self.assertRaises(Exception) as e:
+            # Test a non-existant file
+            with mockRunLogs.BufferLog() as mock:
                 pathTools.copyOrWarn("Test File", "DirDoesntExist", pathDirCopy)
-                self.assertIn("Could not copy", e)
-                self.assertIn("No such file or directory", e)
+                self.assertIn("Could not copy", mock.getStdout())
 
     def test_separateModuleAndAttribute(self):
         self.assertRaises(