From 65574120048c7d6ee0d523d667e8d2f49c45015f Mon Sep 17 00:00:00 2001 From: Chris Keckler Date: Fri, 30 Sep 2022 17:00:33 -0700 Subject: [PATCH] Fix burnup-dependent reaction rates (#911) --- armi/reactor/blocks.py | 33 -------- armi/reactor/components/component.py | 77 ------------------- armi/reactor/composites.py | 105 +++++++++++++++++++++++--- armi/reactor/tests/test_blocks.py | 33 ++++---- armi/reactor/tests/test_components.py | 12 --- armi/reactor/tests/test_composites.py | 12 +++ doc/release/0.2.rst | 1 + 7 files changed, 128 insertions(+), 145 deletions(-) diff --git a/armi/reactor/blocks.py b/armi/reactor/blocks.py index 8aa3822e9..ce25084e4 100644 --- a/armi/reactor/blocks.py +++ b/armi/reactor/blocks.py @@ -1496,39 +1496,6 @@ def getLumpedFissionProductCollection(self): """ return composites.ArmiObject.getLumpedFissionProductCollection(self) - def getReactionRates(self, nucName, nDensity=None): - """ - Parameters - ---------- - nucName - str - nuclide name -- e.g. 'U235' - nDensity - float - number Density - - Returns - ------- - rxnRates : dict - dictionary of reaction rates (rxn/s) for nG, nF, n2n, nA and nP - - Note - ---- - If you set nDensity to 1/CM2_PER_BARN this makes 1 group cross section generation easier - """ - if nDensity is None: - nDensity = self.getNumberDensity(nucName) - try: - return components.component.getReactionRateDict( - nucName, - self.core.lib, - self.p.xsType, - self.getIntegratedMgFlux(), - nDensity, - ) - except AttributeError: - # AttributeError because there was no library because no parent.r -- this is a armiObject without flux so - # send it some zeros - return {"nG": 0, "nF": 0, "n2n": 0, "nA": 0, "nP": 0} - def rotate(self, deg): """Function for rotating a block's spatially varying variables by a specified angle. diff --git a/armi/reactor/components/component.py b/armi/reactor/components/component.py index 0f25d0250..1ee28ed47 100644 --- a/armi/reactor/components/component.py +++ b/armi/reactor/components/component.py @@ -1187,41 +1187,6 @@ def getLumpedFissionProductCollection(self): else: return composites.ArmiObject.getLumpedFissionProductCollection(self) - def getReactionRates(self, nucName, nDensity=None): - """ - Parameters - ---------- - nucName - str - nuclide name -- e.g. 'U235' - nDensity - float - number Density - - Returns - ------- - rxnRates - dict - dictionary of reaction rates (rxn/s) for nG, nF, n2n, nA and nP - - Note - ---- - if you set nDensity to 1/CM2_PER_BARN this makes 1 group cross section generation easier - - """ - if nDensity is None: - nDensity = self.getNumberDensity(nucName) - try: - return getReactionRateDict( - nucName, - self.getAncestorWithFlags(Flags.CORE).lib, - self.parent.p.xsType, - self.getIntegratedMgFlux(), - nDensity, - ) - except (AttributeError, KeyError): - # AttributeError because there was no library because no parent.r -- this is a armiObject without flux so - # send it some zeros - # KeyError because nucName was not in the library - return {"nG": 0, "nF": 0, "n2n": 0, "nA": 0, "nP": 0} - def getMicroSuffix(self): return self.parent.getMicroSuffix() @@ -1245,45 +1210,3 @@ class ShapedComponent(Component): """A component with well-defined dimensions.""" pass - - -def getReactionRateDict(nucName, lib, xsType, mgFlux, nDens): - """ - Parameters - ---------- - nucName - str - nuclide name -- e.g. 'U235', 'PU239', etc. Not to be confused with the nuclide - _label_, see the nucDirectory module for a description of the difference. - lib - isotxs - cross section library - xsType - str - cross section type -- e.g. - 'A' - mgFlux - numpy.nArray - integrated mgFlux (n-cm/s) - nDens - float - number density (at/bn-cm) - - Returns - ------- - rxnRates - dict - dictionary of reaction rates (rxn/s) for nG, nF, n2n, nA and nP - - Note - ---- - assume there is no n3n cross section in ISOTXS - - """ - nucLabel = nuclideBases.byName[nucName].label - key = "{}{}A".format(nucLabel, xsType) - libNuc = lib[key] - rxnRates = {"n3n": 0} - for rxName, mgXSs in [ - ("nG", libNuc.micros.nGamma), - ("nF", libNuc.micros.fission), - ("n2n", libNuc.micros.n2n), - ("nA", libNuc.micros.nalph), - ("nP", libNuc.micros.np), - ]: - rxnRates[rxName] = nDens * sum(mgXSs * mgFlux) - - return rxnRates diff --git a/armi/reactor/composites.py b/armi/reactor/composites.py index 31265b804..462bf64e2 100644 --- a/armi/reactor/composites.py +++ b/armi/reactor/composites.py @@ -47,6 +47,7 @@ from armi.physics.neutronics.fissionProductModel import fissionProductModel from armi.reactor import grids from armi.reactor import parameters + from armi.reactor.flags import Flags, TypeSpec from armi.reactor.parameters import resolveCollections from armi.utils import densityTools @@ -901,7 +902,10 @@ def getMassFrac(self, nucName): return sum(self.getMassFracs().get(nucName, 0.0) for nucName in nuclideNames) def getMicroSuffix(self): - pass + raise NotImplementedError( + f"Cannot get the suffix on {type(self)} objects. Only certain subclasses" + " of composite such as Blocks or Components have the concept of micro suffixes." + ) def _getNuclidesFromSpecifier(self, nucSpec): """ @@ -1492,10 +1496,6 @@ def getAncestorWithFlags(self, typeSpec: TypeSpec, exactMatch=False): armi.composites.ArmiObject the first ancestor up the chain of parents that matches the passed flags - Notes - ----- - This will throw an error if no ancestor can be found that matches the typeSpec - See Also -------- ArmiObject.hasFlags() @@ -3097,15 +3097,58 @@ def getIntegratedMgFlux(self, adjoint=False, gamma=False): ) return integratedMgFlux + def _getReactionRates(self, nucName, nDensity=None): + """ + Parameters + ---------- + nucName : str + nuclide name -- e.g. 'U235' + nDensity : float + number density + + Returns + ------- + rxnRates : dict + dictionary of reaction rates (rxn/s) for nG, nF, n2n, nA and nP + + Notes + ---- + If you set nDensity to 1/CM2_PER_BARN this makes 1 group cross section generation easier + """ + from armi.reactor.blocks import Block + + if nDensity is None: + nDensity = self.getNumberDensity(nucName) + try: + return getReactionRateDict( + nucName, + self.getAncestorWithFlags(Flags.CORE).lib, + self.getAncestor(lambda x: isinstance(x, Block)).getMicroSuffix(), + self.getIntegratedMgFlux(), + nDensity, + ) + except AttributeError: + runLog.warning( + f"Object {self} does not belong to a core and so has no reaction rates.", + single=True, + ) + return {"nG": 0, "nF": 0, "n2n": 0, "nA": 0, "nP": 0} + except KeyError: + runLog.warning( + f"Attempting to get a reaction rate on an isotope not in the lib {nucName}.", + single=True, + ) + return {"nG": 0, "nF": 0, "n2n": 0, "nA": 0, "nP": 0} + def getReactionRates(self, nucName, nDensity=None): """ Get the reaction rates of a certain nuclide on this object. Parameters ---------- - nucName - str + nucName : str nuclide name -- e.g. 'U235' - nDensity - float + nDensity : float number Density Returns @@ -3121,8 +3164,10 @@ def getReactionRates(self, nucName, nDensity=None): """ rxnRates = {"nG": 0, "nF": 0, "n2n": 0, "nA": 0, "nP": 0, "n3n": 0} - for armiObject in self: - for rxName, val in armiObject.getReactionRates( + # not all composite objects are iterable (i.e. components), so in that + # case just examine only the object itself + for armiObject in self.getChildren() or [self]: + for rxName, val in armiObject._getReactionRates( nucName, nDensity=nDensity ).items(): rxnRates[rxName] += val @@ -3292,3 +3337,45 @@ def getDominantMaterial( return samples[maxMatName] return None + + +def getReactionRateDict(nucName, lib, xsSuffix, mgFlux, nDens): + """ + Parameters + ---------- + nucName : str + nuclide name -- e.g. 'U235', 'PU239', etc. Not to be confused with the nuclide + _label_, see the nucDirectory module for a description of the difference. + lib : isotxs + cross section library + xsSuffix : str + cross section suffix, consisting of the type followed by the burnup group, + e.g. 'AB' for the second burnup group of type A + mgFlux : numpy.nArray + integrated mgFlux (n-cm/s) + nDens : float + number density (atom/bn-cm) + + Returns + ------- + rxnRates - dict + dictionary of reaction rates (rxn/s) for nG, nF, n2n, nA and nP + + Notes + ----- + Assume there is no n3n cross section in ISOTXS + """ + nucLabel = nuclideBases.byName[nucName].label + key = "{}{}".format(nucLabel, xsSuffix) + libNuc = lib[key] + rxnRates = {"n3n": 0} + for rxName, mgXSs in [ + ("nG", libNuc.micros.nGamma), + ("nF", libNuc.micros.fission), + ("n2n", libNuc.micros.n2n), + ("nA", libNuc.micros.nalph), + ("nP", libNuc.micros.np), + ]: + rxnRates[rxName] = nDens * sum(mgXSs * mgFlux) + + return rxnRates diff --git a/armi/reactor/tests/test_blocks.py b/armi/reactor/tests/test_blocks.py index 4cc28c963..4e7d49bdb 100644 --- a/armi/reactor/tests/test_blocks.py +++ b/armi/reactor/tests/test_blocks.py @@ -33,6 +33,8 @@ from armi.utils import hexagon, units from armi.utils.units import MOLES_PER_CC_TO_ATOMS_PER_BARN_CM +NUM_PINS_IN_TEST_BLOCK = 217 + def buildSimpleFuelBlock(): """Return a simple block containing fuel, clad, duct, and coolant.""" @@ -79,7 +81,7 @@ def loadTestBlock(cold=True): assemNum = 3 block = blocks.HexBlock("TestHexBlock") block.setType("defaultType") - block.p.nPins = 217 + block.p.nPins = NUM_PINS_IN_TEST_BLOCK assembly = makeTestAssembly(assemNum, 1, r=r) # NOTE: temperatures are supposed to be in C @@ -93,7 +95,7 @@ def loadTestBlock(cold=True): "Thot": hotTempFuel, "od": 0.84, "id": 0.6, - "mult": 217.0, + "mult": NUM_PINS_IN_TEST_BLOCK, } fuel = components.Circle("fuel", "UZr", **fuelDims) @@ -102,7 +104,7 @@ def loadTestBlock(cold=True): "Thot": hotTempCoolant, "od": "fuel.id", "id": 0.3, - "mult": 217.0, + "mult": NUM_PINS_IN_TEST_BLOCK, } bondDims["components"] = {"fuel": fuel} bond = components.Circle("bond", "Sodium", **bondDims) @@ -112,7 +114,7 @@ def loadTestBlock(cold=True): "Thot": hotTempStructure, "od": "bond.id", "id": 0.0, - "mult": 217.0, + "mult": NUM_PINS_IN_TEST_BLOCK, } annularVoidDims["components"] = {"bond": bond} annularVoid = components.Circle("annular void", "Void", **annularVoidDims) @@ -122,7 +124,7 @@ def loadTestBlock(cold=True): "Thot": hotTempStructure, "od": 0.90, "id": 0.85, - "mult": 217.0, + "mult": NUM_PINS_IN_TEST_BLOCK, } innerLiner = components.Circle("inner liner", "Graphite", **innerLinerDims) @@ -131,7 +133,7 @@ def loadTestBlock(cold=True): "Thot": hotTempStructure, "od": "inner liner.id", "id": "fuel.od", - "mult": 217.0, + "mult": NUM_PINS_IN_TEST_BLOCK, } fuelLinerGapDims["components"] = {"inner liner": innerLiner, "fuel": fuel} fuelLinerGap = components.Circle("gap1", "Void", **fuelLinerGapDims) @@ -141,7 +143,7 @@ def loadTestBlock(cold=True): "Thot": hotTempStructure, "od": 0.95, "id": 0.90, - "mult": 217.0, + "mult": NUM_PINS_IN_TEST_BLOCK, } outerLiner = components.Circle("outer liner", "HT9", **outerLinerDims) @@ -150,7 +152,7 @@ def loadTestBlock(cold=True): "Thot": hotTempStructure, "od": "outer liner.id", "id": "inner liner.od", - "mult": 217.0, + "mult": NUM_PINS_IN_TEST_BLOCK, } linerLinerGapDims["components"] = { "outer liner": outerLiner, @@ -163,7 +165,7 @@ def loadTestBlock(cold=True): "Thot": hotTempStructure, "od": 1.05, "id": 0.95, - "mult": 217.0, + "mult": NUM_PINS_IN_TEST_BLOCK, } cladding = components.Circle("clad", "HT9", **claddingDims) @@ -172,7 +174,7 @@ def loadTestBlock(cold=True): "Thot": hotTempStructure, "od": "clad.id", "id": "outer liner.od", - "mult": 217.0, + "mult": NUM_PINS_IN_TEST_BLOCK, } linerCladGapDims["components"] = {"clad": cladding, "outer liner": outerLiner} linerCladGap = components.Circle("gap3", "Void", **linerCladGapDims) @@ -184,7 +186,7 @@ def loadTestBlock(cold=True): "id": 0.0, "axialPitch": 30.0, "helixDiameter": 1.1, - "mult": 217.0, + "mult": NUM_PINS_IN_TEST_BLOCK, } wire = components.Helix("wire", "HT9", **wireDims) @@ -1574,13 +1576,16 @@ def test_getReactionRates(self): r.core.lib = isotxs.readBinary(ISOAA_PATH) block.p.mgFlux = 1 - self.assertEqual( + self.assertAlmostEqual( block.getReactionRates("PU239")["nG"], block.getNumberDensity("PU239") * sum(r.core.lib["PU39AA"].micros.nGamma), ) - with self.assertRaises(KeyError): - block.getReactionRates("PU39") + # the key is invalid, so should get back all zeros + self.assertEqual( + block.getReactionRates("PU39"), + {"nG": 0, "nF": 0, "n2n": 0, "nA": 0, "nP": 0, "n3n": 0}, + ) class Test_NegativeVolume(unittest.TestCase): diff --git a/armi/reactor/tests/test_components.py b/armi/reactor/tests/test_components.py index 8d8fe8b4b..5ce88a202 100644 --- a/armi/reactor/tests/test_components.py +++ b/armi/reactor/tests/test_components.py @@ -20,7 +20,6 @@ import math import unittest -from armi import nuclearDataIO from armi.reactor import components from armi.reactor.components import ( Component, @@ -46,8 +45,6 @@ ComponentType, ) from armi.reactor.components import materials -from armi.reactor.components.component import getReactionRateDict -from armi.tests import ISOAA_PATH from armi.utils import units @@ -1414,15 +1411,6 @@ def test_getEnrichment(self): self.assertAlmostEqual(self.fuel.getEnrichment(), 0.3) -class TestGetReactionRateDict(unittest.TestCase): - def test_getReactionRateDict(self): - lib = nuclearDataIO.isotxs.readBinary(ISOAA_PATH) - rxRatesDict = getReactionRateDict( - nucName="PU239", lib=lib, xsType="A", mgFlux=1, nDens=1 - ) - self.assertEqual(rxRatesDict["nG"], sum(lib["PU39AA"].micros.nGamma)) - - if __name__ == "__main__": # import sys; sys.argv = ['', 'TestMaterialAdjustments.test_adjustMassFrac_U235'] unittest.main() diff --git a/armi/reactor/tests/test_composites.py b/armi/reactor/tests/test_composites.py index f65a5d92b..d9fafbb1c 100644 --- a/armi/reactor/tests/test_composites.py +++ b/armi/reactor/tests/test_composites.py @@ -17,6 +17,7 @@ from copy import deepcopy import unittest +from armi import nuclearDataIO from armi import runLog from armi import settings from armi.nucDirectory import nucDir, nuclideBases @@ -28,11 +29,13 @@ from armi.reactor import composites from armi.reactor import assemblies from armi.reactor.components import basicShapes +from armi.reactor.composites import getReactionRateDict from armi.reactor import grids from armi.reactor.blueprints import assemblyBlueprint from armi.reactor import parameters from armi.reactor.flags import Flags from armi.reactor.tests.test_blocks import loadTestBlock +from armi.tests import ISOAA_PATH class MockBP: @@ -608,5 +611,14 @@ def test_copyParamsFrom(self): self.assertEqual(obj2.p.percentBu, self.obj.p.percentBu) +class TestGetReactionRateDict(unittest.TestCase): + def test_getReactionRateDict(self): + lib = nuclearDataIO.isotxs.readBinary(ISOAA_PATH) + rxRatesDict = getReactionRateDict( + nucName="PU239", lib=lib, xsSuffix="AA", mgFlux=1, nDens=1 + ) + self.assertEqual(rxRatesDict["nG"], sum(lib["PU39AA"].micros.nGamma)) + + if __name__ == "__main__": unittest.main() diff --git a/doc/release/0.2.rst b/doc/release/0.2.rst index 4a91e3fb0..dd256fe69 100644 --- a/doc/release/0.2.rst +++ b/doc/release/0.2.rst @@ -61,6 +61,7 @@ Bug fixes #. Fixed bug in ``armi/reactor/blocks.py::HexBlock::rotatePins`` that failed to modify pinLocation parameter (`#855 `_). #. Fixed bug in ``armi/reactor.py::Core::_applyThermalExpansion`` that failed to call the ``block.completeInitiaLoading`` (`#885 `_). #. Fixed bug where a settings validator would complain that both simple and detailed cycles inputs were being used even if only the detailed input was entered. +#. Fixed a bug where `getReactionRates()` was not accounting for burnup-dependent cross-sections. This changes the signature for the `getReactionRateDict()` function. #. TBD ARMI v0.2.3