Improving getOptimalAssemblyOrientation

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(