diff --git a/armi/physics/fuelCycle/assemblyRotationAlgorithms.py b/armi/physics/fuelCycle/assemblyRotationAlgorithms.py index 30c87ba2e..0230f2fbc 100644 --- a/armi/physics/fuelCycle/assemblyRotationAlgorithms.py +++ b/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): diff --git a/armi/physics/fuelCycle/fuelHandlers.py b/armi/physics/fuelCycle/fuelHandlers.py index 0c1f8c452..abe079d0b 100644 --- a/armi/physics/fuelCycle/fuelHandlers.py +++ b/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" diff --git a/armi/physics/fuelCycle/hexAssemblyFuelMgmtUtils.py b/armi/physics/fuelCycle/hexAssemblyFuelMgmtUtils.py index 9afe2d4c8..c39581548 100644 --- a/armi/physics/fuelCycle/hexAssemblyFuelMgmtUtils.py +++ b/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( diff --git a/armi/physics/fuelCycle/tests/test_assemblyRotationAlgorithms.py b/armi/physics/fuelCycle/tests/test_assemblyRotationAlgorithms.py index bbbdec230..66ea0a6ed 100644 --- a/armi/physics/fuelCycle/tests/test_assemblyRotationAlgorithms.py +++ b/armi/physics/fuelCycle/tests/test_assemblyRotationAlgorithms.py @@ -19,36 +19,315 @@ These algorithms are defined in assemblyRotationAlgorithms.py, but they are used in: ``FuelHandler.outage()``. """ +import copy +import enum +import math +import typing +from unittest import TestCase, mock + +import numpy as np + from armi.physics.fuelCycle import assemblyRotationAlgorithms as rotAlgos from armi.physics.fuelCycle import fuelHandlers +from armi.physics.fuelCycle.hexAssemblyFuelMgmtUtils import ( + getOptimalAssemblyOrientation, +) from armi.physics.fuelCycle.settings import CONF_ASSEM_ROTATION_STATIONARY from armi.physics.fuelCycle.tests.test_fuelHandlers import addSomeDetailAssemblies -from armi.physics.fuelCycle.tests.test_fuelHandlers import FuelHandlerTestHelper +from armi.reactor.assemblies import HexAssembly +from armi.reactor.blocks import HexBlock from armi.reactor.flags import Flags +from armi.reactor.tests import test_reactors -class TestFuelHandlerMgmtTools(FuelHandlerTestHelper): - def test_buReducingAssemblyRotation(self): - fh = fuelHandlers.FuelHandler(self.o) - hist = self.o.getInterface("history") - newSettings = {CONF_ASSEM_ROTATION_STATIONARY: True} +class MockFuelHandler(fuelHandlers.FuelHandler): + """Implements the entire interface but with empty methods.""" + + def chooseSwaps(self, *args, **kwargs): + pass + + +class _PinLocations(enum.IntEnum): + """Zero-indexed locations for specific points of interest. + + If a data vector has an entry to all ``self.N_PINS=169`` pins in the test model, + then ``data[PIN_LOCATIONS.UPPER_RIGHT_VERTEX]`` will access the data for the pin + along the upper right 60 symmetry line. Since we're dealing with rotations here, it + does not need to literally be the pin at the vertex. Just along the symmetry line + to help explain tests. + + The use case here is setting the pin or burnup array to be a constant value, but + using a single max or minimum value to determine rotation. + """ + + CENTER = 0 + UPPER_RIGHT_VERTEX = 1 + UPPER_LEFT_VERTEX = 2 + DUE_LEFT_VERTEX = 3 + LOWER_LEFT_VERTEX = 4 + LOWER_RIGHT_VERTEX = 5 + DUE_RIGHT_VERTEX = 6 + + +class ShuffleAndRotateTestHelper(TestCase): + """Fixture class to assist in testing rotation of assemblies via the fuel handler.""" + + N_PINS = 169 + + def setUp(self): + self.o, self.r = test_reactors.loadTestReactor() + self.r.core.locateAllAssemblies() + + @staticmethod + def ensureBlockHasSpatialGrid(b: HexBlock): + """If ``b`` does not have a spatial grid, auto create one.""" + if b.spatialGrid is None: + b.getPinPitch = mock.Mock(return_value=1.1) + b.autoCreateSpatialGrids() + + def setAssemblyPinBurnups(self, a: HexAssembly, burnups: np.ndarray): + """Prepare the assembly that will be shuffled and rotated.""" + peakBu = burnups.max() + for b in a.getChildrenWithFlags(Flags.FUEL): + self.ensureBlockHasSpatialGrid(b) + b.p.percentBuPeak = peakBu + for c in b.getChildrenWithFlags(Flags.FUEL): + c.p.pinPercentBu = burnups + + def setAssemblyPinPowers(self, a: HexAssembly, pinPowers: np.ndarray): + """Prep the assembly that existed at the site a shuffled assembly will occupy.""" + for b in a.getChildrenWithFlags(Flags.FUEL): + self.ensureBlockHasSpatialGrid(b) + b.p.linPowByPin = pinPowers + + def powerWithMinValue(self, minIndex: int) -> np.ndarray: + """Create a vector of pin powers with a minimum value at a given index.""" + data = np.ones(self.N_PINS) + data[minIndex] = 0 + return data + + def burnupWithMaxValue(self, maxIndex: int) -> np.ndarray: + """Create a vector of pin burnups with a maximum value at a given index.""" + data = np.zeros(self.N_PINS) + data[maxIndex] = 50 + return data + + def compareMockedToExpectedRotation( + self, nRotations: int, mRotate: mock.Mock, msg: typing.Optional[str] = None + ): + """Helper function to check the mocked rotate and compare against expected rotation.""" + expectedRadians = nRotations * math.pi / 3 + (actualRadians,) = mRotate.call_args.args + self.assertAlmostEqual(actualRadians, expectedRadians, msg=msg) + + +class TestOptimalAssemblyRotation(ShuffleAndRotateTestHelper): + """Test the burnup dependent assembly rotation methods.""" + + def setUp(self): + super().setUp() + self.assembly: HexAssembly = self.r.core.getFirstAssembly(Flags.FUEL) + + def test_flatPowerNoRotation(self): + """If all pin powers are identical, no rotation is suggested.""" + burnups = self.burnupWithMaxValue(_PinLocations.UPPER_LEFT_VERTEX) + powers = np.ones_like(burnups) + self.setAssemblyPinBurnups(self.assembly, burnups) + self.setAssemblyPinPowers(self.assembly, powers) + rot = getOptimalAssemblyOrientation(self.assembly, self.assembly) + self.assertEqual(rot, 0) + + def test_maxBurnupAtCenterNoRotation(self): + """If max burnup pin is at the center, no rotation is suggested.""" + burnups = self.burnupWithMaxValue(_PinLocations.CENTER) + powers = np.zeros_like(burnups) + self.setAssemblyPinBurnups(self.assembly, burnups) + self.setAssemblyPinPowers(self.assembly, powers) + rot = getOptimalAssemblyOrientation(self.assembly, self.assembly) + self.assertEqual(rot, 0) + + def test_oppositeRotation(self): + """Test a 180 degree rotation is suggested when the max burnup pin is opposite the lowest power pin. + + Use the second ring of the hexagon because it's easier to write out pin locations + and check work. + + .. test:: Test the burnup equalizing rotation algorithm. + :id: T_ARMI_ROTATE_HEX_BURNUP + :tests: R_ARMI_ROTATE_HEX_BURNUP + :acceptance_criteria: After rotating a hexagonal assembly, confirm the pin with the highest burnup is + in the same sector as pin with the lowest power in the high burnup pin's ring. + + Notes + ----- + Use zero-indexed pin location not pin ID to assign burnups and powers. Since + we have a single component, ``Block.p.linPowByPin[i] <-> Component.p.pinPercentBu[i]`` + """ + shuffledAssembly = self.assembly + previousAssembly = copy.deepcopy(shuffledAssembly) + pairs = ( + (_PinLocations.DUE_RIGHT_VERTEX, _PinLocations.DUE_LEFT_VERTEX), + (_PinLocations.UPPER_LEFT_VERTEX, _PinLocations.LOWER_RIGHT_VERTEX), + (_PinLocations.UPPER_RIGHT_VERTEX, _PinLocations.LOWER_LEFT_VERTEX), + (_PinLocations.DUE_LEFT_VERTEX, _PinLocations.DUE_RIGHT_VERTEX), + (_PinLocations.LOWER_RIGHT_VERTEX, _PinLocations.UPPER_LEFT_VERTEX), + (_PinLocations.LOWER_LEFT_VERTEX, _PinLocations.UPPER_RIGHT_VERTEX), + ) + for startPin, oppositePin in pairs: + powers = self.powerWithMinValue(oppositePin) + burnups = self.burnupWithMaxValue(startPin) + self.setAssemblyPinBurnups(shuffledAssembly, burnups) + self.setAssemblyPinPowers(previousAssembly, powers) + rot = getOptimalAssemblyOrientation(shuffledAssembly, previousAssembly) + # 180 degrees is three 60 degree rotations + self.assertEqual(rot, 3, msg=f"{startPin=} :: {oppositePin=}") + + def test_noBlocksWithBurnup(self): + """Require at least one block to have burnup.""" + with self.assertRaisesRegex(ValueError, "Error finding max burnup"): + getOptimalAssemblyOrientation(self.assembly, self.assembly) + + def test_mismatchPinPowersAndLocations(self): + """Require pin powers and locations to be have the same length.""" + powers = np.arange(self.N_PINS + 1) + burnups = np.arange(self.N_PINS) + self.setAssemblyPinBurnups(self.assembly, burnups) + self.setAssemblyPinPowers(self.assembly, powers) + with self.assertRaisesRegex( + ValueError, "Inconsistent pin powers and number of pins" + ): + getOptimalAssemblyOrientation(self.assembly, self.assembly) + + +class TestFuelHandlerMgmtTools(ShuffleAndRotateTestHelper): + def test_buRotationWithFreshFeed(self): + """Test that rotation works if a new assembly is swapped with fresh fuel. + + Fresh feed assemblies will not exist in the reactor, and various checks that + try to the "previous" assembly's location can fail. + """ + newSettings = { + "fluxRecon": True, + "assemblyRotationAlgorithm": "buReducingAssemblyRotation", + } self.o.cs = self.o.cs.modified(newSettings=newSettings) - assem = self.o.r.core.getFirstAssembly(Flags.FUEL) - # apply dummy pin-level data to allow intelligent rotation - for b in assem.getBlocks(Flags.FUEL): - b.initializePinLocations() - b.p.percentBuMaxPinLocation = 10 - b.p.percentBuMax = 5 - b.p.linPowByPin = list(reversed(range(b.getNumPins()))) + fresh = self.r.core.createFreshFeed(self.o.cs) + self.assertEqual(fresh.lastLocationLabel, HexAssembly.LOAD_QUEUE) + fh = MockFuelHandler(self.o) + fh.chooseSwaps = mock.Mock(side_effect=lambda _: fh.moved.append(fresh)) - addSomeDetailAssemblies(hist, [assem]) - rotNum = b.getRotationNum() - rotAlgos.buReducingAssemblyRotation(fh) - self.assertNotEqual(b.getRotationNum(), rotNum) + with mock.patch( + "armi.physics.fuelCycle.assemblyRotationAlgorithms.getOptimalAssemblyOrientation", + ) as p: + fh.outage() + # The only moved assembly was most recently outside the core so we have no need to rotate + # Make sure our fake chooseSwaps added the fresh assembly to the moved assemblies + fh.chooseSwaps.assert_called_once() + p.assert_not_called() + + def test_buRotationWithStationaryRotation(self): + """Test that the burnup equalizing rotation algorithm works on non-shuffled assemblies.""" + newSettings = { + CONF_ASSEM_ROTATION_STATIONARY: True, + "fluxRecon": True, + "assemblyRotationAlgorithm": "buReducingAssemblyRotation", + } + self.o.cs = self.o.cs.modified(newSettings=newSettings) + + # Grab two assemblies that were not moved. One of which will have the detailed information + # needed for rotation + detailedAssem, coarseAssem = self.o.r.core.getChildrenWithFlags(Flags.FUEL)[:2] + self.setAssemblyPinBurnups(detailedAssem, burnups=np.arange(self.N_PINS)) + self.setAssemblyPinPowers(detailedAssem, pinPowers=np.arange(self.N_PINS)) + detailedAssem.rotate = mock.Mock() + coarseAssem.rotate = mock.Mock() + + fh = MockFuelHandler(self.o) + + with mock.patch( + "armi.physics.fuelCycle.assemblyRotationAlgorithms.getOptimalAssemblyOrientation", + return_value=5, + ) as p: + fh.outage() + p.assert_called_once_with(detailedAssem, detailedAssem) + # Assembly with detailed pin powers and pin burnups will be rotated + detailedAssem.rotate.assert_called_once() + self.compareMockedToExpectedRotation(5, detailedAssem.rotate) + # Assembly without pin level data will not be rotated + coarseAssem.rotate.assert_not_called() + + def test_rotateInShuffleQueue(self): + """Test for expected behavior when multiple assemblies are shuffled and rotated in one outage. + + Examine the behavior of three assemblies: ``first -> second -> third`` + + 1. ``first`` is moved to the location of ``second`` and rotated by comparing + ``first`` burnup against ``second`` pin powers. + 2. ``second`` is moved to the location of ``third`` and rotated by comparing + ``second`` burnup against ``third`` pin powers. + + where: + + * ``first`` burnup is maximized in the upper left direction. + * ``second`` pin power is minimized along the lower left direction. + * ``second`` burnup is maximized in the upper right direction. + * ``third`` pin power is minimized in the direct right direction. + + We should expect: + + 1. ``first`` is rotated from upper left to lower left => two 60 degree CCW rotations. + 2. ``second`` is rotated from upper right to direct right => five 60 degree CCW rotations. + """ + newSettings = { + CONF_ASSEM_ROTATION_STATIONARY: False, + "fluxRecon": True, + "assemblyRotationAlgorithm": "buReducingAssemblyRotation", + } + self.o.cs = self.o.cs.modified(newSettings=newSettings) + + first, second, third = self.r.core.getChildrenWithFlags(Flags.FUEL)[:3] + + firstBurnups = self.burnupWithMaxValue(_PinLocations.UPPER_LEFT_VERTEX) + self.setAssemblyPinBurnups(first, firstBurnups) + + secondPowers = self.powerWithMinValue(_PinLocations.LOWER_LEFT_VERTEX) + self.setAssemblyPinPowers(second, pinPowers=secondPowers) + + secondBurnups = self.burnupWithMaxValue(_PinLocations.UPPER_RIGHT_VERTEX) + self.setAssemblyPinBurnups(second, burnups=secondBurnups) + + thirdPowers = self.powerWithMinValue(_PinLocations.DUE_RIGHT_VERTEX) + self.setAssemblyPinPowers(third, thirdPowers) + + # Set the shuffling sequence + # first -> second + # second -> third + second.lastLocationLabel = first.getLocation() + third.lastLocationLabel = second.getLocation() + + first.rotate = mock.Mock() + second.rotate = mock.Mock() + third.rotate = mock.Mock() + + fh = MockFuelHandler(self.o) + fh.chooseSwaps = mock.Mock( + side_effect=lambda _: fh.moved.extend([second, third]) + ) + fh.outage() + + first.rotate.assert_called_once() + self.compareMockedToExpectedRotation(2, first.rotate, "First") + second.rotate.assert_called_once() + self.compareMockedToExpectedRotation(5, second.rotate, "Second") + third.rotate.assert_not_called() + + +class SimpleRotationTests(ShuffleAndRotateTestHelper): + """Test the simple rotation where assemblies are rotated a fixed amount.""" def test_simpleAssemblyRotation(self): - """Test rotating assemblies 120 degrees.""" + """Test rotating assemblies 120 degrees with two rotation events.""" fh = fuelHandlers.FuelHandler(self.o) newSettings = {CONF_ASSEM_ROTATION_STATIONARY: True} self.o.cs = self.o.cs.modified(newSettings=newSettings) diff --git a/armi/physics/fuelCycle/tests/test_fuelHandlers.py b/armi/physics/fuelCycle/tests/test_fuelHandlers.py index 8fd0f1f4b..5248aa4de 100644 --- a/armi/physics/fuelCycle/tests/test_fuelHandlers.py +++ b/armi/physics/fuelCycle/tests/test_fuelHandlers.py @@ -41,9 +41,7 @@ from armi.reactor.tests import test_reactors from armi.reactor.zones import Zone from armi.settings import caseSettings -from armi.tests import ArmiTestHelper -from armi.tests import mockRunLogs -from armi.tests import TEST_ROOT +from armi.tests import TEST_ROOT, ArmiTestHelper, mockRunLogs from armi.utils import directoryChangers @@ -121,6 +119,7 @@ def setUp(self): self.refAssembly = copy.deepcopy(self.assembly) self.directoryChanger.open() + self.r.core.locateAllAssemblies() def tearDown(self): # clean up the test @@ -222,6 +221,8 @@ def test_outage(self, mockChooseSwaps): self.assertEqual(len(fh.moved), 0) def test_outageEdgeCase(self): + """Check that an error is raised if the list of moved assemblies is invalid.""" + class MockFH(fuelHandlers.FuelHandler): def chooseSwaps(self, factor=1.0): self.moved = [None] diff --git a/armi/physics/fuelCycle/tests/test_utils.py b/armi/physics/fuelCycle/tests/test_utils.py new file mode 100644 index 000000000..a03eee700 --- /dev/null +++ b/armi/physics/fuelCycle/tests/test_utils.py @@ -0,0 +1,169 @@ +# Copyright 2024 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import copy +from unittest import TestCase + +import numpy as np + +from armi.physics.fuelCycle import utils +from armi.reactor.blocks import Block +from armi.reactor.components import Circle +from armi.reactor.flags import Flags +from armi.reactor.grids import IndexLocation, MultiIndexLocation + + +class FuelCycleUtilsTests(TestCase): + """Tests for geometry indifferent fuel cycle routines.""" + + N_PINS = 169 + + def setUp(self): + self.block = Block("test block") + self.fuel = Circle( + "test pin", + material="UO2", + Tinput=20, + Thot=20, + mult=self.N_PINS, + id=0.0, + od=1.0, + ) + + clad = Circle( + "clad", + material="HT9", + Tinput=20, + Thot=300, + id=1.0, + od=1.1, + ) + self.block.add(self.fuel) + self.block.add(clad) + # Force no fuel flags + self.fuel.p.flags = Flags.PIN + + def test_maxBurnupLocationFromComponents(self): + """Test that the ``Component.p.pinPercentBu`` parameter can reveal max burnup location.""" + self.fuel.spatialLocator = MultiIndexLocation(None) + locations = [] + for i in range(self.N_PINS): + loc = IndexLocation(i, 0, 0, None) + self.fuel.spatialLocator.append(loc) + locations.append(loc) + self.fuel.p.pinPercentBu = np.ones(self.N_PINS, dtype=float) + + # Pick an arbitrary index for the pin with the most burnup + maxBuIndex = self.N_PINS // 3 + self.fuel.p.pinPercentBu[maxBuIndex] *= 2 + expectedLoc = locations[maxBuIndex] + actual = utils.maxBurnupLocator(self.block) + self.assertEqual(actual, expectedLoc) + + def test_singleLocatorWithBurnup(self): + """Test that a single component with burnup can be used to find the highest burnup.""" + freeComp = Circle( + "free fuel", material="UO2", Tinput=200, Thot=200, id=0, od=1, mult=1 + ) + freeComp.spatialLocator = IndexLocation(2, 4, 0, None) + freeComp.p.pinPercentBu = [ + 0.01, + ] + loc = utils.maxBurnupLocator([freeComp]) + self.assertIs(loc, freeComp.spatialLocator) + + def test_maxBurnupLocatorWithNoBurnup(self): + """Ensure we catch an error if no burnup is found across components.""" + with self.assertRaisesRegex(ValueError, "No burnups found"): + utils.maxBurnupLocator([]) + + def test_maxBurnupLocatorMismatchedData(self): + """Ensure pin burnup and locations must agree.""" + freeComp = Circle( + "free fuel", material="UO2", Tinput=200, Thot=200, id=0, od=1, mult=1 + ) + freeComp.spatialLocator = IndexLocation(2, 4, 0, None) + freeComp.p.pinPercentBu = [ + 0.01, + 0.02, + ] + with self.assertRaisesRegex(ValueError, "Pin burnup.*pin locations.*differ"): + utils.maxBurnupLocator([freeComp]) + + def test_assemblyHasPinPower(self): + """Test the ability to check if an assembly has fuel pin powers.""" + fakeAssem = [self.block] + # No fuel blocks, no pin power on blocks => no pin powers + self.assertFalse(utils.assemblyHasFuelPinBurnup(fakeAssem)) + + # Yes fuel blocks, no pin power on blocks => no pin powers + self.block.p.flags |= Flags.FUEL + self.assertFalse(utils.assemblyHasFuelPinPowers(fakeAssem)) + + # Yes fuel blocks, yes pin power on blocks => yes pin powers + self.block.p.linPowByPin = np.arange(self.N_PINS, dtype=float) + self.assertTrue(utils.assemblyHasFuelPinPowers(fakeAssem)) + + # Yes fuel blocks, yes pin power assigned but all zeros => no pin powers + self.block.p.linPowByPin = np.zeros(self.N_PINS, dtype=float) + self.assertFalse(utils.assemblyHasFuelPinPowers(fakeAssem)) + + def test_assemblyHasPinBurnups(self): + """Test the ability to check if an assembly has fuel pin burnup.""" + fakeAssem = [self.block] + # No fuel components => no assembly burnups + self.assertFalse(self.block.getChildrenWithFlags(Flags.FUEL)) + self.assertFalse(utils.assemblyHasFuelPinBurnup(fakeAssem)) + # No fuel with burnup => no assembly burnups + self.block.p.flags |= Flags.FUEL + self.fuel.p.flags |= Flags.FUEL + self.assertFalse(utils.assemblyHasFuelPinBurnup(fakeAssem)) + # Fuel pin has burnup => yes assembly burnup + self.fuel.p.pinPercentBu = np.arange(self.N_PINS, dtype=float) + self.assertTrue(utils.assemblyHasFuelPinBurnup(fakeAssem)) + # Fuel pin has empty burnup => no assembly burnup + self.fuel.p.pinPercentBu = np.zeros(self.N_PINS) + self.assertFalse(utils.assemblyHasFuelPinBurnup(fakeAssem)) + # Yes burnup but no fuel flags => no assembly burnup + self.fuel.p.flags ^= Flags.FUEL + self.assertFalse(self.fuel.hasFlags(Flags.FUEL)) + self.fuel.p.pinPercentBu = np.arange(self.N_PINS, dtype=float) + self.assertFalse(utils.assemblyHasFuelPinBurnup(fakeAssem)) + + def test_maxBurnupBlock(self): + """Test the ability to find maximum burnup block in an assembly.""" + reflector = Block("reflector") + assem = [reflector, self.block] + self.block.p.percentBuPeak = 40 + expected = utils.maxBurnupBlock(assem) + self.assertIs(expected, self.block) + + # add a new block with more burnup higher up the stack + hotter = copy.deepcopy(self.block) + hotter.p.percentBuPeak *= 2 + expected = utils.maxBurnupBlock( + [reflector, self.block, hotter, self.block, reflector] + ) + self.assertIs(expected, hotter) + + def test_maxBurnupBlockNoBlocks(self): + """Ensure a more helpful error is provided for empty sequence.""" + with self.assertRaisesRegex(ValueError, "Error finding max burnup"): + utils.maxBurnupBlock([]) + + def test_maxBurnupBlockNoBurnup(self): + """Ensure that we will not return a block with zero burnup.""" + self.block.p.percentBuPeak = 0.0 + with self.assertRaisesRegex(ValueError, "Error finding max burnup"): + utils.maxBurnupBlock([self.block]) diff --git a/armi/physics/fuelCycle/utils.py b/armi/physics/fuelCycle/utils.py new file mode 100644 index 000000000..e0f9c14c3 --- /dev/null +++ b/armi/physics/fuelCycle/utils.py @@ -0,0 +1,142 @@ +# Copyright 2024 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""Geometric agnostic routines that are useful for fuel cycle analysis on pin-type reactors.""" + +import operator +import typing + +import numpy as np + +from armi import runLog +from armi.reactor.flags import Flags +from armi.reactor.grids import IndexLocation, MultiIndexLocation + +if typing.TYPE_CHECKING: + from armi.reactor.blocks import Block + from armi.reactor.components import Component + + +def assemblyHasFuelPinPowers(a: typing.Iterable["Block"]) -> bool: + """Determine if an assembly has pin powers. + + These are necessary for determining rotation and may or may + not be present on all assemblies. + + Parameters + ---------- + a : Assembly + Assembly in question + + Returns + ------- + bool + If at least one fuel block in the assembly has pin powers. + """ + # Avoid using Assembly.getChildrenWithFlags(Flags.FUEL) + # because that creates an entire list where we may just need the first + # fuel block + fuelBlocks = filter(lambda b: b.hasFlags(Flags.FUEL), a) + return any(b.hasFlags(Flags.FUEL) and np.any(b.p.linPowByPin) for b in fuelBlocks) + + +def assemblyHasFuelPinBurnup(a: typing.Iterable["Block"]) -> bool: + """Determine if an assembly has pin burnups. + + These are necessary for determining rotation and may or may not + be present on all assemblies. + + Parameters + ---------- + a : Assembly + Assembly in question + + Returns + ------- + bool + If a block with pin burnup was found. + + Notes + ----- + Checks if any `Component.p.pinPercentBu`` is set and contains non-zero data + on a fuel component in the block. + """ + # Avoid using Assembly.getChildrenWithFlags(Flags.FUEL) + # because that creates an entire list where we may just need the first + # fuel block. Same for avoiding Block.getChildrenWithFlags. + hasFuelFlags = lambda o: o.hasFlags(Flags.FUEL) + for b in filter(hasFuelFlags, a): + for c in filter(hasFuelFlags, b): + if np.any(c.p.pinPercentBu): + return True + return False + + +def maxBurnupLocator( + children: typing.Iterable["Component"], +) -> IndexLocation: + """Find the location of the pin with highest burnup by looking at components. + + Parameters + ---------- + children : iterable[Component] + Iterator over children with a spatial locator and ``pinPercentBu`` parameter + + Returns + ------- + IndexLocation + Location of the pin with the highest burnup. + + Raises + ------ + ValueError + If no children have burnup, or the burnup and locators differ. + """ + maxBu = 0 + maxLocation = None + withBurnupAndLocs = filter( + lambda c: c.spatialLocator is not None and c.p.pinPercentBu is not None, + children, + ) + for child in withBurnupAndLocs: + pinBu = child.p.pinPercentBu + if isinstance(child.spatialLocator, MultiIndexLocation): + locations = child.spatialLocator + else: + locations = [child.spatialLocator] + if len(locations) != pinBu.size: + raise ValueError( + f"Pin burnup (n={len(locations)}) and pin locations (n={pinBu.size}) " + f"on {child} differ: {locations=} :: {pinBu=}" + ) + myMaxIX = pinBu.argmax() + myMaxBu = pinBu[myMaxIX] + if myMaxBu > maxBu: + maxBu = myMaxBu + maxLocation = locations[myMaxIX] + if maxLocation is not None: + return maxLocation + raise ValueError("No burnups found!") + + +def maxBurnupBlock(a: typing.Iterable["Block"]) -> "Block": + """Find the block that contains the pin with the highest burnup.""" + buGetter = operator.attrgetter("p.percentBuPeak") + # Discard any blocks with zero burnup + blocksWithBurnup = filter(buGetter, a) + try: + return max(blocksWithBurnup, key=buGetter) + except Exception as ee: + msg = f"Error finding max burnup block from {a}" + runLog.error(msg) + raise ValueError(msg) from ee diff --git a/doc/release/0.5.rst b/doc/release/0.5.rst index 91a786eea..0b92c954d 100644 --- a/doc/release/0.5.rst +++ b/doc/release/0.5.rst @@ -81,6 +81,7 @@ API Changes #. History Tracker: "detail assemblies" are now fuel and control assemblies. (`PR#1990 `_) #. Removing ``Block.breakFuelComponentsIntoIndividuals()``. (`PR#1990 `_) #. Moving ``getPuMoles`` from blocks.py up to composites.py. (`PR#1990 `_) +#. Requiring ``buReducingAssemblyRotation`` and ``getOptimalAssemblyOrientation`` to have pin-level burnup. (`PR#2019 `_) Bug Fixes ---------