Allowing creation of partially heterogeneous assemblies for 1D XS mod…

…el (#1949)

armi/nuclearDataIO/nuclearFileMetadata.py

@@ -252,4 +252,4 @@ def _getSkippedKeys(self, other, selfContainer, otherContainer, mergedData):
 
 
 class NuclideMetadata(_Metadata):
-    """Simple dictionary for providing metadata about how to read/write a nuclde to/from a file."""
+    """Simple dictionary for providing metadata about how to read/write a nuclide to/from a file."""

armi/nuclearDataIO/tests/test_xsLibraries.py

@@ -127,7 +127,7 @@ def test_mergeFailsWithNonIsotxsFiles(self):
             os.remove(dummyFileName)
 
         with TemporaryDirectoryChanger():
-            dummyFileName = "ISOtopics.txt"
+            dummyFileName = "ISO[]"
             with open(dummyFileName, "w") as file:
                 file.write(
                     "This is a file that starts with the letters 'ISO' but will"

armi/nuclearDataIO/xsLibraries.py

@@ -127,7 +127,7 @@ def getISOTXSLibrariesToMerge(xsLibrarySuffix, xsLibFileNames):
         isosWithSuffix = [
             iso
             for iso in isosToMerge
-            if re.match(f".*ISO[A-Z]{{2}}F?{xsLibrarySuffix}$", iso)
+            if re.match(f".*ISO[A-Za-z]{{2}}F?{xsLibrarySuffix}$", iso)
         ]
         isosToMerge = [
             iso
@@ -193,7 +193,7 @@ def mergeXSLibrariesInWorkingDirectory(
     for xsLibFilePath in sorted(xsLibFiles):
         try:
             # get XS ID from the cross section library name
-            xsID = re.search("ISO([A-Z0-9]{2})", xsLibFilePath).group(1)
+            xsID = re.search("ISO([A-Z0-9a-z]{2})", xsLibFilePath).group(1)
         except AttributeError:
             # if glob has matched something that is not actually an ISOXX file,
             # the .group() call will fail

armi/physics/neutronics/crossSectionGroupManager.py

@@ -63,6 +63,7 @@
 from armi.physics.neutronics.const import CONF_CROSS_SECTION
 from armi.reactor import flags
 from armi.reactor.components import basicShapes
+from armi.reactor.converters.blockConverters import stripComponents
 from armi.reactor.flags import Flags
 from armi.utils.units import TRACE_NUMBER_DENSITY
 
@@ -623,7 +624,11 @@ def _getAverageComponentNucs(self, components, bWeights):
             weight = bWeight * c.getArea()
             totalWeight += weight
             densities += weight * np.array(c.getNuclideNumberDensities(allNucNames))
-        return allNucNames, densities / totalWeight
+        if totalWeight > 0.0:
+            weightedDensities = densities / totalWeight
+        else:
+            weightedDensities = np.zeros_like(densities)
+        return allNucNames, weightedDensities
 
     def _orderComponentsInGroup(self, repBlock):
         """Order the components based on dimension and material type within the representative
@@ -648,6 +653,64 @@ def _getNucTempHelper(self):
         return nvt, nv
 
 
+class CylindricalComponentsDuctHetAverageBlockCollection(
+    CylindricalComponentsAverageBlockCollection
+):
+    """
+    Creates a representative block for the purpose of cross section generation with a one-
+    dimensional cylindrical model where all material inside the duct is homogenized.
+
+    .. impl:: Create partially heterogeneous representative blocks.
+        :id: I_ARMI_XSGM_CREATE_REPR_BLOCKS2
+        :implements: R_ARMI_XSGM_CREATE_REPR_BLOCKS
+
+        This class constructs representative blocks based on a volume-weighted average using
+        cylindrical blocks from an existing block list. Inheriting functionality from the abstract
+        :py:class:`Reactor <armi.physics.neutronics.crossSectionGroupManager.BlockCollection>`
+        object, this class will construct representative blocks using averaged parameters of all
+        blocks in the given collection. Number density averages are computed at a component level.
+        Nuclide temperatures from a median block-average temperature are used and the average burnup
+        is evaluated across all blocks in the block list.
+
+        The average nuclide temperatures are calculated only for the homogenized region inside of
+        the duct. For the non-homogenized regions, the MC2 writer uses the component temperatures.
+
+    Notes
+    -----
+    The representative block for this collection is the same as the parent. The only difference between
+    the two collection types is that this collection calculates average nuclide temperatures based only
+    on the components that are inside of the duct.
+    """
+
+    def _getNewBlock(self):
+        newBlock = copy.deepcopy(self._selectCandidateBlock())
+        newBlock.name = "1D_CYL_DUCT_HET_AVG_" + newBlock.getMicroSuffix()
+        return newBlock
+
+    def _makeRepresentativeBlock(self):
+        """Build a representative fuel block based on component number densities."""
+        self.calcAvgNuclideTemperatures()
+        return CylindricalComponentsAverageBlockCollection._makeRepresentativeBlock(
+            self
+        )
+
+    def _getNucTempHelper(self):
+        """All candidate blocks are used in the average."""
+        nvt = np.zeros(len(self.allNuclidesInProblem))
+        nv = np.zeros(len(self.allNuclidesInProblem))
+        for block in self.getCandidateBlocks():
+            wt = self.getWeight(block)
+            # remove the duct and intercoolant from the block before
+            # calculating average nuclide temps
+            newBlock, _mixtureFlags = stripComponents(block, Flags.DUCT)
+            nvtBlock, nvBlock = getBlockNuclideTemperatureAvgTerms(
+                newBlock, self.allNuclidesInProblem
+            )
+            nvt += nvtBlock * wt
+            nv += nvBlock * wt
+        return nvt, nv
+
+
 class SlabComponentsAverageBlockCollection(BlockCollection):
     """
     Creates a representative 1D slab block.
@@ -789,7 +852,11 @@ def _getAverageComponentNucs(self, components, bWeights):
             weight = bWeight * c.getArea()
             totalWeight += weight
             densities += weight * np.array(c.getNuclideNumberDensities(allNucNames))
-        return allNucNames, densities / totalWeight
+        if totalWeight > 0.0:
+            weightedDensities = densities / totalWeight
+        else:
+            weightedDensities = np.zeros_like(densities)
+        return allNucNames, weightedDensities
 
     def _orderComponentsInGroup(self, repBlock):
         """Order the components based on dimension and material type within the representative block."""
@@ -1520,6 +1587,9 @@ def updateNuclideTemperatures(self, blockCollectionByXsGroup=None):
 FLUX_WEIGHTED_AVERAGE_BLOCK_COLLECTION = "FluxWeightedAverage"
 SLAB_COMPONENTS_BLOCK_COLLECTION = "ComponentAverage1DSlab"
 CYLINDRICAL_COMPONENTS_BLOCK_COLLECTION = "ComponentAverage1DCylinder"
+CYLINDRICAL_COMPONENTS_DUCT_HET_BLOCK_COLLECTION = (
+    "ComponentAverage1DCylinderDuctHeterogeneous"
+)
 
 # Mapping between block collection string constants and their
 # respective block collection classes.
@@ -1529,12 +1599,17 @@ def updateNuclideTemperatures(self, blockCollectionByXsGroup=None):
     FLUX_WEIGHTED_AVERAGE_BLOCK_COLLECTION: FluxWeightedAverageBlockCollection,
     SLAB_COMPONENTS_BLOCK_COLLECTION: SlabComponentsAverageBlockCollection,
     CYLINDRICAL_COMPONENTS_BLOCK_COLLECTION: CylindricalComponentsAverageBlockCollection,
+    CYLINDRICAL_COMPONENTS_DUCT_HET_BLOCK_COLLECTION: CylindricalComponentsDuctHetAverageBlockCollection,
 }
 
 
 def blockCollectionFactory(xsSettings, allNuclidesInProblem):
     """Build a block collection based on user settings and input."""
     blockRepresentation = xsSettings.blockRepresentation
+    if (
+        blockRepresentation == CYLINDRICAL_COMPONENTS_BLOCK_COLLECTION
+    ) and xsSettings.ductHeterogeneous:
+        blockRepresentation = CYLINDRICAL_COMPONENTS_DUCT_HET_BLOCK_COLLECTION
     validBlockTypes = xsSettings.validBlockTypes
     averageByComponent = xsSettings.averageByComponent
     return BLOCK_COLLECTIONS[blockRepresentation](

armi/physics/neutronics/crossSectionSettings.py

@@ -49,6 +49,7 @@
 CONF_HOMOGBLOCK = "useHomogenizedBlockComposition"
 CONF_INTERNAL_RINGS = "numInternalRings"
 CONF_MERGE_INTO_CLAD = "mergeIntoClad"
+CONF_MERGE_INTO_FUEL = "mergeIntoFuel"
 CONF_MESH_PER_CM = "meshSubdivisionsPerCm"
 CONF_REACTION_DRIVER = "nuclideReactionDriver"
 CONF_XSID = "xsID"
@@ -57,6 +58,8 @@
 CONF_COMPONENT_AVERAGING = "averageByComponent"
 CONF_XS_MAX_ATOM_NUMBER = "xsMaxAtomNumber"
 CONF_MIN_DRIVER_DENSITY = "minDriverDensity"
+CONF_DUCT_HETEROGENEOUS = "ductHeterogeneous"
+CONF_TRACE_ISOTOPE_THRESHOLD = "traceIsotopeThreshold"
 
 
 class XSGeometryTypes(Enum):
@@ -136,6 +139,7 @@ def getStr(cls, typeSpec: Enum):
         CONF_XSID,
         CONF_GEOM,
         CONF_MERGE_INTO_CLAD,
+        CONF_MERGE_INTO_FUEL,
         CONF_DRIVER,
         CONF_HOMOGBLOCK,
         CONF_INTERNAL_RINGS,
@@ -149,6 +153,8 @@ def getStr(cls, typeSpec: Enum):
         CONF_XS_PRIORITY,
         CONF_XS_MAX_ATOM_NUMBER,
         CONF_MIN_DRIVER_DENSITY,
+        CONF_DUCT_HETEROGENEOUS,
+        CONF_TRACE_ISOTOPE_THRESHOLD,
     },
     XSGeometryTypes.getStr(XSGeometryTypes.TWO_DIMENSIONAL_HEX): {
         CONF_XSID,
@@ -186,6 +192,7 @@ def getStr(cls, typeSpec: Enum):
         vol.Optional(CONF_INTERNAL_RINGS): vol.Coerce(int),
         vol.Optional(CONF_EXTERNAL_RINGS): vol.Coerce(int),
         vol.Optional(CONF_MERGE_INTO_CLAD): [str],
+        vol.Optional(CONF_MERGE_INTO_FUEL): [str],
         vol.Optional(CONF_XS_FILE_LOCATION): [str],
         vol.Optional(CONF_EXTERNAL_FLUX_FILE_LOCATION): str,
         vol.Optional(CONF_MESH_PER_CM): vol.Coerce(float),
@@ -194,6 +201,8 @@ def getStr(cls, typeSpec: Enum):
         vol.Optional(CONF_XS_MAX_ATOM_NUMBER): vol.Coerce(int),
         vol.Optional(CONF_MIN_DRIVER_DENSITY): vol.Coerce(float),
         vol.Optional(CONF_COMPONENT_AVERAGING): bool,
+        vol.Optional(CONF_DUCT_HETEROGENEOUS): bool,
+        vol.Optional(CONF_TRACE_ISOTOPE_THRESHOLD): vol.Coerce(float),
     }
 )
 
@@ -410,6 +419,12 @@ class XSModelingOptions:
         and is sometimes used to merge a "gap" or low-density region into
         a "clad" region to avoid numerical issues.
 
+    mergeIntoFuel : list of str
+        This is a lattice physics configuration option that is a list of component
+        names to merge into a "fuel" component. This is highly-design specific
+        and is sometimes used to merge a "gap" or low-density region into
+        a "fuel" region to avoid numerical issues.
+
     meshSubdivisionsPerCm : float
         This is a lattice physics configuration option that can be used to control
         subregion meshing of the representative block in 1D problems.
@@ -419,7 +434,7 @@ class XSModelingOptions:
         no others will allocate to it. This is useful for time balancing when you
         have one task that takes much longer than the others.
 
-    xsPriority:
+    xsPriority: int
         The priority of the mpi tasks that results from this xsID. Lower priority
         will execute first. starting longer jobs first is generally more efficient.
 
@@ -429,10 +444,31 @@ class XSModelingOptions:
         (e.g., fission products) as a depletion product of an isotope with a much
         smaller atomic number.
 
+    averageByComponent: bool
+        Controls whether the representative block averaging is performed on a
+        component-by-component basis or on the block as a whole. If True, the
+        resulting representative block will have component compositions that
+        largely reflect those of the underlying blocks in the collection. If
+        False, the number densities of some nuclides in the individual
+        components may not be reflective of those of the underlying components
+        due to the block number density "dehomogenization".
+
     minDriverDensity: float
         The minimum number density for nuclides included in driver material for a 1D
         lattice physics model.
 
+    ductHeterogeneous : bool
+        This is a lattice physics configuration option used to enable a partially
+        heterogeneous approximation for a 1D cylindrical model. Everything inside of the
+        duct will be treated as homogeneous.
+
+    traceIsotopeThreshold : float
+        This is a lattice physics configuration option used to enable a separate 0D fuel
+        cross section calculation for trace fission products when using a 1D cross section
+        model. This can significantly reduce the memory and run time required for the 1D
+        model. The setting takes a float value that represents the number density cutoff
+        for isotopes to be considered "trace". If no value is provided, the default is 0.0.
+
     Notes
     -----
     Not all default attributes may be useful for your specific application and you may
@@ -458,12 +494,15 @@ def __init__(
         numInternalRings=None,
         numExternalRings=None,
         mergeIntoClad=None,
+        mergeIntoFuel=None,
         meshSubdivisionsPerCm=None,
         xsExecuteExclusive=None,
         xsPriority=None,
         xsMaxAtomNumber=None,
         averageByComponent=False,
         minDriverDensity=0.0,
+        ductHeterogeneous=False,
+        traceIsotopeThreshold=0.0,
     ):
         self.xsID = xsID
         self.geometry = geometry
@@ -482,10 +521,13 @@ def __init__(
         self.numInternalRings = numInternalRings
         self.numExternalRings = numExternalRings
         self.mergeIntoClad = mergeIntoClad
+        self.mergeIntoFuel = mergeIntoFuel
         self.meshSubdivisionsPerCm = meshSubdivisionsPerCm
         self.xsMaxAtomNumber = xsMaxAtomNumber
         self.minDriverDensity = minDriverDensity
         self.averageByComponent = averageByComponent
+        self.ductHeterogeneous = ductHeterogeneous
+        self.traceIsotopeThreshold = traceIsotopeThreshold
         # these are related to execution
         self.xsExecuteExclusive = xsExecuteExclusive
         self.xsPriority = xsPriority
@@ -668,12 +710,15 @@ def setDefaults(self, blockRepresentation, validBlockTypes):
                 CONF_GEOM: self.geometry,
                 CONF_DRIVER: "",
                 CONF_MERGE_INTO_CLAD: ["gap"],
+                CONF_MERGE_INTO_FUEL: [],
                 CONF_MESH_PER_CM: 1.0,
                 CONF_INTERNAL_RINGS: 0,
                 CONF_EXTERNAL_RINGS: 1,
                 CONF_HOMOGBLOCK: False,
                 CONF_BLOCK_REPRESENTATION: crossSectionGroupManager.CYLINDRICAL_COMPONENTS_BLOCK_COLLECTION,
                 CONF_BLOCKTYPES: validBlockTypes,
+                CONF_DUCT_HETEROGENEOUS: False,
+                CONF_TRACE_ISOTOPE_THRESHOLD: 0.0,
             }
         elif self.geometry == XSGeometryTypes.getStr(
             XSGeometryTypes.TWO_DIMENSIONAL_HEX

armi/physics/neutronics/latticePhysics/latticePhysicsWriter.py

@@ -110,9 +110,12 @@ def __init__(
         self.xsId = representativeBlock.getMicroSuffix()
         self.xsSettings = self.cs[CONF_CROSS_SECTION][self.xsId]
         self.mergeIntoClad = self.xsSettings.mergeIntoClad
+        self.mergeIntoFuel = self.xsSettings.mergeIntoFuel
         self.driverXsID = self.xsSettings.driverID
         self.numExternalRings = self.xsSettings.numExternalRings
         self.criticalBucklingSearchActive = self.xsSettings.criticalBuckling
+        self.ductHeterogeneous = self.xsSettings.ductHeterogeneous
+        self.traceIsotopeThreshold = self.xsSettings.traceIsotopeThreshold
 
         self.executeExclusive = self.xsSettings.xsExecuteExclusive
         self.priority = self.xsSettings.xsPriority
@@ -266,8 +269,12 @@ def _getAllNuclidesByCategory(self, component=None):
                 continue  # skip LFPs here but add individual FPs below.
 
             if isinstance(subjectObject, components.Component):
-                # Heterogeneous number densities and temperatures
-                nucTemperatureInC = subjectObject.temperatureInC
+                if self.ductHeterogeneous and "Homogenized" in subjectObject.name:
+                    # Nuclide temperatures representing heterogeneous model component temperatures
+                    nucTemperatureInC = self._getAvgNuclideTemperatureInC(nucName)
+                else:
+                    # Heterogeneous number densities and temperatures
+                    nucTemperatureInC = subjectObject.temperatureInC
             else:
                 # Homogeneous number densities and temperatures
                 nucTemperatureInC = self._getAvgNuclideTemperatureInC(nucName)

armi/physics/neutronics/tests/test_crossSectionSettings.py

@@ -118,6 +118,8 @@ def test_homogeneousXsDefaultSettingAssignment(self):
         self.assertNotIn("YA", xsModel)
         self.assertEqual(xsModel["YA"].geometry, "0D")
         self.assertEqual(xsModel["YA"].criticalBuckling, True)
+        self.assertEqual(xsModel["YA"].ductHeterogeneous, False)
+        self.assertEqual(xsModel["YA"].traceIsotopeThreshold, 0.0)
 
     def test_setDefaultSettingsByLowestBuGroupHomogeneous(self):
         # Initialize some micro suffix in the cross sections
@@ -177,14 +179,23 @@ def test_optionalKey(self):
         """Test that optional key shows up with default value."""
         cs = settings.Settings()
         xsModel = XSSettings()
-        da = XSModelingOptions("DA", geometry="1D cylinder", meshSubdivisionsPerCm=1.0)
+        da = XSModelingOptions(
+            "DA",
+            geometry="1D cylinder",
+            meshSubdivisionsPerCm=1.0,
+            ductHeterogeneous=True,
+            traceIsotopeThreshold=1.0e-5,
+        )
         xsModel["DA"] = da
         xsModel.setDefaults(
             cs[CONF_XS_BLOCK_REPRESENTATION],
             cs[CONF_DISABLE_BLOCK_TYPE_EXCLUSION_IN_XS_GENERATION],
         )
         self.assertEqual(xsModel["DA"].mergeIntoClad, ["gap"])
         self.assertEqual(xsModel["DA"].meshSubdivisionsPerCm, 1.0)
+        self.assertEqual(xsModel["DA"].ductHeterogeneous, True)
+        self.assertEqual(xsModel["DA"].traceIsotopeThreshold, 1.0e-5)
+        self.assertEqual(xsModel["DA"].mergeIntoFuel, [])
 
     def test_badCrossSections(self):
         with self.assertRaises(TypeError):