Skip to content

Commit

Permalink
Fix _swapFluxParam for non-uniform assemblies (#577)
Browse files Browse the repository at this point in the history
  • Loading branch information
@john-science
john-science committed Feb 16, 2022
1 parent 4316eb7 commit b3cf31e
Show file tree
Hide file tree
Showing 13 changed files with 728 additions and 140 deletions.
97 changes: 72 additions & 25 deletions armi/physics/fuelCycle/fuelHandlers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@
from armi.utils import directoryChangers, pathTools
from armi import utils
from armi.utils import plotting
from armi.utils.mathematics import resampleStepwise

runLog = logging.getLogger(__name__)

Expand Down Expand Up @@ -1360,12 +1361,12 @@ def swapAssemblies(self, a1, a2):
--------
dischargeSwap : swap assemblies where one is outside the core and the other is inside
"""

if a1 is None or a2 is None:
runLog.warning(
"Cannot swap None assemblies. Check your findAssembly results. Skipping swap"
)
return

runLog.extra("Swapping {} with {}.".format(a1, a2))
# add assemblies into the moved location
for a in [a1, a2]:
Expand Down Expand Up @@ -1460,29 +1461,78 @@ def _swapFluxParam(self, incoming, outgoing):
Assumes assemblies have the same block structure. If not, blocks will be swapped one-for-one until
the shortest one has been used up and then the process will truncate.
"""
if len(incoming) != len(outgoing):
runLog.warning(
"{0} and {1} have different numbers of blocks. Flux swapping (for XS weighting) will "
"be questionable".format(incoming, outgoing)
# Find the block-based mesh points for each assembly
meshIn = self.r.core.findAllAxialMeshPoints([incoming], False)
meshOut = self.r.core.findAllAxialMeshPoints([outgoing], False)

# If the assembly mesh points don't match, the swap won't be easy
if meshIn != meshOut:
runLog.debug(
"{0} and {1} have different meshes, resampling.".format(
incoming, outgoing
)
)

# grab the current values for incoming and outgoing
fluxIn = [b.p.flux for b in incoming]
mgFluxIn = [b.p.mgFlux for b in incoming]
powerIn = [b.p.power for b in incoming]
fluxOut = [b.p.flux for b in outgoing]
mgFluxOut = [b.p.mgFlux for b in outgoing]
powerOut = [b.p.power for b in outgoing]

# resample incoming to outgoing, and vice versa
fluxOutNew = resampleStepwise(meshIn, fluxIn, meshOut)
mgFluxOutNew = resampleStepwise(meshIn, mgFluxIn, meshOut)
powerOutNew = resampleStepwise(meshIn, powerIn, meshOut, avg=False)
fluxInNew = resampleStepwise(meshOut, fluxOut, meshIn)
mgFluxInNew = resampleStepwise(meshOut, mgFluxOut, meshIn)
powerInNew = resampleStepwise(meshOut, powerOut, meshIn, avg=False)

# load the new outgoing values into place
for b, flux, mgFlux, power in zip(
outgoing, fluxOutNew, mgFluxOutNew, powerOutNew
):
b.p.flux = flux
b.p.mgFlux = mgFlux
b.p.power = power
b.p.pdens = power / b.getVolume()

# load the new incoming values into place
for b, flux, mgFlux, power in zip(
incoming, fluxInNew, mgFluxInNew, powerInNew
):
b.p.flux = flux
b.p.mgFlux = mgFlux
b.p.power = power
b.p.pdens = power / b.getVolume()

return

# Since the axial mesh points match, do the simple swap
for bi, (bIncoming, bOutgoing) in enumerate(zip(incoming, outgoing)):
if (
bi not in self.cs["stationaryBlocks"]
): # stationary blocks are already swapped.
incomingFlux = bIncoming.p.flux
incomingMgFlux = bIncoming.p.mgFlux
incomingPower = bIncoming.p.power
outgoingFlux = bOutgoing.p.flux
outgoingMgFlux = bOutgoing.p.mgFlux
outgoingPower = bOutgoing.p.power
if outgoingFlux > 0.0:
bIncoming.p.flux = outgoingFlux
bIncoming.p.mgFlux = outgoingMgFlux
bIncoming.p.power = outgoingPower
if incomingFlux > 0.0:
bOutgoing.p.flux = incomingFlux
bOutgoing.p.mgFlux = incomingMgFlux
bOutgoing.p.power = incomingPower
if bi in self.cs["stationaryBlocks"]:
# stationary blocks are already swapped
continue

incomingFlux = bIncoming.p.flux
incomingMgFlux = bIncoming.p.mgFlux
incomingPower = bIncoming.p.power
outgoingFlux = bOutgoing.p.flux
outgoingMgFlux = bOutgoing.p.mgFlux
outgoingPower = bOutgoing.p.power

if outgoingFlux > 0.0:
bIncoming.p.flux = outgoingFlux
bIncoming.p.mgFlux = outgoingMgFlux
bIncoming.p.power = outgoingPower
bIncoming.p.pdens = outgoingPower / bIncoming.getVolume()

if incomingFlux > 0.0:
bOutgoing.p.flux = incomingFlux
bOutgoing.p.mgFlux = incomingMgFlux
bOutgoing.p.power = incomingPower
bOutgoing.p.pdens = incomingPower / bOutgoing.getVolume()

def swapCascade(self, assemList):
"""
Expand All @@ -1499,7 +1549,6 @@ def swapCascade(self, assemList):
[A <- B <- C <- D]
"""

# first check for duplicates
for assem in assemList:
if assemList.count(assem) != 1:
Expand Down Expand Up @@ -1592,7 +1641,6 @@ def readMoves(fname):
makeShuffleReport : writes the file that is read here.
"""

try:
f = open(fname)
except:
Expand Down Expand Up @@ -1906,7 +1954,6 @@ def doRepeatShuffle(
-----
This is a helper function for repeatShufflePattern
"""

moved = []

# shuffle all of the load chain assemblies (These include discharges to SFP
Expand Down
177 changes: 168 additions & 9 deletions armi/physics/fuelCycle/tests/test_fuelHandlers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,7 @@ def tearDown(self):

self.directoryChanger.close()

def test_FindHighBu(self):
def test_findHighBu(self):
loc = self.r.core.spatialGrid.getLocatorFromRingAndPos(5, 4)
a = self.r.core.childrenByLocator[loc]
# set burnup way over 1.0, which is otherwise the highest bu in the core
Expand All @@ -133,7 +133,7 @@ def test_FindHighBu(self):
)
self.assertIs(a, a1)

def test_Width(self):
def test_width(self):
"""Tests the width capability of findAssembly."""
fh = fuelHandlers.FuelHandler(self.o)
assemsByRing = collections.defaultdict(list)
Expand Down Expand Up @@ -224,10 +224,8 @@ def test_Width(self):
"The lowest power ring returned is {0}. It should be {1}".format(ring, 1),
)

def test_FindMany(self):
r"""
Tests the findMany and type aspects of the fuel handler
"""
def test_findMany(self):
"""Tests the findMany and type aspects of the fuel handler"""
fh = fuelHandlers.FuelHandler(self.o)

igniters = fh.findAssembly(typeSpec=Flags.IGNITER | Flags.FUEL, findMany=True)
Expand Down Expand Up @@ -256,9 +254,7 @@ def test_FindMany(self):
)

def test_findInSFP(self):
r"""
Tests ability to pull from the spent fuel pool.
"""
"""Tests ability to pull from the spent fuel pool"""
fh = fuelHandlers.FuelHandler(self.o)
spent = fh.findAssembly(
findMany=True,
Expand Down Expand Up @@ -531,6 +527,169 @@ def test_buildEqRingSchedule(self):
locSchedule = fh.buildEqRingSchedule([2, 1])
self.assertEqual(locSchedule, ["002-001", "002-002", "001-001"])

def test_swapFluxParamSameLength(self):
"""
Test the _swapFluxParams method for the usual case,
where each of the input assembles have the same number of assemblies,
on the same mesh
"""
# grab the assemblies
assems = self.r.core.getAssemblies(Flags.FEED)
self.assertEqual(len(assems), 14)

for a in assems:
self.assertEqual(len(a.getBlocks()), 5)

# make two copies of an arbitraty assembly
a1 = copy.deepcopy(list(assems)[1])
a2 = copy.deepcopy(list(assems)[1])
blocks1 = list(a1.getBlocks())
blocks2 = list(a2.getBlocks())
self.assertEqual(len(blocks1), 5)
self.assertEqual(len(blocks2), 5)
self.assertEqual(blocks1[3].p.height, 25)
self.assertEqual(blocks2[3].p.height, 25)

# 1. alter the values of a single block in assembly 2
b2 = list(blocks2)[1]
b2.p.flux = b2.p.flux * 2
b2.p.power = 1000
b2.p.pdens = b2.p.power / b2.getVolume()

# grab the power before the swap
power1 = sum([b.p.power for b in a1.getBlocks()])
power2 = sum([b.p.power for b in a2.getBlocks()])

# 2. validate the situation is as you'd expect
self.assertEqual(list(a1.getBlocks())[1].p.flux, 50000000000.0)
self.assertEqual(list(a2.getBlocks())[1].p.flux, 100000000000.0)
self.assertEqual(list(a1.getBlocks())[1].p.power, 0.0)
self.assertEqual(list(a2.getBlocks())[1].p.power, 1000.0)
self.assertEqual(list(a1.getBlocks())[1].p.pdens, 0.0)
self.assertGreater(list(a2.getBlocks())[1].p.pdens, 0.0)

# 3. do the swap
fh = fuelHandlers.FuelHandler(self.o)
fh._swapFluxParam(a1, a2)

# 4. validate the swap worked
self.assertEqual(list(a1.getBlocks())[1].p.flux, 100000000000.0)
self.assertEqual(list(a2.getBlocks())[1].p.flux, 50000000000.0)
self.assertEqual(list(a1.getBlocks())[1].p.power, 1000.0)
self.assertEqual(list(a2.getBlocks())[1].p.power, 0.0)
self.assertGreater(list(a1.getBlocks())[1].p.pdens, 0.0)
self.assertEqual(list(a2.getBlocks())[1].p.pdens, 0.0)
self.assertEqual(sum([b.p.power for b in a1.getBlocks()]), power2)
self.assertEqual(sum([b.p.power for b in a2.getBlocks()]), power1)

def test_swapFluxParamDifferentLengths(self):
"""
Test the _swapFluxParams method for the less common, and more complicated case,
where the input assembles have different numbers of blocks, potentially on
wildly different point meshes.
"""
# grab the assemblies
assems = self.r.core.getAssemblies(Flags.FEED)

# make two copies of an arbitraty assembly
a1 = copy.deepcopy(list(assems)[1])
a2 = copy.deepcopy(list(assems)[1])
height2 = 25.0
self.assertEqual(list(a1.getBlocks())[3].p.height, height2)
self.assertEqual(list(a2.getBlocks())[3].p.height, height2)

# grab the blocks from the second assembly
blocks2 = list(a2.getBlocks())
self.assertEqual(len(blocks2), 5)

# grab a single block from the second assembly, to be altered
b2 = list(blocks2)[1]
self.assertEqual(b2.p.height, height2)
self.assertEqual(b2.p.flux, 50000000000.0)
self.assertIsNone(b2.p.mgFlux)
self.assertEqual(b2.p.power, 0.0)
self.assertEqual(b2.p.pdens, 0.0)
volume2 = 6074.356
self.assertAlmostEqual(b2.getVolume(), volume2, delta=0.1)

# split the the block into two of half the heights
b20 = copy.deepcopy(b2)
b21 = copy.deepcopy(b2)
b20.setHeight(height2 / 2)
b21.setHeight(height2 / 2)
self.assertAlmostEqual(b20.getVolume(), volume2 / 2, delta=0.1)
self.assertAlmostEqual(b21.getVolume(), volume2 / 2, delta=0.1)

# give the two new (smaller) blocks some power/pdens
b20.p.power = 1000
b21.p.power = 2000
b20.p.pdens = b20.p.power / b20.getVolume()
b21.p.pdens = b21.p.power / b21.getVolume()
self.assertEqual(b20.p.power, 1000.0)
self.assertEqual(b21.p.power, 2000.0)
self.assertAlmostEqual(b20.p.pdens, 0.3292, delta=0.1)
self.assertAlmostEqual(b21.p.pdens, 0.6585, delta=0.1)

# give the second assembly the new blocks
a2.removeAll()
a2.setChildren([blocks2[0]] + [b20, b21] + blocks2[2:])

# validate the situation is as you'd expect
self.assertEqual(len(a1.getBlocks()), 5)
self.assertEqual(len(a2.getBlocks()), 6)

# validate the power before the swap
power1 = [b.p.power for b in a1.getBlocks()]
power2 = [b.p.power for b in a2.getBlocks()]

self.assertEqual(power1, [0, 0, 0, 0, 0])
self.assertEqual(power2, [0, 1000, 2000, 0, 0, 0])

# validate the power density before the swap
for b in a1.getBlocks():
self.assertEqual(b.p.pdens, 0.0)

pdens2i = [0, 0.32925299379047496, 0.6585059875809499, 0, 0, 0]
for i, b in enumerate(a2.getBlocks()):
self.assertAlmostEqual(b.p.pdens, pdens2i[i], msg=i)

# validate the flux before the swap
for b in a1.getBlocks():
self.assertEqual(b.p.flux, 50000000000.0)

for b in a2.getBlocks():
self.assertEqual(b.p.flux, 50000000000.0)

# do the swap, using averages
fh = fuelHandlers.FuelHandler(self.o)
fh._swapFluxParam(a1, a2)

# grab the power after the swap
power1f = [b.p.power for b in a1.getBlocks()]
power2f = [b.p.power for b in a2.getBlocks()]

# validate the swap worked
self.assertEqual(len(a1.getBlocks()), 5)
self.assertEqual(len(a2.getBlocks()), 6)

self.assertEqual(power1f, [0, 3000, 0, 0, 0])
self.assertEqual(power2f, [0, 0, 0, 0, 0, 0])

# validate the power density after the swap
pdens1f = [0, 0.4938794906857124, 0, 0, 0]
for i, b in enumerate(a1.getBlocks()):
self.assertAlmostEqual(b.p.pdens, pdens1f[i], msg=i)

for i, b in enumerate(a2.getBlocks()):
self.assertAlmostEqual(b.p.pdens, 0, msg=i)

# validate the flux after the swap
for b in a1.getBlocks():
self.assertEqual(b.p.flux, 50000000000.0)

for b in a2.getBlocks():
self.assertEqual(b.p.flux, 50000000000.0)


class TestFuelPlugin(unittest.TestCase):
"""Tests that make sure the plugin is being discovered well."""
Expand Down
Loading

0 comments on commit b3cf31e

Please sign in to comment.