Accounting for symmetry in volume-integrated parameters (#2017)

armi/bookkeeping/tests/test_historyTracker.py

@@ -184,7 +184,6 @@ def test_historyParameters(self):
             params[param] = []
             for ts, years in enumerate(timesInYears):
                 cycle, node = utils.getCycleNodeFromCumulativeNode(ts, self.o.cs)
-
                 params[param].append(
                     hti.getBlockHistoryVal(bName, param, (cycle, node))
                 )
@@ -196,7 +195,8 @@ def test_historyParameters(self):
         # verify the power parameter is retrievable from the history
         self.assertEqual(o.cs["power"], 1000000000.0)
         self.assertAlmostEqual(params["power"][0], 360, delta=0.1)
-        self.assertEqual(params["power"][0], params["power"][1])
+        # assembly was moved to the central location with 1/3rd symmetry
+        self.assertEqual(params["power"][0] / 3, params["power"][1])
 
         # verify the power density parameter is retrievable from the history
         self.assertAlmostEqual(params["pdens"][0], 0.0785, delta=0.001)

armi/physics/fuelCycle/fuelHandlers.py

@@ -38,6 +38,7 @@
 from armi.physics.fuelCycle.settings import CONF_ASSEMBLY_ROTATION_ALG
 from armi.reactor.flags import Flags
 from armi.utils.customExceptions import InputError
+from armi.reactor.parameters import ParamLocation
 
 
 class FuelHandler:
@@ -212,10 +213,18 @@ def _compareAssem(candidate, current):
     @staticmethod
     def _getParamMax(a, paramName, blockLevelMax=True):
         """Get parameter with Block-level maximum."""
+        multiplier = a.getSymmetryFactor()
+        if multiplier != 1:
+            # handle special case: volume-integrated parameters where symmetry factor is not 1
+            isVolumeIntegrated = (
+                a.getBlocks()[0].p.paramDefs[paramName].location
+                == ParamLocation.VOLUME_INTEGRATED
+            )
+            multiplier = a.getSymmetryFactor() if isVolumeIntegrated else 1.0
         if blockLevelMax:
-            return a.getChildParamValues(paramName).max()
+            return a.getChildParamValues(paramName).max() * multiplier
 
-        return a.p[paramName]
+        return a.p[paramName] * multiplier
 
     def findAssembly(
         self,

armi/physics/fuelCycle/tests/test_fuelHandlers.py

@@ -261,8 +261,9 @@ def test_width(self):
         for ring, power in zip(range(1, 8), range(10, 80, 10)):
             aList = assemsByRing[ring]
             for a in aList:
+                sf = a.getSymmetryFactor()  # center assembly is only 1/3rd in the core
                 for b in a:
-                    b.p.power = power
+                    b.p.power = power / sf
 
         paramName = "power"
         # 1 ring outer and inner from ring 3

armi/reactor/assemblies.py

@@ -22,6 +22,7 @@
 import pickle
 from random import randint
 from typing import ClassVar, Optional, Type
+from collections.abc import Iterable
 
 import numpy as np
 from scipy import interpolate
@@ -206,13 +207,34 @@ def insert(self, index, obj):
 
     def moveTo(self, locator):
         """Move an assembly somewhere else."""
+        oldSymmetryFactor = self.getSymmetryFactor()
         composites.Composite.moveTo(self, locator)
         if self.lastLocationLabel != self.DATABASE:
             self.p.numMoves += 1
             self.p.daysSinceLastMove = 0.0
         self.parent.childrenByLocator[locator] = self
         # symmetry may have changed (either moving on or off of symmetry line)
         self.clearCache()
+        self.scaleParamsToNewSymmetryFactor(oldSymmetryFactor)
+
+    def scaleParamsToNewSymmetryFactor(self, oldSymmetryFactor):
+        scalingFactor = oldSymmetryFactor / self.getSymmetryFactor()
+        if scalingFactor == 1:
+            return
+
+        volIntegratedParamsToScale = self.getBlocks()[0].p.paramDefs.atLocation(
+            ParamLocation.VOLUME_INTEGRATED
+        )
+        for b in self.getBlocks():
+            for param in volIntegratedParamsToScale:
+                name = param.name
+                if b.p[name] is None or isinstance(b.p[name], str):
+                    continue
+                elif isinstance(b.p[name], Iterable):
+                    b.p[name] = [value * scalingFactor for value in b.p[name]]
+                else:
+                    # numpy array or other
+                    b.p[name] = b.p[name] * scalingFactor
 
     def getNum(self):
         """Return unique integer for this assembly."""

armi/reactor/blockParameters.py

@@ -81,24 +81,6 @@ def getBlockParameterDefinitions():
             description="Ratio of fissile mass to heavy metal mass at block-level",
         )
 
-        pb.defParam(
-            "molesHmBOL",
-            units=f"{units.MOLES}",
-            description="Total number of atoms of heavy metal at BOL assuming a full assembly",
-        )
-
-        pb.defParam(
-            "massHmBOL",
-            units=units.GRAMS,
-            description="Mass of heavy metal at BOL",
-        )
-
-        pb.defParam(
-            "initialB10ComponentVol",
-            units=f"{units.CM}^3",
-            description="cc's of un-irradiated, cold B10 containing component (includes full volume if any B10)",
-        )
-
         pb.defParam(
             "molesHmBOLByPin",
             units=f"{units.MOLES}",
@@ -108,12 +90,6 @@ def getBlockParameterDefinitions():
             location=ParamLocation.CHILDREN,
         )
 
-        pb.defParam(
-            "molesHmNow",
-            units=f"{units.MOLES}",
-            description="Total number of atoms of heavy metal",
-        )
-
         pb.defParam(
             "newDPA",
             units=units.DPA,
@@ -161,6 +137,37 @@ def getBlockParameterDefinitions():
             categories=["cumulative"],
         )
 
+    with pDefs.createBuilder(
+        default=0.0, location=ParamLocation.VOLUME_INTEGRATED, categories=["depletion"]
+    ) as pb:
+
+        pb.defParam(
+            "molesHmNow",
+            units=f"{units.MOLES}",
+            description="Total number of atoms of heavy metal",
+        )
+
+        pb.defParam(
+            "molesHmBOL",
+            units=f"{units.MOLES}",
+            description="Total number of atoms of heavy metal at BOL.",
+        )
+
+        pb.defParam(
+            "massHmBOL",
+            units=units.GRAMS,
+            description="Mass of heavy metal at BOL",
+        )
+
+        pb.defParam(
+            "initialB10ComponentVol",
+            units=f"{units.CM}^3",
+            description=(
+                "cc's of un-irradiated, cold B10 containing component "
+                "(includes full volume of any components with B10)"
+            ),
+        )
+
     pDefs.add(
         Parameter(
             name="depletionMatrix",

armi/reactor/blocks.py

@@ -790,21 +790,16 @@ def completeInitialLoading(self, bolBlock=None):
 
         self.p.enrichmentBOL = self.getFissileMassEnrich()
         massHmBOL = 0.0
-        sf = self.getSymmetryFactor()
         for child in self:
-            # multiplying by sf ends up cancelling out the symmetry factor used in
-            # Component.getMass(). So massHmBOL does not respect the symmetry factor.
-            hmMass = child.getHMMass() * sf
+            hmMass = child.getHMMass()
             massHmBOL += hmMass
             # Components have the following parameters but not every composite will
             # massHmBOL, molesHmBOL, puFrac
             if isinstance(child, components.Component):
                 child.p.massHmBOL = hmMass
-                # to stay consistent with massHmBOL, molesHmBOL and puFrac should be
-                # independent of sf. As such, the need to be scaled by 1/sf.
-                child.p.molesHmBOL = child.getHMMoles() / sf
+                child.p.molesHmBOL = child.getHMMoles()
                 child.p.puFrac = (
-                    self.getPuMoles() / sf / child.p.molesHmBOL
+                    self.getPuMoles() / child.p.molesHmBOL
                     if child.p.molesHmBOL > 0.0
                     else 0.0
                 )

armi/reactor/composites.py

@@ -1947,29 +1947,18 @@ def getHMMass(self):
 
     def getHMMoles(self):
         """
-        Get the number of moles of heavy metal in this object in full symmetry.
+        Get the number of moles of heavy metal in this object.
 
         Notes
         -----
-        If an object is on a symmetry line, the number of moles will be scaled up by the
-        symmetry factor. This is done because this is typically used for tracking
-        burnup, and BOL moles are computed in full objects too so there are no
-        complications as things move on and off of symmetry lines.
-
-        Warning
-        -------
-        getHMMoles is different than every other get mass call since it multiplies by
-        symmetry factor but getVolume() on the block level divides by symmetry factor
-        causing them to cancel out.
-
-        This was needed so that HM moles mass did not change based on if the
-        block/assembly was on a symmetry line or not.
+        If an object is on a symmetry line, the volume reported by getVolume
+        is reduced to reflect that the block is not wholly within the reactor. This
+        reduction in volume reduces the reported HM moles.
         """
         return (
             self.getHMDens()
             / units.MOLES_PER_CC_TO_ATOMS_PER_BARN_CM
             * self.getVolume()
-            * self.getSymmetryFactor()
         )
 
     def getHMDens(self):
@@ -3129,12 +3118,7 @@ def getPuMoles(self):
         nucNames = [nuc.name for nuc in elements.byZ[94].nuclides]
         puN = sum(self.getNuclideNumberDensities(nucNames))
 
-        return (
-            puN
-            / units.MOLES_PER_CC_TO_ATOMS_PER_BARN_CM
-            * self.getVolume()
-            * self.getSymmetryFactor()
-        )
+        return puN / units.MOLES_PER_CC_TO_ATOMS_PER_BARN_CM * self.getVolume()
 
 
 class StateRetainer:

armi/reactor/tests/test_assemblies.py

@@ -204,7 +204,6 @@ def setUp(self):
         )
 
         self.assembly = makeTestAssembly(NUM_BLOCKS, self.assemNum, r=self.r)
-        self.r.core.add(self.assembly)
 
         # Use these if they are needed
         self.blockParams = {
@@ -267,6 +266,7 @@ def setUp(self):
             self.assembly.add(b)
             self.blockList.append(b)
 
+        self.r.core.add(self.assembly)
         self.assembly.calculateZCoords()
 
     def test_isOnWhichSymmetryLine(self):
@@ -345,6 +345,35 @@ def test_moveTo(self):
         cur = self.assembly.spatialLocator
         self.assertEqual(cur, ref)
 
+    def test_scaleParamsWhenMoved(self):
+        """Volume integrated parameters must be scaled when an assembly is placed on a core boundary."""
+        blockParams = {
+            # volume integrated parameters
+            "massHmBOL": 9.0,
+            "molesHmBOL": np.array([[1, 2, 3], [4, 5, 6]]),  # ndarray for testing
+            "adjMgFlux": [1, 2, 3],  # Should normally be an ndarray, list for testing
+            "lastMgFlux": "foo",  # Should normally be an ndarray, str for testing
+        }
+        for b in self.assembly.getBlocks(Flags.FUEL):
+            b.p.update(blockParams)
+
+        i, j = grids.HexGrid.getIndicesFromRingAndPos(1, 1)
+        locator = self.r.core.spatialGrid[i, j, 0]
+        self.assertEqual(self.assembly.getSymmetryFactor(), 1)
+        self.assembly.moveTo(locator)
+        self.assertEqual(self.assembly.getSymmetryFactor(), 3)
+        for b in self.assembly.getBlocks(Flags.FUEL):
+            # float
+            assert_allclose(b.p["massHmBOL"] / blockParams["massHmBOL"], 1 / 3)
+            # np.ndarray
+            assert_allclose(b.p["molesHmBOL"] / blockParams["molesHmBOL"], 1 / 3)
+            # list
+            assert_allclose(
+                np.array(b.p["adjMgFlux"]) / np.array(blockParams["adjMgFlux"]), 1 / 3
+            )
+            # string
+            self.assertEqual(b.p["lastMgFlux"], blockParams["lastMgFlux"])
+
     def test_getName(self):
         cur = self.assembly.getName()
         ref = self.name

armi/reactor/tests/test_blocks.py

@@ -1230,34 +1230,24 @@ def test_completeInitialLoading(self, mock_sf):
 
         sf = self.block.getSymmetryFactor()
         cur = self.block.p.molesHmBOL
-        ref = (
-            self.block.getHMDens()
-            / MOLES_PER_CC_TO_ATOMS_PER_BARN_CM
-            * height
-            * area
-            * sf
-        )
+        ref = self.block.getHMDens() / MOLES_PER_CC_TO_ATOMS_PER_BARN_CM * height * area
         self.assertAlmostEqual(cur, ref, places=12)
 
         totalHMMass = 0.0
         for c in self.block:
             nucs = c.getNuclides()
             hmNucs = [nuc for nuc in nucs if nucDir.isHeavyMetal(nuc)]
             hmNDens = {hmNuc: c.getNumberDensity(hmNuc) for hmNuc in hmNucs}
-            hmMass = densityTools.calculateMassDensity(hmNDens) * c.getVolume()
+            # use sf to account for only a 1/sf portion of the component being in the block
+            hmMass = densityTools.calculateMassDensity(hmNDens) * c.getVolume() / sf
             totalHMMass += hmMass
             if hmMass:
-                # hmMass does not need to account for sf since what's calculated in blocks.completeInitialLoading
-                # ends up cancelling out sf
                 self.assertAlmostEqual(c.p.massHmBOL, hmMass, places=12)
-                # since sf is cancelled out in massHmBOL, there needs to be a factor 1/sf here to cancel out the
-                # factor of sf in getHMMoles.
                 self.assertAlmostEqual(
                     c.p.molesHmBOL,
                     sum(ndens for ndens in hmNDens.values())
                     / units.MOLES_PER_CC_TO_ATOMS_PER_BARN_CM
-                    * c.getVolume()
-                    / sf,
+                    * c.getVolume(),
                     places=12,
                 )
             else:

armi/tests/tutorials/data_model.ipynb

@@ -294,14 +294,18 @@
     "mgFluxBase = np.arange(5)\n",
     "def setFakePower(core):\n",
     "    for a in core:\n",
+    "        sf = a.getSymmetryFactor()\n",
     "        for b in a:\n",
     "            vol = b.getVolume()\n",
     "            coords = b.spatialLocator.getGlobalCoordinates()\n",
     "            r = np.linalg.norm(abs(coords-center))\n",
     "            fuelFlag = 10 if b.isFuel() else 1.0\n",
-    "            b.p.power = peakPower / r**2 * fuelFlag\n",
+    "            # Use the symmetry factor to account for the central assembly being split\n",
+    "            b.p.power = peakPower / r**2 * fuelFlag / sf\n",
     "            b.p.pdens = b.p.power/vol\n",
     "            b.p.mgFlux = mgFluxBase*b.p.pdens\n",
+    "            if b.isFuel():\n",
+    "                print(b.p.power, b.getLocation())\n",
     "setFakePower(core)"
    ]
   },