Repository: terrapower/armi
Pull request #1412 (2023-09-14 22:56:29)

Representative block component temperature fix
## What is the change?

Number densities are set on a representative block to match the average number densities of all blocks in the `BlockCollection`. `block.setNumberDensities()` will smear the mass of a nuclide evenly over all components that contain the nuclide.

This PR aims to avoid the smearing when possible by setting the number densities on a component-by-component basis. This is only possible if all of the blocks in the `BlockCollection` have the same number of components and types of each component.

## Why is the change being made?

If tracking a large number of fission products, fuel components can often have trace amounts of isotopes that are typically associated with cladding, structural, or coolant material. When `block.setNumberDensities()` is called in this situation, the structural and coolant mass are smeared into the fuel components. A `BlockCollection` uses `block.setNumberDensities()` to set the number density on the representative block for that collection.

When the lattice physics writer later queries the representative block to get the average temperature of a structural nuclide, e.g., Fe-56, it performs a mass-weighted temperature of all the Fe-56 in the block. Since Fe-56 is in the fuel component of the representative block, the fuel component temperature will contribute to the average nuclide temperature of Fe-56. This is artificial and undesirable.

The solution here is to prevent the masses from being incorrectly smeared amongst the components by averaging number densities at the component level instead of the block level. The solution will not work if the `BlockCollection` has blocks of different types with different components; in this case, the code will revert to the old behavior.

<!-- MANDATORY: Explain why the change is necessary -->
<!-- Optional: Link to any related GitHub Issues -->

This addresses one problem arising from the issue described in #1374, but this PR alone is not enough to close that issue. 

---

## Checklist

<!--
    You (the pull requester) should put an `x` in the boxes below you have completed.
    If you're unsure about any of them, don't hesitate to ask. We're here to help!
    Learn what a "good PR" looks like here:
    https://terrapower.github.io/armi/developer/tooling.html#good-pull-requests
-->

- [x] This PR has only one purpose or idea.
- [x] Tests have been added/updated to verify that the new/changed code works.

<!-- Check the code quality -->

- [x] The code style follows [good practices](https://terrapower.github.io/armi/developer/standards_and_practices.html).
- [x] The commit message(s) follow [good practices](https://terrapower.github.io/armi/developer/tooling.html).

<!-- Check the project-level cruft -->

- [x] The [release notes](https://terrapower.github.io/armi/release/index.html) (location `doc/release/0.X.rst`) are up-to-date with any important changes.
- [x] The documentation is still up-to-date in the `doc` folder.
- [x] The dependencies are still up-to-date in `setup.py`.


keckler: Is this the problem case that #1374 was hinting at?
john-science: @mgjarrett Have you tested to see what affect this has on our downstream modeling?  If you'd like, we can talk about that offline.
mgjarrett: > Is this the problem case that #1374 was hinting at?

This is one instance of it. I can think of at least one other (uniform mesh converter), but there could also be more.
mgjarrett: I updated one docstring and somehow all of the tests started failing. Is there some GH actions problem?
john-science: > It looks like matplotlib 3.8.0 was released in the last couple days...

Yeah, that fix is in ARMI `main`. It looks like your PR's CI is running again.
albeanth: > If tracking a large number of fission products, fuel components can often have trace amounts of isotopes that are typically associated with cladding, structural, or coolant material. 

I might get the answer to this during my review, but what causes this?
jakehader: @albeanth this is related to my implementation changes to the fission product model with explicit fission products. Basically if anything is depletable then it gets all trace nuclides since the depletion chain isn't generically represented in the framework and is offloaded to a code like ORIGEN. To fix this we can make some assumptions about what should go into fuel versus other components but the most comprehensive way would be to have a detailed decay and transmutation chain implemented in the framework that could be traversed to ensure fission products are not assigned to clad, coolant, and other structural components.
mgjarrett: @albeanth I added several unit tests and updated some of the docstrings and comments
john-science: @mgjarrett The CI has been red on this PR for a week or two, so I hadn't been looking at it.  (My eyes tend to skip over PRs where the CI is red, it's not on purpose.)

It turns out it was just a merge conflict on the release notes. Easy.
john-science: @mgjarrett Since this branch is in your forked repo, I can fixed the `black` error:

```bash
would reformat armi/physics/neutronics/tests/test_crossSectionManager.py
````

Patch:
diff --git a/armi/physics/neutronics/crossSectionGroupManager.py b/armi/physics/neutronics/crossSectionGroupManager.py
index 3c816217..11a7fb53 100644
--- a/armi/physics/neutronics/crossSectionGroupManager.py
+++ b/armi/physics/neutronics/crossSectionGroupManager.py
@@ -120,10 +120,13 @@ class BlockCollection(list):
     This is a list with special methods.
     """
 
-    def __init__(self, allNuclidesInProblem, validBlockTypes=None):
+    def __init__(
+        self, allNuclidesInProblem, validBlockTypes=None, averageByComponent=False
+    ):
         list.__init__(self)
         self.allNuclidesInProblem = allNuclidesInProblem
         self.weightingParam = None
+        self.averageByComponent = averageByComponent
 
         # allowed to be independent of fuel component temperatures b/c Doppler
         self.avgNucTemperatures = {}
@@ -319,7 +322,17 @@ class AverageBlockCollection(BlockCollection):
         newBlock = self._getNewBlock()
         lfpCollection = self._getAverageFuelLFP()
         newBlock.setLumpedFissionProducts(lfpCollection)
-        newBlock.setNumberDensities(self._getAverageNumberDensities())
+        # check if components are similar
+        if self._performAverageByComponent():
+            # set number densities and temperatures on a component basis
+            for compIndex, c in enumerate(sorted(newBlock.getComponents())):
+                c.setNumberDensities(
+                    self._getAverageComponentNumberDensities(compIndex)
+                )
+                c.temperatureInC = self._getAverageComponentTemperature(compIndex)
+        else:
+            newBlock.setNumberDensities(self._getAverageNumberDensities())
+
         newBlock.p.percentBu = self._calcWeightedBurnup()
         self.calcAvgNuclideTemperatures()
         return newBlock
@@ -359,6 +372,96 @@ class AverageBlockCollection(BlockCollection):
             nv += nvBlock * wt
         return nvt, nv
 
+    def _getAverageComponentNumberDensities(self, compIndex):
+        """
+        Get weighted average number densities of a component in the collection.
+
+        Returns
+        -------
+        numberDensities : dict
+            nucName, ndens data (atoms/bn-cm)
+        """
+        nuclides = self.allNuclidesInProblem
+        blocks = self.getCandidateBlocks()
+        weights = numpy.array([self.getWeight(b) for b in blocks])
+        weights /= weights.sum()  # normalize by total weight
+        components = [sorted(b.getComponents())[compIndex] for b in blocks]
+        ndens = weights.dot([c.getNuclideNumberDensities(nuclides) for c in components])
+        return dict(zip(nuclides, ndens))
+
+    def _getAverageComponentTemperature(self, compIndex):
+        """
+        Get weighted average component temperature for the collection
+
+        Notes
+        -----
+        Weighting is both by the block weight within the collection and the relative mass of the component.
+        The block weight is already scaled by the block volume, so we need to pull that out of the block
+        weighting because it would effectively be double-counted in the component mass. b.getHeight()
+        is proportional to block volume, so it is used here as a computationally cheaper proxy for scaling
+        by block volume.
+
+        Returns
+        -------
+        numberDensities : dict
+            nucName, ndens data (atoms/bn-cm)
+        """
+        blocks = self.getCandidateBlocks()
+        weights = numpy.array([self.getWeight(b) / b.getHeight() for b in blocks])
+        weights /= weights.sum()  # normalize by total weight
+        components = [sorted(b.getComponents())[compIndex] for b in blocks]
+        weightedAvgComponentMass = sum(
+            w * c.getMass() for w, c in zip(weights, components)
+        )
+        if weightedAvgComponentMass == 0.0:
+            # if there is no component mass (e.g., gap), do a regular average
+            return numpy.mean(numpy.array([c.temperatureInC for c in components]))
+        else:
+            return (
+                weights.dot(
+                    numpy.array([c.temperatureInC * c.getMass() for c in components])
+                )
+                / weightedAvgComponentMass
+            )
+
+    def _performAverageByComponent(self):
+        """
+        Check if block collection averaging can/should be performed by component
+
+        If the components of blocks in the collection are similar and the user
+        has requested component-level averaging, return True.
+        Otherwise, return False.
+        """
+        if not self.averageByComponent:
+            return False
+        else:
+            return self._checkBlockSimilarity()
+
+    def _checkBlockSimilarity(self):
+        """
+        Check if blocks in the collection have similar components
+
+        If the components of blocks in the collection are similar and the user
+        has requested component-level averaging, return True.
+        Otherwise, return False.
+        """
+        cFlags = dict()
+        for b in self.getCandidateBlocks():
+            cFlags[b] = [c.p.flags for c in sorted(b.getComponents())]
+        refB = b
+        refFlags = cFlags[refB]
+        for b, compFlags in cFlags.items():
+            for c, refC in zip(compFlags, refFlags):
+                if c != refC:
+                    runLog.warning(
+                        "Non-matching block in AverageBlockCollection!\n"
+                        f"{refC} component flags in {refB} does not match {c} in {b}.\n"
+                        f"Number densities will be smeared in representative block."
+                    )
+                    return False
+        else:
+            return True
+
 
 def getBlockNuclideTemperatureAvgTerms(block, allNucNames):
     """
@@ -370,22 +473,22 @@ def getBlockNuclideTemperatureAvgTerms(block, allNucNames):
     as trace values at the proper component temperatures.
     """
 
-    def getNumberDensityWithTrace(component, nucName):
-        # needed to make sure temperature of 0-density nuclides in fuel get fuel temperature
-        try:
-            dens = component.p.numberDensities[nucName] or TRACE_NUMBER_DENSITY
-        except KeyError:
-            dens = 0.0
-        return dens
+    def getNumberDensitiesWithTrace(component, allNucNames):
+        """
+        Needed to make sure temperature of 0-density nuclides in fuel get fuel temperature
+        """
+        return [
+            component.p.numberDensities[nucName] or TRACE_NUMBER_DENSITY
+            if nucName in component.p.numberDensities
+            else 0.0
+            for nucName in allNucNames
+        ]
 
     vol = block.getVolume()
     components, volFracs = zip(*block.getVolumeFractions())
     # D = CxN matrix of number densities
     ndens = numpy.array(
-        [
-            [getNumberDensityWithTrace(c, nucName) for nucName in allNucNames]
-            for c in components
-        ]
+        [getNumberDensitiesWithTrace(c, allNucNames) for c in components]
     )
     temperatures = numpy.array(
         [c.temperatureInC for c in components]
@@ -1314,6 +1417,9 @@ def blockCollectionFactory(xsSettings, allNuclidesInProblem):
     """Build a block collection based on user settings and input."""
     blockRepresentation = xsSettings.blockRepresentation
     validBlockTypes = xsSettings.validBlockTypes
+    averageByComponent = xsSettings.averageByComponent
     return BLOCK_COLLECTIONS[blockRepresentation](
-        allNuclidesInProblem, validBlockTypes=validBlockTypes
+        allNuclidesInProblem,
+        validBlockTypes=validBlockTypes,
+        averageByComponent=averageByComponent,
     )
diff --git a/armi/physics/neutronics/crossSectionSettings.py b/armi/physics/neutronics/crossSectionSettings.py
index 92c58019..5cdbdf55 100644
--- a/armi/physics/neutronics/crossSectionSettings.py
+++ b/armi/physics/neutronics/crossSectionSettings.py
@@ -54,6 +54,7 @@ CONF_REACTION_DRIVER = "nuclideReactionDriver"
 CONF_XSID = "xsID"
 CONF_XS_EXECUTE_EXCLUSIVE = "xsExecuteExclusive"
 CONF_XS_PRIORITY = "xsPriority"
+CONF_COMPONENT_AVERAGING = "averageByComponent"
 CONF_XS_MAX_ATOM_NUMBER = "xsMaxAtomNumber"
 CONF_MIN_DRIVER_DENSITY = "minDriverDensity"
 
@@ -113,6 +114,7 @@ _VALID_INPUTS_BY_GEOMETRY_TYPE = {
         CONF_BLOCKTYPES,
         CONF_BLOCK_REPRESENTATION,
         CONF_EXTERNAL_FLUX_FILE_LOCATION,
+        CONF_COMPONENT_AVERAGING,
         CONF_XS_EXECUTE_EXCLUSIVE,
         CONF_XS_PRIORITY,
         CONF_XS_MAX_ATOM_NUMBER,
@@ -124,6 +126,7 @@ _VALID_INPUTS_BY_GEOMETRY_TYPE = {
         CONF_BLOCKTYPES,
         CONF_BLOCK_REPRESENTATION,
         CONF_EXTERNAL_FLUX_FILE_LOCATION,
+        CONF_COMPONENT_AVERAGING,
         CONF_XS_EXECUTE_EXCLUSIVE,
         CONF_XS_PRIORITY,
         CONF_XS_MAX_ATOM_NUMBER,
@@ -141,6 +144,7 @@ _VALID_INPUTS_BY_GEOMETRY_TYPE = {
         CONF_BLOCKTYPES,
         CONF_BLOCK_REPRESENTATION,
         CONF_EXTERNAL_FLUX_FILE_LOCATION,
+        CONF_COMPONENT_AVERAGING,
         CONF_XS_EXECUTE_EXCLUSIVE,
         CONF_XS_PRIORITY,
         CONF_XS_MAX_ATOM_NUMBER,
@@ -156,6 +160,7 @@ _VALID_INPUTS_BY_GEOMETRY_TYPE = {
         CONF_EXTERNAL_RINGS,
         CONF_BLOCK_REPRESENTATION,
         CONF_EXTERNAL_FLUX_FILE_LOCATION,
+        CONF_COMPONENT_AVERAGING,
         CONF_XS_EXECUTE_EXCLUSIVE,
         CONF_XS_PRIORITY,
         CONF_XS_MAX_ATOM_NUMBER,
@@ -188,6 +193,7 @@ _SINGLE_XS_SCHEMA = vol.Schema(
         vol.Optional(CONF_XS_PRIORITY): vol.Coerce(float),
         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,
     }
 )
 
@@ -456,6 +462,7 @@ class XSModelingOptions:
         xsExecuteExclusive=None,
         xsPriority=None,
         xsMaxAtomNumber=None,
+        averageByComponent=False,
         minDriverDensity=0.0,
     ):
         self.xsID = xsID
@@ -478,6 +485,7 @@ class XSModelingOptions:
         self.meshSubdivisionsPerCm = meshSubdivisionsPerCm
         self.xsMaxAtomNumber = xsMaxAtomNumber
         self.minDriverDensity = minDriverDensity
+        self.averageByComponent = averageByComponent
         # these are related to execution
         self.xsExecuteExclusive = xsExecuteExclusive
         self.xsPriority = xsPriority
@@ -688,6 +696,7 @@ class XSModelingOptions:
 
         defaults[CONF_XS_EXECUTE_EXCLUSIVE] = False
         defaults[CONF_XS_PRIORITY] = 5
+        defaults[CONF_COMPONENT_AVERAGING] = False
 
         for attrName, defaultValue in defaults.items():
             currentValue = getattr(self, attrName)
diff --git a/doc/release/0.2.rst b/doc/release/0.2.rst
index 62bed0fc..10dabfde 100644
--- a/doc/release/0.2.rst
+++ b/doc/release/0.2.rst
@@ -12,6 +12,7 @@ What's new in ARMI
 #. ARMI now mandates ``ruff`` linting. (`PR#1419 <https://github.com/terrapower/armi/pull/1419>`_)
 #. Removed all old ARMI requirements, to start the work fresh. (`PR#1438 <https://github.com/terrapower/armi/pull/1438>`_)
 #. Downgrading Draft PRs as policy. (`PR#1444 <https://github.com/terrapower/armi/pull/1444>`_)
+#. Attempt to set representative block number densities by component if possible. (`PR#1412 <https://github.com/terrapower/armi/pull/1412>`_)
 #. TBD
 
 Bug fixes
@@ -308,4 +309,3 @@ Bug fixes
 #. Fixed order of operations issue in rotatePins
 #. Fixed incorrect multiplicity for non-grid block components
 #. Many additional bugfixes and cleanups (see PR list)
-


Test patch:
diff --git a/armi/physics/neutronics/tests/test_crossSectionManager.py b/armi/physics/neutronics/tests/test_crossSectionManager.py
index 2282651c..f6c1a040 100644
--- a/armi/physics/neutronics/tests/test_crossSectionManager.py
+++ b/armi/physics/neutronics/tests/test_crossSectionManager.py
@@ -21,6 +21,7 @@ from io import BytesIO
 import copy
 import os
 import unittest
+from unittest.mock import MagicMock
 
 from six.moves import cPickle
 
@@ -43,8 +44,9 @@ from armi.physics.neutronics.settings import (
 )
 from armi.reactor.blocks import HexBlock
 from armi.reactor.flags import Flags
-from armi.reactor.tests import test_reactors
+from armi.reactor.tests import test_reactors, test_blocks
 from armi.tests import TEST_ROOT
+from armi.tests import mockRunLogs
 from armi.utils import units
 from armi.utils.directoryChangers import TemporaryDirectoryChanger
 
@@ -97,28 +99,238 @@ class TestBlockCollectionMedian(unittest.TestCase):
 
 
 class TestBlockCollectionAverage(unittest.TestCase):
-    def setUp(self):
+    @classmethod
+    def setUpClass(cls):
         fpFactory = test_lumpedFissionProduct.getDummyLFPFile()
-        self.blockList = makeBlocks(5)
-        for bi, b in enumerate(self.blockList):
+        cls.blockList = makeBlocks(5)
+        for bi, b in enumerate(cls.blockList):
             b.setType("fuel")
             b.p.percentBu = bi / 4.0 * 100
             b.setLumpedFissionProducts(fpFactory.createLFPsFromFile())
+            # put some trace Fe-56 and Na-23 into the fuel
+            # zero out all fuel nuclides except U-235 (for mass-weighting of component temperature)
+            fuelComp = b.getComponent(Flags.FUEL)
+            for nuc in fuelComp.getNuclides():
+                b.setNumberDensity(nuc, 0.0)
             b.setNumberDensity("U235", bi)
+            fuelComp.setNumberDensity("FE56", 1e-15)
+            fuelComp.setNumberDensity("NA23", 1e-15)
             b.p.gasReleaseFraction = bi * 2 / 8.0
+            for c in b:
+                if c.hasFlags(Flags.FUEL):
+                    c.temperatureInC = 600.0 + bi
+                elif c.hasFlags([Flags.CLAD, Flags.DUCT, Flags.WIRE]):
+                    c.temperatureInC = 500.0 + bi
+                elif c.hasFlags([Flags.BOND, Flags.COOLANT, Flags.INTERCOOLANT]):
+                    c.temperatureInC = 400.0 + bi
+
+    def setUp(self):
         self.bc = AverageBlockCollection(
             self.blockList[0].r.blueprints.allNuclidesInProblem
         )
         self.bc.extend(self.blockList)
+        self.bc.averageByComponent = True
+
+    def test_performAverageByComponent(self):
+        """
+        Check the averageByComponent attribute
+        """
+        self.bc._checkBlockSimilarity = MagicMock(return_value=True)
+        self.assertTrue(self.bc._performAverageByComponent())
+        self.bc.averageByComponent = False
+        self.assertFalse(self.bc._performAverageByComponent())
+
+    def test_checkBlockSimilarity(self):
+        """
+        Check the block similarity test
+        """
+        self.assertTrue(self.bc._checkBlockSimilarity())
+        self.bc.append(test_blocks.loadTestBlock())
+        self.assertFalse(self.bc._checkBlockSimilarity())
 
     def test_createRepresentativeBlock(self):
+        """
+        Test creation of a representative block
+        """
         avgB = self.bc.createRepresentativeBlock()
         self.assertNotIn(avgB, self.bc)
-        # 0 + 1/4 + 2/4 + 3/4 + 4/4 =
-        # (0 + 1 + 2 + 3 + 4 ) / 5 = 10/5 = 2.0
+        # (0 + 1 + 2 + 3 + 4) / 5 = 10/5 = 2.0
         self.assertAlmostEqual(avgB.getNumberDensity("U235"), 2.0)
+        # (0 + 1/4 + 2/4 + 3/4 + 4/4) / 5 * 100.0 = 50.0
         self.assertEqual(avgB.p.percentBu, 50.0)
 
+        # check that a new block collection of the representative block has right temperatures
+        # this is required for Doppler coefficient calculations
+        newBc = AverageBlockCollection(
+            self.blockList[0].r.blueprints.allNuclidesInProblem
+        )
+        newBc.append(avgB)
+        newBc.calcAvgNuclideTemperatures()
+        self.assertAlmostEqual(newBc.avgNucTemperatures["U235"], 603.0)
+        self.assertAlmostEqual(newBc.avgNucTemperatures["FE56"], 502.0)
+        self.assertAlmostEqual(newBc.avgNucTemperatures["NA23"], 402.0)
+
+    def test_createRepresentativeBlockDissimilar(self):
+        """
+        Test creation of a representative block from a collection with dissimilar blocks
+        """
+
+        uniqueBlock = test_blocks.loadTestBlock()
+        uniqueBlock.p.percentBu = 50.0
+        fpFactory = test_lumpedFissionProduct.getDummyLFPFile()
+        uniqueBlock.setLumpedFissionProducts(fpFactory.createLFPsFromFile())
+        uniqueBlock.setNumberDensity("U235", 2.0)
+        uniqueBlock.p.gasReleaseFraction = 1.0
+        for c in uniqueBlock:
+            if c.hasFlags(Flags.FUEL):
+                c.temperatureInC = 600.0
+            elif c.hasFlags([Flags.CLAD, Flags.DUCT, Flags.WIRE]):
+                c.temperatureInC = 500.0
+            elif c.hasFlags([Flags.BOND, Flags.COOLANT, Flags.INTERCOOLANT]):
+                c.temperatureInC = 400.0
+        self.bc.append(uniqueBlock)
+
+        with mockRunLogs.BufferLog() as mock:
+            avgB = self.bc.createRepresentativeBlock()
+            self.assertIn(
+                "Non-matching block in AverageBlockCollection", mock.getStdout()
+            )
+
+        self.assertNotIn(avgB, self.bc)
+        # (0 + 1 + 2 + 3 + 4 + 2) / 6.0 = 12/6 = 2.0
+        self.assertAlmostEqual(avgB.getNumberDensity("U235"), 2.0)
+        # (0 + 1/4 + 2/4 + 3/4 + 4/4) / 5 * 100.0 = 50.0
+        self.assertAlmostEqual(avgB.p.percentBu, 50.0)
+
+        # U35 has different average temperature because blocks have different U235 content
+        newBc = AverageBlockCollection(
+            self.blockList[0].r.blueprints.allNuclidesInProblem
+        )
+        newBc.append(avgB)
+        newBc.calcAvgNuclideTemperatures()
+        # temps expected to be proportional to volume-fraction weighted temperature
+        # this is a non-physical result, but it demonstrates a problem that exists in the code
+        # when dissimilar blocks are put together in a BlockCollection
+        structureVolume = sum(
+            c.getVolume()
+            for c in avgB.getComponents([Flags.CLAD, Flags.DUCT, Flags.WIRE])
+        )
+        fuelVolume = avgB.getComponent(Flags.FUEL).getVolume()
+        coolantVolume = sum(
+            c.getVolume()
+            for c in avgB.getComponents([Flags.BOND, Flags.COOLANT, Flags.INTERCOOLANT])
+        )
+        expectedIronTemp = (structureVolume * 500.0 + fuelVolume * 600.0) / (
+            structureVolume + fuelVolume
+        )
+        expectedSodiumTemp = (coolantVolume * 400.0 + fuelVolume * 600.0) / (
+            coolantVolume + fuelVolume
+        )
+        self.assertAlmostEqual(newBc.avgNucTemperatures["U235"], 600.0)
+        self.assertAlmostEqual(newBc.avgNucTemperatures["FE56"], expectedIronTemp)
+        self.assertAlmostEqual(newBc.avgNucTemperatures["NA23"], expectedSodiumTemp)
+
+
+class TestComponentAveraging(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        fpFactory = test_lumpedFissionProduct.getDummyLFPFile()
+        cls.blockList = makeBlocks(3)
+        for bi, b in enumerate(cls.blockList):
+            b.setType("fuel")
+            b.setLumpedFissionProducts(fpFactory.createLFPsFromFile())
+            # put some trace Fe-56 and Na-23 into the fuel
+            # zero out all fuel nuclides except U-235 (for mass-weighting of component temperature)
+            for nuc in b.getNuclides():
+                b.setNumberDensity(nuc, 0.0)
+            b.setNumberDensity("U235", bi)
+            b.setNumberDensity("FE56", bi / 2.0)
+            b.setNumberDensity("NA23", bi / 3.0)
+            for c in b:
+                if c.hasFlags(Flags.FUEL):
+                    c.temperatureInC = 600.0 + bi
+                elif c.hasFlags([Flags.CLAD, Flags.DUCT, Flags.WIRE]):
+                    c.temperatureInC = 500.0 + bi
+                elif c.hasFlags([Flags.BOND, Flags.COOLANT, Flags.INTERCOOLANT]):
+                    c.temperatureInC = 400.0 + bi
+
+    def setUp(self):
+        self.bc = AverageBlockCollection(
+            self.blockList[0].r.blueprints.allNuclidesInProblem
+        )
+        blockCopies = [copy.deepcopy(b) for b in self.blockList]
+        self.bc.extend(blockCopies)
+
+    def test_getAverageComponentNumberDensities(self):
+        """
+        Test component number density averaging
+        """
+        # becaue of the way densities are set up, the middle block (index 1 of 0-2) component densities are equivalent to the average
+        b = self.bc[1]
+        for compIndex, c in enumerate(b.getComponents()):
+            avgDensities = self.bc._getAverageComponentNumberDensities(compIndex)
+            compDensities = c.getNumberDensities()
+            for nuc in c.getNuclides():
+                self.assertAlmostEqual(
+                    compDensities[nuc],
+                    avgDensities[nuc],
+                    msg=f"{nuc} density {compDensities[nuc]} not equal to {avgDensities[nuc]}!",
+                )
+
+    def test_getAverageComponentTemperature(self):
+        """
+        Test mass-weighted component temperature averaging
+        """
+        b = self.bc[0]
+        massWeightedIncrease = 5.0 / 3.0
+        baseTemps = [600, 400, 500, 500, 400, 500, 400]
+        expectedTemps = [t + massWeightedIncrease for t in baseTemps]
+        for compIndex, c in enumerate(b.getComponents()):
+            avgTemp = self.bc._getAverageComponentTemperature(compIndex)
+            self.assertAlmostEqual(
+                expectedTemps[compIndex],
+                avgTemp,
+                msg=f"{c} avg temperature {avgTemp} not equal to expected {expectedTemps[compIndex]}!",
+            )
+
+    def test_getAverageComponentTemperatureVariedWeights(self):
+        """
+        Test mass-weighted component temperature averaging with variable weights
+        """
+        # make up a fake weighting with power param
+        self.bc.weightingParam = "power"
+        for i, b in enumerate(self.bc):
+            b.p.power = i
+        weightedIncrease = 1.8
+        baseTemps = [600, 400, 500, 500, 400, 500, 400]
+        expectedTemps = [t + weightedIncrease for t in baseTemps]
+        for compIndex, c in enumerate(b.getComponents()):
+            avgTemp = self.bc._getAverageComponentTemperature(compIndex)
+            self.assertAlmostEqual(
+                expectedTemps[compIndex],
+                avgTemp,
+                msg=f"{c} avg temperature {avgTemp} not equal to expected {expectedTemps[compIndex]}!",
+            )
+
+    def test_getAverageComponentTemperatureNoMass(self):
+        """
+        Test component temperature averaging when the components have no mass
+        """
+        for b in self.bc:
+            for nuc in b.getNuclides():
+                b.setNumberDensity(nuc, 0.0)
+
+        unweightedIncrease = 1.0
+        baseTemps = [600, 400, 500, 500, 400, 500, 400]
+        expectedTemps = [t + unweightedIncrease for t in baseTemps]
+        for compIndex, c in enumerate(b.getComponents()):
+            avgTemp = self.bc._getAverageComponentTemperature(compIndex)
+            self.assertAlmostEqual(
+                expectedTemps[compIndex],
+                avgTemp,
+                msg=f"{c} avg temperature {avgTemp} not equal to expected {expectedTemps[compIndex]}!",
+            )
+
 
 class TestBlockCollectionComponentAverage(unittest.TestCase):
     r"""tests for ZPPR 1D XS gen cases."""
@@ -442,10 +654,11 @@ class TestBlockCollectionComponentAverage1DCylinder(unittest.TestCase):
 
 
 class TestBlockCollectionFluxWeightedAverage(unittest.TestCase):
-    def setUp(self):
+    @classmethod
+    def setUpClass(cls):
         fpFactory = test_lumpedFissionProduct.getDummyLFPFile()
-        self.blockList = makeBlocks(5)
-        for bi, b in enumerate(self.blockList):
+        cls.blockList = makeBlocks(5)
+        for bi, b in enumerate(cls.blockList):
             b.setType("fuel")
             b.p.percentBu = bi / 4.0 * 100
             b.setLumpedFissionProducts(fpFactory.createLFPsFromFile())
@@ -453,6 +666,7 @@ class TestBlockCollectionFluxWeightedAverage(unittest.TestCase):
             b.p.gasReleaseFraction = bi * 2 / 8.0
             b.p.flux = bi + 1
 
+    def setUp(self):
         self.bc = FluxWeightedAverageBlockCollection(
             self.blockList[0].r.blueprints.allNuclidesInProblem
         )