Fix component density issue (#841)

- fixed an Issue (#819) where components initialized through the constructor would have
  different density than those created through blueprints.

armi/reactor/blueprints/blockBlueprint.py

@@ -119,16 +119,6 @@ def construct(
                 materialInput, componentDesign
             )
             c = componentDesign.construct(blueprint, filteredMaterialInput)
-            if cs["inputHeightsConsideredHot"]:
-                if "group" in c.name:
-                    for component in c:
-                        component.applyHotHeightDensityReduction()
-                        componentBlueprint.insertDepletableNuclideKeys(
-                            component, blueprint
-                        )
-                else:
-                    c.applyHotHeightDensityReduction()
-                    componentBlueprint.insertDepletableNuclideKeys(c, blueprint)
             components[c.name] = c
             if spatialGrid:
                 componentLocators = gridDesign.getMultiLocator(

armi/reactor/blueprints/tests/test_blockBlueprints.py

@@ -339,20 +339,9 @@ def test_densityConsistentWithComponentConstructor(self):
             delta=biggerDelta,
         )
 
-        # This should be equal, but block construction calls applyHotHeightDensityReduction
-        # while programmatic construction allows components to exist in a state where
-        # their density is not consistent with material density.
-        self.assertNotAlmostEqual(
+        self.assertAlmostEqual(
             clad.getMassDensity(),
             programaticClad.getMassDensity(),
-            delta=biggerDelta,
-        )
-        # its off by a factor of thermal expansion
-        self.assertNotAlmostEqual(
-            clad.getMassDensity(),
-            programaticClad.getMassDensity()
-            * programaticClad.getThermalExpansionFactor(),
-            delta=biggerDelta,
         )
 
 

armi/reactor/components/component.py

@@ -329,13 +329,14 @@ def setProperties(self, properties):
 
     def applyMaterialMassFracsToNumberDensities(self):
         """
-        Set number densities for the component based on material mass fractions using hot temperatures.
+        Set the hot number densities for the component based on material mass fractions/density.
 
         Notes
         -----
-        - the density returned accounts for the radial expansion of the component
+        - the density returned accounts for the expansion of the component
           due to the difference in self.inputTemperatureInC and self.temperatureInC
-        - axial expansion effects are not included here.
+        - After the expansion, the density of the component should reflect the 3d
+          density of the material
 
         See Also
         --------
@@ -348,10 +349,12 @@ def applyMaterialMassFracsToNumberDensities(self):
         self.p.numberDensities = densityTools.getNDensFromMasses(
             density, self.material.p.massFrac
         )
+        self.applyHotHeightDensityReduction()
 
     def applyHotHeightDensityReduction(self):
         """
-        Adjust number densities to account for prescribed hot block heights (axial expansion).
+        Adjust number densities to account for hot block heights (axial expansion)
+        (crucial for preserving 3D density).
 
         Notes
         -----
@@ -363,6 +366,9 @@ def applyHotHeightDensityReduction(self):
         --------
         self.applyMaterialMassFracsToNumberDensities
         """
+        # this is the same as getThermalExpansionFactor but doesn't fail
+        # on non-fluid materials that have 0 or undefined thermal expansion
+        # (we don't want materials to fail on  __init__ which calls this)
         axialExpansionFactor = 1.0 + self.material.linearExpansionFactor(
             self.temperatureInC, self.inputTemperatureInC
         )

armi/reactor/converters/axialExpansionChanger.py

@@ -125,6 +125,23 @@ def setAssembly(self, a, setFuel=True):
         self.expansionData = ExpansionData(a, setFuel)
         self._isTopDummyBlockPresent()
 
+    def applyColdHeightMassIncrease(self):
+        """
+        Increase component mass because they are declared at cold dims
+
+        Notes
+        -----
+        A cold 1 cm tall component will have more mass that a component with the
+        same mass/length as a component with a hot height of 1 cm. This should be
+        called when the setting `inputHeightsConsideredHot` is used. This basically
+        undoes component.applyHotHeightDensityReduction
+        """
+        for c in self.linked.a.getComponents():
+            axialExpansionFactor = 1.0 + c.material.linearExpansionFactor(
+                c.temperatureInC, c.inputTemperatureInC
+            )
+            c.changeNDensByFactor(axialExpansionFactor)
+
     def _isTopDummyBlockPresent(self):
         """determines if top most block of assembly is a dummy block
 
@@ -175,9 +192,10 @@ def axiallyExpandAssembly(self, thermal: bool = False):
             # if ib == 0, leave block bottom = 0.0
             if ib > 0:
                 b.p.zbottom = self.linked.linkedBlocks[b][0].p.ztop
-            # if not in the dummy block, get expansion factor, do alignment, and modify block
-            if ib < (numOfBlocks - 1):
+            isDummyBlock = ib == (numOfBlocks - 1)
+            if not isDummyBlock:
                 for c in b:
+
                     growFrac = self.expansionData.getExpansionFactor(c)
                     runLog.debug(
                         msg="      Component {0}, growFrac = {1:.4e}".format(
@@ -222,8 +240,9 @@ def axiallyExpandAssembly(self, thermal: bool = False):
             # see also b.setHeight()
             # - the above not chosen due to call to calculateZCoords
             oldComponentVolumes = [c.getVolume() for c in b]
-            oldHeight = b.getHeight()
+            oldHeight = b.p.height
             b.p.height = b.p.ztop - b.p.zbottom
+
             _checkBlockHeight(b)
             _conserveComponentMass(b, oldHeight, oldComponentVolumes)
             # set block mid point and redo mesh

armi/reactor/converters/tests/test_axialExpansionChanger.py

@@ -758,20 +758,21 @@ class TestInputHeightsConsideredHot(unittest.TestCase):
     """verify thermal expansion for process loading of core"""
 
     def setUp(self):
-        """provide the base case"""
+        """This test uses a different armiRun.yaml than the default"""
+
         _o, r = loadTestReactor(
             os.path.join(TEST_ROOT, "detailedAxialExpansion"),
-            {"inputHeightsConsideredHot": True},
+            customSettings={"inputHeightsConsideredHot": True},
         )
         self.stdAssems = [a for a in r.core.getAssemblies()]
 
         _oCold, rCold = loadTestReactor(
             os.path.join(TEST_ROOT, "detailedAxialExpansion"),
-            {"inputHeightsConsideredHot": False},
+            customSettings={"inputHeightsConsideredHot": False},
         )
         self.testAssems = [a for a in rCold.core.getAssemblies()]
 
-    def test_coldAssemblyHeight(self):
+    def test_coldAssemblyExpansion(self):
         """block heights are cold and should be expanded
 
         Notes
@@ -801,6 +802,14 @@ def test_coldAssemblyHeight(self):
                     checkColdBlockHeight(bStd, bExp, self.assertEqual, "the same")
                 else:
                     checkColdBlockHeight(bStd, bExp, self.assertNotEqual, "different")
+                if bStd.hasFlags(Flags.FUEL):
+                    # fuel mass should grow because heights are considered cold heights
+                    # and a cold 1 cm column has more mass than a hot 1 cm column
+                    if not isinstance(
+                        bStd.getComponent(Flags.FUEL).material, custom.Custom
+                    ):
+                        # custom materials don't expand
+                        self.assertGreater(bExp.getMass("U235"), bStd.getMass("U235"))
 
 
 def checkColdBlockHeight(bStd, bExp, assertType, strForAssertion):

armi/reactor/converters/tests/test_blockConverter.py

@@ -55,7 +55,6 @@ def _perturbTemps(self, block, cName, tCold, tHot):
         "Give the component different ref and hot temperatures than in test_Blocks."
         c = block.getComponent(Flags.fromString(cName))
         c.refTemp, c.refHot = tCold, tHot
-        c.applyHotHeightDensityReduction()
         c.setTemperature(tHot)
         return block
 

armi/reactor/reactors.py

@@ -215,7 +215,6 @@ def __init__(self, name):
         self._circularRingPitch = 1.0
         self._automaticVariableMesh = False
         self._minMeshSizeRatio = 0.15
-        self._inputHeightsConsideredHot = True
         self._detailedAxialExpansion = False
 
     def setOptionsFromCs(self, cs):
@@ -227,7 +226,6 @@ def setOptionsFromCs(self, cs):
         self._circularRingPitch = cs["circularRingPitch"]
         self._automaticVariableMesh = cs["automaticVariableMesh"]
         self._minMeshSizeRatio = cs["minMeshSizeRatio"]
-        self._inputHeightsConsideredHot = cs["inputHeightsConsideredHot"]
         self._detailedAxialExpansion = cs["detailedAxialExpansion"]
 
     def __getstate__(self):
@@ -2338,6 +2336,12 @@ def _applyThermalExpansion(
         for a in assems:
             if not a.hasFlags(Flags.CONTROL):
                 axialExpChngr.setAssembly(a)
+                # this doesn't get applied to control assems, so CR will be interpreted
+                # as hot. This should be conservative because the control rods will
+                # be modeled as slightly shorter with the correct hot density. Density
+                # is more important than height, so we are forcing density to be correct
+                # since we can't do axial expansion (yet)
+                axialExpChngr.applyColdHeightMassIncrease()
                 axialExpChngr.expansionData.computeThermalExpansionFactors()
                 axialExpChngr.axiallyExpandAssembly(thermal=True)
         # resolve axially disjoint mesh (if needed)

armi/reactor/tests/test_blocks.py

@@ -2267,7 +2267,6 @@ def test_coldMass(self):
         and hot height.
         """
         fuel = self.b.getComponent(Flags.FUEL)
-        fuel.applyHotHeightDensityReduction()
         # set ref (input/cold) temperature.
         Thot = fuel.temperatureInC
         Tcold = fuel.inputTemperatureInC