Conserve mass at component level for a.setBlockMesh() (#1665)

* Conserve mass at component level for a.setBlockMesh()

a.setBlockMesh() had an option, `conserveMassFlag="auto"`, which would
try to automatically conserve mass of the fuel while allowing other
components to change mass with changing block heights. This conservation
was based upon getting the nuclides associated with the fuel component,
and then scaling the density of those nuclides based on the change in
block height using b.setNumberDensities().

This method does not work when the fuel has nuclides that are also in
other components of the block, like cladding or structure. In this case,
the mass of those components will also be conserved, which is not the
desired behavior.

Instead, a.setBlockMesh() needs to conserve mass by setting number
densities directly on the components where mass needs to be conserved,
instead of setting the densities on the full block.

armi/reactor/assemblies.py

@@ -32,6 +32,7 @@
 from armi.reactor import composites
 from armi.reactor import grids
 from armi.reactor.flags import Flags
+from armi.materials.material import Fluid
 from armi.reactor.parameters import ParamLocation
 
 
@@ -673,7 +674,7 @@ def makeAxialSnapList(self, refAssem=None, refMesh=None, force=False):
 
     def _shouldMassBeConserved(self, belowFuelColumn, b):
         """
-        Determine from a rule set if the mass of a block should be conserved during axial expansion.
+        Determine from a rule set if the mass of a block component should be conserved during axial expansion.
 
         Parameters
         ----------
@@ -689,8 +690,8 @@ def _shouldMassBeConserved(self, belowFuelColumn, b):
         conserveMass : boolean
             Should the mass be conserved in this block
 
-        adjustList : list of nuclides
-            What nuclides should have their mass conserved (if any)
+        conserveComponents : list of components
+            What components should have their mass conserved (if any)
 
         belowFuelColumn : boolean
             Update whether the block is above or below a fuel column
@@ -703,32 +704,31 @@ def _shouldMassBeConserved(self, belowFuelColumn, b):
         if b.hasFlags(Flags.FUEL):
             # fuel block
             conserveMass = True
-            adjustList = b.getComponent(Flags.FUEL).getNuclides()
+            conserveComponents = b.getComponents(Flags.FUEL)
         elif self.hasFlags(Flags.FUEL):
             # non-fuel block of a fuel assembly.
             if belowFuelColumn:
                 # conserve mass of everything below the fuel so as to not invalidate
                 # grid-plate dose calcs.
                 conserveMass = True
-                adjustList = b.getNuclides()
-                # conserve mass of everything except coolant.
-                coolant = b.getComponent(Flags.COOLANT)
-                coolantList = coolant.getNuclides() if coolant else []
-                for nuc in coolantList:
-                    if nuc in adjustList:
-                        adjustList.remove(nuc)
+                # conserve mass of everything except fluids.
+                conserveComponents = [
+                    comp
+                    for comp in b.getComponents()
+                    if not isinstance(comp.material, Fluid)
+                ]
             else:
                 # plenum or above block in fuel assembly. don't conserve mass.
                 conserveMass = False
-                adjustList = None
+                conserveComponents = []
         else:
             # non fuel block in non-fuel assem. Don't conserve mass.
             conserveMass = False
-            adjustList = None
+            conserveComponents = []
 
-        return conserveMass, adjustList
+        return conserveMass, conserveComponents
 
-    def setBlockMesh(self, blockMesh, conserveMassFlag=False, adjustList=None):
+    def setBlockMesh(self, blockMesh, conserveMassFlag=False):
         """
         Snaps the axial mesh points of this assembly to correspond with the reference mesh.
 
@@ -756,7 +756,14 @@ def setBlockMesh(self, blockMesh, conserveMassFlag=False, adjustList=None):
         Parameters
         ----------
         blockMesh : iterable
-            a list of floats describing the upper mesh points of each block in cm.
+            A list of floats describing the upper mesh points of each block in cm.
+
+        conserveMassFlag : bool or str
+            Option for how to treat mass conservation when the block mesh changes.
+            Conservation of mass for fuel components is enabled by
+            conserveMassFlag="auto". If not auto, a boolean value should be
+            passed. The default is False, which does not conserve any masses.
+            True conserves mass for all components.
 
         See Also
         --------
@@ -796,15 +803,19 @@ def setBlockMesh(self, blockMesh, conserveMassFlag=False, adjustList=None):
                 return
 
             if conserveMassFlag == "auto":
-                conserveMass, adjustList = self._shouldMassBeConserved(
+                conserveMass, conserveComponents = self._shouldMassBeConserved(
                     belowFuelColumn, b
                 )
             else:
                 conserveMass = conserveMassFlag
-
-            b.setHeight(
-                newTop - zBottom, conserveMass=conserveMass, adjustList=adjustList
-            )
+                conserveComponents = b.getComponents()
+
+            oldBlockHeight = b.getHeight()
+            b.setHeight(newTop - zBottom)
+            if conserveMass:
+                heightRatio = oldBlockHeight / b.getHeight()
+                for c in conserveComponents:
+                    c.changeNDensByFactor(heightRatio)
             zBottom = newTop
 
         self.calculateZCoords()

armi/reactor/tests/test_assemblies.py

@@ -1163,6 +1163,12 @@ def test_snapAxialMeshToReferenceConservingMassBasedOnBlockIgniter(self):
 
         igniterFuel = self.r.core.childrenByLocator[grid[0, 0, 0]]
         # gridplate, fuel, fuel, fuel, plenum
+        for b in igniterFuel.getBlocks(Flags.FUEL):
+            fuelComp = b.getComponent(Flags.FUEL)
+            # add isotopes from clad and coolant to fuel component to test mass conservation
+            # mass should only be conserved within fuel component, not over the whole block
+            fuelComp.setNumberDensity("FE56", 1e-10)
+            fuelComp.setNumberDensity("NA23", 1e-10)
         b = igniterFuel[0]
         coolantNucs = b.getComponent(Flags.COOLANT).getNuclides()
         coolMass = 0
@@ -1175,6 +1181,8 @@ def test_snapAxialMeshToReferenceConservingMassBasedOnBlockIgniter(self):
         igniterHMMass1 = b.getHMMass()
         igniterZircMass1 = b.getMass("ZR")
         igniterFuelBlockMass = b.getMass()
+        igniterDuctMass = b.getComponent(Flags.DUCT).getMass()
+        igniterCoolMass = b.getComponent(Flags.COOLANT).getMass()
 
         coolMass = 0
         b = igniterFuel[4]
@@ -1199,6 +1207,8 @@ def test_snapAxialMeshToReferenceConservingMassBasedOnBlockIgniter(self):
         b = igniterFuel[1]
         igniterHMMass1AfterExpand = b.getHMMass()
         igniterZircMass1AfterExpand = b.getMass("ZR")
+        igniterDuctMassAfterExpand = b.getComponent(Flags.DUCT).getMass()
+        igniterCoolMassAfterExpand = b.getComponent(Flags.COOLANT).getMass()
 
         coolMass = 0
         b = igniterFuel[4]
@@ -1209,6 +1219,14 @@ def test_snapAxialMeshToReferenceConservingMassBasedOnBlockIgniter(self):
         self.assertAlmostEqual(igniterMassGrid, igniterMassGridAfterExpand, 7)
         self.assertAlmostEqual(igniterHMMass1, igniterHMMass1AfterExpand, 7)
         self.assertAlmostEqual(igniterZircMass1, igniterZircMass1AfterExpand, 7)
+        # demonstrate that the duct and coolant mass are not conserved.
+        # number density stays constant, mass is scaled by ratio of new to old height
+        self.assertAlmostEqual(
+            igniterDuctMass, igniterDuctMassAfterExpand * 25.0 / 26.0, 7
+        )
+        self.assertAlmostEqual(
+            igniterCoolMass, igniterCoolMassAfterExpand * 25.0 / 26.0, 7
+        )
         # Note the masses are linearly different by the amount that the plenum shrunk
         self.assertAlmostEqual(
             igniterPlenumMass, igniterPlenumMassAfterExpand * 75 / 67.0, 7
@@ -1233,6 +1251,8 @@ def test_snapAxialMeshToReferenceConservingMassBasedOnBlockIgniter(self):
         igniterHMMass1AfterShrink = b.getHMMass()
         igniterZircMass1AfterShrink = b.getMass("ZR")
         igniterFuelBlockMassAfterShrink = b.getMass()
+        igniterDuctMassAfterShrink = b.getComponent(Flags.DUCT).getMass()
+        igniterCoolMassAfterShrink = b.getComponent(Flags.COOLANT).getMass()
 
         coolMass = 0
         b = igniterFuel[4]
@@ -1245,6 +1265,8 @@ def test_snapAxialMeshToReferenceConservingMassBasedOnBlockIgniter(self):
         self.assertAlmostEqual(igniterHMMass1, igniterHMMass1AfterShrink, 7)
         self.assertAlmostEqual(igniterZircMass1, igniterZircMass1AfterShrink, 7)
         self.assertAlmostEqual(igniterFuelBlockMass, igniterFuelBlockMassAfterShrink, 7)
+        self.assertAlmostEqual(igniterDuctMass, igniterDuctMassAfterShrink, 7)
+        self.assertAlmostEqual(igniterCoolMass, igniterCoolMassAfterShrink, 7)
         self.assertAlmostEqual(igniterPlenumMass, igniterPlenumMassAfterShrink, 7)
 
     def test_snapAxialMeshToReferenceConservingMassBasedOnBlockShield(self):

doc/release/0.3.rst

@@ -8,6 +8,7 @@ Release Date: TBD
 
 New Features
 ------------
+#. Conserve mass by component in assembly.setBlockMesh(). (`PR#1665 <https://github.com/terrapower/armi/pull/1665>`_)
 #. TBD
 
 API Changes