Enabling axial expansion with detailed depletion (#1954)

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,

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(

armi/reactor/composites.py

@@ -1575,6 +1575,9 @@ def changeNDensByFactor(self, factor):
             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):
         """

armi/reactor/converters/axialExpansionChanger/axialExpansionChanger.py

@@ -348,11 +348,7 @@ def axiallyExpandAssembly(self):
                             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 @@ def manageCoreMesh(self, r):
         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()

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.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 @@ def test_thermalExpansionContractionConservation_simple(self):
         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 @@ def test_thermalExpansionContractionConservation_simple(self):
             # 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 @@ def _checkNDens(self, prevNDen, newNDens, ratio):
                 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 @@ def _getComponentMassAndNDens(a):
                 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 @@ def setUp(self):
         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."""

armi/reactor/tests/test_components.py

@@ -697,9 +697,11 @@ def test_getNumberDensities(self):
     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()

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