Improving getOptimalAssemblyOrientation (#2019)

Improving  getOptimalAssemblyOrientation

The algorithm goes as follows:

1. Get all the pin powers and `IndexLocation`s from the block at the
   previous location
2. Obtain the `IndexLocation` of the pin with the highest burnup
3. For each possible rotation,
   - Find the new location with `HexGrid.rotateIndex`
   - Find the index where that location occurs in previous locations
   - Find the previous power at that location
4. Return the rotation with the lowest previous power

This algorithm assumes a few things:

1. `len(pinLocations) == len(pinPowers)` in both cases. This may make
   sense, but we've found some cases where this assumption breaks. Not
   even edge cases, like the C5G7 LWR benchmark.
2. Your assembly has at least 60 degree symmetry of fuel pins and
   powers. This means if we find a fuel pin with high burnup and rotate
   it 60 degrees, there should be another fuel pin at that lattice site.
   This is mostly a safe assumption since many hexagonal reactors have
   at least 60 degree symmetry of fuel pin layout. This assumption holds
   if you have a full hexagonal lattice of fuel pins as well.

---------

Co-authored-by: John Stilley <[email protected]>

Author: drewj-tp

armi/physics/fuelCycle/assemblyRotationAlgorithms.py

@@ -22,12 +22,18 @@
 .. warning:: Nothing should go in this file, but rotation algorithms.
 """
 import math
+from collections import defaultdict
 
 from armi import runLog
 from armi.physics.fuelCycle.hexAssemblyFuelMgmtUtils import (
     getOptimalAssemblyOrientation,
 )
 from armi.physics.fuelCycle.settings import CONF_ASSEM_ROTATION_STATIONARY
+from armi.physics.fuelCycle.utils import (
+    assemblyHasFuelPinBurnup,
+    assemblyHasFuelPinPowers,
+)
+from armi.reactor.assemblies import Assembly
 
 
 def _rotationNumberToRadians(rot: int) -> float:
@@ -49,42 +55,47 @@ def buReducingAssemblyRotation(fh):
     simpleAssemblyRotation : an alternative rotation algorithm
     """
     runLog.info("Algorithmically rotating assemblies to minimize burnup")
-    numRotated = 0
-    hist = fh.o.getInterface("history")
-    for aPrev in fh.moved:  # much more convenient to loop through aPrev first
+    # Store how we should rotate each assembly but don't perform the rotation just yet
+    # Consider assembly A is shuffled to a new location and rotated.
+    # Now, assembly B is shuffled to where assembly A used to be. We need to consider the
+    # power profile of A prior to it's rotation to understand the power profile B may see.
+    rotations: dict[int, list[Assembly]] = defaultdict(list)
+    for aPrev in fh.moved:
+        # If the assembly was out of the core, it will not have pin powers.
+        # No rotation information to be gained.
+        if aPrev.lastLocationLabel in Assembly.NOT_IN_CORE:
+            continue
         aNow = fh.r.core.getAssemblyWithStringLocation(aPrev.lastLocationLabel)
+        # An assembly in the SFP could have burnup but if it's coming from the load
+        # queue it's totally fresh. Skip a check over all pins in the model
+        if aNow.lastLocationLabel == Assembly.LOAD_QUEUE:
+            continue
         # no point in rotation if there's no pin detail
-        if aNow in hist.getDetailAssemblies():
-            _rotateByComparingLocations(aNow, aPrev)
-            numRotated += 1
+        if assemblyHasFuelPinPowers(aPrev) and assemblyHasFuelPinBurnup(aNow):
+            rot = getOptimalAssemblyOrientation(aNow, aPrev)
+            rotations[rot].append(aNow)
 
     if fh.cs[CONF_ASSEM_ROTATION_STATIONARY]:
-        for a in hist.getDetailAssemblies():
-            if a not in fh.moved:
-                _rotateByComparingLocations(a, a)
-                numRotated += 1
-
-    runLog.info("Rotated {0} assemblies".format(numRotated))
-
-
-def _rotateByComparingLocations(aNow, aPrev):
-    """Rotate an assembly based on its previous location.
-
-    Parameters
-    ----------
-    aNow : Assembly
-        Assembly to be rotated
-    aPrev : Assembly
-        Assembly that previously occupied the location of this assembly.
-        If ``aNow`` has not been moved, this should be ``aNow``
-    """
-    rot = getOptimalAssemblyOrientation(aNow, aPrev)
-    radians = _rotationNumberToRadians(rot)
-    aNow.rotate(radians)
-    (ring, pos) = aNow.spatialLocator.getRingPos()
-    runLog.important(
-        "Rotating Assembly ({0},{1}) to Orientation {2}".format(ring, pos, rot)
-    )
+        for a in filter(
+            lambda asm: asm not in fh.moved
+            and assemblyHasFuelPinPowers(asm)
+            and assemblyHasFuelPinBurnup(asm),
+            fh.r.core,
+        ):
+            rot = getOptimalAssemblyOrientation(a, a)
+            rotations[rot].append(a)
+
+    nRotations = 0
+    for rot, assems in filter(lambda item: item[0], rotations.items()):
+        # Radians used for the actual rotation. But a neater degrees print out is nice for logs
+        radians = _rotationNumberToRadians(rot)
+        degrees = round(math.degrees(radians), 3)
+        for a in assems:
+            runLog.important(f"Rotating assembly {a} {degrees} CCW.")
+            a.rotate(radians)
+            nRotations += 1
+
+    runLog.info(f"Rotated {nRotations} assemblies.")
 
 
 def simpleAssemblyRotation(fh):

armi/physics/fuelCycle/fuelHandlers.py

@@ -25,9 +25,9 @@
 This module also handles repeat shuffles when doing a restart.
 """
 # ruff: noqa: F401
+import inspect
 import os
 import re
-import warnings
 
 import numpy as np
 
@@ -114,10 +114,7 @@ def outage(self, factor=1.0):
             # The user can choose the algorithm method name directly in the settings
             if hasattr(rotAlgos, self.cs[CONF_ASSEMBLY_ROTATION_ALG]):
                 rotationMethod = getattr(rotAlgos, self.cs[CONF_ASSEMBLY_ROTATION_ALG])
-                try:
-                    rotationMethod()
-                except TypeError:
-                    rotationMethod(self)
+                rotationMethod(self)
             else:
                 raise RuntimeError(
                     "FuelHandler {0} does not have a rotation algorithm called {1}.\n"

armi/physics/fuelCycle/hexAssemblyFuelMgmtUtils.py

@@ -19,108 +19,113 @@
 -----
 We are keeping these in ARMI even if they appear unused internally.
 """
+import math
+import typing
+
 import numpy as np
 
 from armi import runLog
-from armi.reactor.flags import Flags
-from armi.utils.hexagon import getIndexOfRotatedCell
+from armi.physics.fuelCycle.utils import maxBurnupBlock, maxBurnupLocator
 from armi.utils.mathematics import findClosest
 
+if typing.TYPE_CHECKING:
+    from armi.reactor.assemblies import HexAssembly
+
 
-def getOptimalAssemblyOrientation(a, aPrev):
+def getOptimalAssemblyOrientation(a: "HexAssembly", aPrev: "HexAssembly") -> int:
     """
-    Get optimal assembly orientation/rotation to minimize peak burnup.
+    Get optimal hex assembly orientation/rotation to minimize peak burnup.
 
-    Notes
-    -----
-    Works by placing the highest-BU pin in the location (of 6 possible locations) with lowest
-    expected pin power. We evaluated "expected pin power" based on the power distribution in
-    aPrev, which is the previous assembly located here. If aPrev has no pin detail, then we must use its
-    corner fast fluxes to make an estimate.
+    .. impl:: Provide an algorithm for rotating hexagonal assemblies to equalize burnup
+        :id: I_ARMI_ROTATE_HEX_BURNUP
+        :implements: R_ARMI_ROTATE_HEX_BURNUP
 
     Parameters
     ----------
     a : Assembly object
         The assembly that is being rotated.
-
     aPrev : Assembly object
         The assembly that previously occupied this location (before the last shuffle).
-
-        If the assembly "a" was not shuffled, then "aPrev" = "a".
-
-        If "aPrev" has pin detail, then we will determine the orientation of "a" based on
-        the pin powers of "aPrev" when it was located here.
-
-        If "aPrev" does NOT have pin detail, then we will determine the orientation of "a" based on
-        the corner fast fluxes in "aPrev" when it was located here.
+        If the assembly "a" was not shuffled, it's sufficient to pass ``a``.
 
     Returns
     -------
-    rot : int
-        An integer from 0 to 5 representing the "orientation" of the assembly.
-        This orientation is relative to the current assembly orientation.
-        rot = 0 corresponds to no rotation.
-        rot represents the number of pi/3 counterclockwise rotations for the default orientation.
+    int
+        An integer from 0 to 5 representing the number of pi/3 (60 degree) counterclockwise
+        rotations from where ``a`` is currently oriented to the "optimal" orientation
 
-    Examples
-    --------
-    >>> getOptimalAssemblyOrientation(a, aPrev)
-    4
+    Raises
+    ------
+    ValueError
+        If there is insufficient information to determine the rotation of ``a``. This could
+        be due to a lack of fuel blocks or parameters like ``linPowByPin``.
 
-    See Also
-    --------
-    rotateAssemblies : calls this to figure out how to rotate
+    Notes
+    -----
+    Works by placing the highest-burnup pin in the location (of 6 possible locations) with lowest
+    expected pin power. We evaluated "expected pin power" based on the power distribution in
+    ``aPrev``, the previous assembly located where ``a`` is going. The algorithm goes as follows.
+
+    1. Get all the pin powers and ``IndexLocation``s from the block at the previous location/timenode.
+    2. Obtain the ``IndexLocation`` of the pin with the highest burnup in the current assembly.
+    3. For each possible rotation,
+        - Find the new location with ``HexGrid.rotateIndex``
+        - Find the index where that location occurs in previous locations
+        - Find the previous power at that location
+    4. Return the rotation with the lowest previous power
+
+    This algorithm assumes a few things.
+
+    1. ``len(HexBlock.getPinCoordinates()) == len(HexBlock.p.linPowByPin)`` and,
+       by extension, ``linPowByPin[i]`` is found at ``getPinCoordinates()[i]``.
+    2. Your assembly has at least 60 degree symmetry of fuel pins and
+       powers. This means if we find a fuel pin and rotate it 60 degrees, there should
+       be another fuel pin at that lattice site. This is mostly a safe assumption
+       since many hexagonal reactors have at least 60 degree symmetry of fuel pin layout.
+       This assumption holds if you have a full hexagonal lattice of fuel pins as well.
+    3. Fuel pins in ``a`` have similar locations in ``aPrev``. This is mostly a safe
+       assumption in that most fuel assemblies have similar layouts so it's plausible
+       that if ``a`` has a fuel pin at ``(1, 0, 0)``, so does ``aPrev``.
     """
-    # determine whether or not aPrev had pin details
-    fuelPrev = aPrev.getFirstBlock(Flags.FUEL)
-    if fuelPrev:
-        aPrevDetailFlag = fuelPrev.p.pinLocation[4] is not None
+    maxBuBlock = maxBurnupBlock(a)
+    if maxBuBlock.spatialGrid is None:
+        msg = f"Block {maxBuBlock} in {a} does not have a spatial grid. Cannot rotate."
+        runLog.error(msg)
+        raise ValueError(msg)
+    maxBuPinLocation = maxBurnupLocator(maxBuBlock)
+    # No need to rotate if max burnup pin is the center
+    if maxBuPinLocation.i == 0 and maxBuPinLocation.j == 0:
+        return 0
+
+    if aPrev is not a:
+        blockAtPreviousLocation = aPrev[a.index(maxBuBlock)]
     else:
-        aPrevDetailFlag = False
-
-    rot = 0  # default: no rotation
-    # First get pin index of maximum BU in this assembly.
-    _maxBuAssem, maxBuBlock = a.getMaxParam("percentBuMax", returnObj=True)
-    if maxBuBlock is None:
-        # no max block. They're all probably zero
-        return rot
-
-    # start at 0 instead of 1
-    maxBuPinIndexAssem = int(maxBuBlock.p.percentBuMaxPinLocation - 1)
-    bIndexMaxBu = a.index(maxBuBlock)
-
-    if maxBuPinIndexAssem == 0:
-        # Don't bother rotating if the highest-BU pin is the central pin. End this method.
-        return rot
-    else:
-        # transfer percentBuMax rotated pin index to non-rotated pin index
-        if aPrevDetailFlag:
-            # aPrev has pin detail
-            # Determine which of 6 possible rotated pin indices had the lowest power when aPrev was here.
-            prevAssemPowHereMIN = float("inf")
-
-            for possibleRotation in range(6):
-                index = getIndexOfRotatedCell(maxBuPinIndexAssem, possibleRotation)
-                # get pin power at this index in the previously assembly located here
-                # power previously at rotated index
-                prevAssemPowHere = aPrev[bIndexMaxBu].p.linPowByPin[index - 1]
-
-                if prevAssemPowHere is not None:
-                    runLog.debug(
-                        "Previous power in rotation {0} where pinLoc={1} is {2:.4E} W/cm"
-                        "".format(possibleRotation, index, prevAssemPowHere)
-                    )
-                    if prevAssemPowHere < prevAssemPowHereMIN:
-                        prevAssemPowHereMIN = prevAssemPowHere
-                        rot = possibleRotation
-        else:
-            raise ValueError(
-                "Cannot perform detailed rotation analysis without pin-level "
-                "flux information."
-            )
-
-        runLog.debug("Best relative rotation is {0}".format(rot))
-        return rot
+        blockAtPreviousLocation = maxBuBlock
+
+    previousLocations = blockAtPreviousLocation.getPinLocations()
+    previousPowers = blockAtPreviousLocation.p.linPowByPin
+    if len(previousLocations) != len(previousPowers):
+        msg = (
+            f"Inconsistent pin powers and number of pins in {blockAtPreviousLocation}. "
+            f"Found {len(previousLocations)} locations but {len(previousPowers)} powers."
+        )
+        runLog.error(msg)
+        raise ValueError(msg)
+
+    ringPowers = {
+        (loc.i, loc.j): p for loc, p in zip(previousLocations, previousPowers)
+    }
+
+    targetGrid = blockAtPreviousLocation.spatialGrid
+    candidateRotation = 0
+    candidatePower = ringPowers.get((maxBuPinLocation.i, maxBuPinLocation.j), math.inf)
+    for rot in range(1, 6):
+        candidateLocation = targetGrid.rotateIndex(maxBuPinLocation, rot)
+        newPower = ringPowers.get((candidateLocation.i, candidateLocation.j), math.inf)
+        if newPower < candidatePower:
+            candidateRotation = rot
+            candidatePower = newPower
+    return candidateRotation
 
 
 def buildRingSchedule(