Repository: terrapower/armi
Pull request #2109 (2025-04-02 23:20:26)

Store number densities in numpy arrays instead of dictionary
## What is the change? Why is it being made?

<!-- MANDATORY: Describe the change -->

Store component number densities in two numpy arrays instead of a dictionary of `[str: float]`.

`c.p.nuclides`: array of byte strings (`dtype='S6'`) containing nuclide names
`c.p.numberDensities`: array of floats (`dtype='np.float64'`) containing the number densities for the corresponding nuclides in `c.p.nuclides`

There is a tradeoff here between speed (dictionary lookup is very fast) and memory footprint (dictionary is approximately 10x-20x larger in memory than two numpy arrays). This tradeoff is being made because some full-size reactor model objects constructed in downstream applications have become too large, causing runs to either crash or consume a disproportionate amount of computational resources.

## SCR Information

<!-- MANDATORY: uncomment one-and-only-one of these -->
Change Type: features
<!-- Change Type: fixes -->
<!-- Change Type: trivial -->
<!-- Change Type: docs -->

<!-- MANDATORY: Describe the change in one sentence -->
One-Sentence Description: Change component number density storage from dictionary to numpy arrays.

<!-- MANDATORY: Describe any impact on the requirements, all on one line -->
One-line Impact on Requirements: NA


---

## Checklist

<!--
    The pull request author should check the box if the condition is met OR if it does not apply.
-->

- [x] This PR has only [one purpose or idea](https://terrapower.github.io/armi/developer/tooling.html#one-idea-one-pr).
- [x] [Tests](https://terrapower.github.io/armi/developer/tooling.html#test-it) have been added/updated to verify any new/changed code.
- [x] The [documentation](https://terrapower.github.io/armi/developer/tooling.html#document-it) is still up-to-date in the `doc` folder.
- [x] The code style follows [good practices](https://terrapower.github.io/armi/developer/standards_and_practices.html).
- [x] The dependencies are still up-to-date in `pyproject.toml`.


john-science: Bump!

This is still interesting to me.

Patch:
diff --git a/armi/bookkeeping/db/compareDB3.py b/armi/bookkeeping/db/compareDB3.py
index c963b9f6..48daf928 100644
--- a/armi/bookkeeping/db/compareDB3.py
+++ b/armi/bookkeeping/db/compareDB3.py
@@ -360,6 +360,14 @@ def _diffSpecialData(
                 diffResults.addDiff(compName, paramName, np.inf, np.inf, np.inf)
                 return
 
+            if dSrc.dtype.type == np.bytes_ or dRef.dtype.type == np.bytes_:
+                # data is byte strings; can't be diffed like numbers
+                if np.array_equal(dSrc, dRef):
+                    diffResults.addDiff(name, name, 0.0, 0.0, 0.0)
+                else:
+                    diffResults.addDiff(name, name, np.inf, np.inf, np.inf)
+                return
+
             # Make sure not to try to compare empty arrays. Numpy is mediocre at these;
             # they are super degenerate and cannot participate in concatenation.
             if 0 not in dSrc.shape:
diff --git a/armi/bookkeeping/db/database.py b/armi/bookkeeping/db/database.py
index 4355a875..90aeff2b 100644
--- a/armi/bookkeeping/db/database.py
+++ b/armi/bookkeeping/db/database.py
@@ -1021,7 +1021,8 @@ class Database:
 
                 data = np.array(pDef.serializer.unpack(data, dataSet.attrs[_SERIALIZER_VERSION], attrs))
 
-            if data.dtype.type is np.bytes_:
+            # nuclides are a special case where we want to keep in np.bytes_ format
+            if data.dtype.type is np.bytes_ and paramName != "nuclides":
                 data = np.char.decode(data)
 
             if attrs.get("specialFormatting", False):
@@ -1044,18 +1045,21 @@ class Database:
                 runLog.error(msg)
                 raise ValueError(msg)
 
-            # iterating of np is not fast...
-            for c, val, linkedDim in itertools.zip_longest(comps, unpackedData, linkedDims, fillvalue=""):
-                try:
-                    if linkedDim != "":
-                        c.p[paramName] = linkedDim
-                    else:
-                        c.p[paramName] = val
-                except AssertionError as ae:
-                    # happens when a param was deprecated but being loaded from old DB
-                    runLog.warning(
-                        f"{str(ae)}\nSkipping load of invalid param `{paramName}` (possibly loading from old DB)\n"
-                    )
+            if paramName == "numberDensities" and attrs.get("dict", False):
+                Database._applyComponentNumberDensitiesMigration(comps, unpackedData)
+            else:
+                # iterating of np is not fast...
+                for c, val, linkedDim in itertools.zip_longest(comps, unpackedData, linkedDims, fillvalue=""):
+                    try:
+                        if linkedDim != "":
+                            c.p[paramName] = linkedDim
+                        else:
+                            c.p[paramName] = val
+                    except AssertionError as ae:
+                        # happens when a param was deprecated but being loaded from old DB
+                        runLog.warning(
+                            f"{str(ae)}\nSkipping load of invalid param `{paramName}` (possibly loading from old DB)\n"
+                        )
 
     def getHistoryByLocation(
         self,
@@ -1443,6 +1447,22 @@ class Database:
 
         return resolved
 
+    @staticmethod
+    def _applyComponentNumberDensitiesMigration(comps, unpackedData):
+        """
+        Special migration from <= v0.5.1 component numberDensities parameter data type.
+
+        old format: dict[str: float]
+        new format: two numpy arrays
+        - nuclides = np.array(dtype="S6")
+        - numberDensities = np.array(dtype=np.float64)
+        """
+        for c, ndensDict in zip(comps, unpackedData):
+            nuclides = np.array(list(ndensDict.keys()), dtype="S6")
+            numberDensities = np.array(list(ndensDict.values()), dtype=np.float64)
+            c.p["nuclides"] = nuclides
+            c.p["numberDensities"] = numberDensities
+
 
 def packSpecialData(
     arrayData: [np.ndarray, JaggedArray], paramName: str
diff --git a/armi/nucDirectory/nuclideBases.py b/armi/nucDirectory/nuclideBases.py
index 59686b9f..79918d66 100644
--- a/armi/nucDirectory/nuclideBases.py
+++ b/armi/nucDirectory/nuclideBases.py
@@ -19,29 +19,25 @@ framework and applications.
     :id: I_ARMI_ND_ISOTOPES0
     :implements: R_ARMI_ND_ISOTOPES
 
-    The :py:mod:`nuclideBases <armi.nucDirectory.nuclideBases>` module defines
-    the :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>`
-    class which is used to organize and store metadata about each nuclide. The
-    metadata is read from ``nuclides.dat`` file in the ARMI resources folder,
-    which contains metadata for 4,614 isotopes. The module also contains classes
-    for special types of nuclides, including :py:class:`DummyNuclideBase
+    The :py:mod:`nuclideBases <armi.nucDirectory.nuclideBases>` module defines the
+    :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>` class which is used to
+    organize and store metadata about each nuclide. The metadata is read from ``nuclides.dat`` file
+    in the ARMI resources folder, which contains metadata for 4,614 isotopes. The module also
+    contains classes for special types of nuclides, including :py:class:`DummyNuclideBase
     <armi.nucDirectory.nuclideBases.DummyNuclideBase>` for dummy nuclides,
-    :py:class:`LumpNuclideBase
-    <armi.nucDirectory.nuclideBases.LumpNuclideBase>`, for lumped fission
+    :py:class:`LumpNuclideBase <armi.nucDirectory.nuclideBases.LumpNuclideBase>`, for lumped fission
     product nuclides, and :py:class:`NaturalNuclideBase
-    <armi.nucDirectory.nuclideBases.NaturalNuclideBase>` for when data is given
-    collectively for an element at natural abundance rather than for individual
-    isotopes.
+    <armi.nucDirectory.nuclideBases.NaturalNuclideBase>` for when data is given collectively for an
+    element at natural abundance rather than for individual isotopes.
 
-    The :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>`
-    provides a data structure for information about a single nuclide, including
-    the atom number, atomic weight, element, isomeric state, half-life, and
-    name.
+    The :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>` provides a data
+    structure for information about a single nuclide, including the atom number, atomic weight,
+    element, isomeric state, half-life, and name.
 
-    The :py:mod:`nuclideBases <armi.nucDirectory.nuclideBases>` module provides
-    a factory and associated functions for instantiating the
-    :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>` objects
-    and building the global nuclide dictionaries, including:
+    The :py:mod:`nuclideBases <armi.nucDirectory.nuclideBases>` module provides a factory and
+    associated functions for instantiating the
+    :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>` objects and building the
+    global nuclide dictionaries, including:
 
     * ``instances`` (list of nuclides)
     * ``byName`` (keyed by name, e.g., ``U235``)
@@ -54,8 +50,7 @@ framework and applications.
     * ``byMcnpId`` (keyed by MCNP ID, e.g., ``92235``)
     * ``byAAAZZZSId`` (keyed by AAAZZZS, e.g., ``2350920``)
 
-The nuclide class structure is outlined :ref:`here
-<nuclide-bases-class-diagram>`.
+The nuclide class structure is outlined :ref:`here <nuclide-bases-class-diagram>`.
 
 .. _nuclide-bases-class-diagram:
 
@@ -104,15 +99,13 @@ from armi import context, runLog
 from armi.nucDirectory import transmutations
 from armi.utils.units import HEAVY_METAL_CUTOFF_Z
 
-# used to prevent multiple applications of burn chains, which would snowball
-# unphysically. This is a bit of a crutch for the global state that is the nuclide
-# directory.
+# Used to prevent multiple applications of burn chains, which would snowball unphysically. This is a
+# bit of a crutch for the global state that is the nuclide directory.
 burnChainImposed = False
 
 instances = []
-# The elements must be imported after the instances list is established
-# to allow for simultaneous initialization of the nuclides and elements
-# together to maintain self-consistency.
+# The elements must be imported after the instances list is established to allow for simultaneous
+# initialization of the nuclides and elements together to maintain self-consistency.
 from armi.nucDirectory import elements  # noqa: E402
 
 # Dictionary of INuclides by the INuclide.name for fast indexing
@@ -293,25 +286,21 @@ class INuclide(NuclideInterface):
         Indicates excitement, 1 is more excited than 0.
 
     abundance : float
-        Isotopic fraction of a naturally occurring nuclide. The sum of all nuclide
-        abundances for a naturally occurring element should be 1.0. This is atom
-        fraction, not mass fraction.
+        Isotopic fraction of a naturally occurring nuclide. The sum of all nuclide abundances for a
+        naturally occurring element should be 1.0. This is atom fraction, not mass fraction.
 
     name : str
         ARMI's unique name for the given nuclide.
 
     label : str
-        ARMI's unique 4 character label for the nuclide.
-        These are not human readable, but do not lose any information.
-        The label is effectively the
-        :attr:`Element.symbol `armi.nucDirectory.elements.Element.symbol`
-        padded to two characters, plus the mass number (A) in base-26 (0-9, A-Z).
-        Additional support for meta-states is provided by adding 100 * the state
-        to the mass number (A).
+        ARMI's unique 4 character label for the nuclide. These are not human readable, but do not
+        lose any information. The label is effectively the
+        :attr:`Element.symbol `armi.nucDirectory.elements.Element.symbol` padded to two characters,
+        plus the mass number (A) in base-26 (0-9, A-Z). Additional support for meta-states is
+        provided by adding 100 * the state to the mass number (A).
 
     nuSF : float
-        Neutrons released per spontaneous fission.
-        This should probably be moved at some point.
+        Neutrons released per spontaneous fission. This should probably be moved at some point.
     """
 
     fissile = ["U235", "PU239", "PU241", "AM242M", "CM244", "U233"]
@@ -436,8 +425,8 @@ class INuclide(NuclideInterface):
     def getDecay(self, decayType):
         """Get a :py:class:`~armi.nucDirectory.transmutations.DecayMode`.
 
-        Retrieve the first :py:class:`~armi.nucDirectory.transmutations.DecayMode`
-        matching the specified decType.
+        Retrieve the first :py:class:`~armi.nucDirectory.transmutations.DecayMode` matching the
+        specified decType.
 
         Parameters
         ----------
@@ -510,15 +499,12 @@ class NuclideBase(INuclide, IMcnpNuclide):
         :id: I_ARMI_ND_ISOTOPES1
         :implements: R_ARMI_ND_ISOTOPES
 
-        The :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>`
-        class provides a data structure for information about a single nuclide,
-        including the atom number, atomic weight, element, isomeric state,
-        half-life, and name. The class contains static methods for creating an
-        internal ARMI name or label for a nuclide. There are instance methods
-        for generating the nuclide ID for external codes, e.g. MCNP or Serpent,
-        and retrieving the nuclide ID for MC\ :sup:`2`-2 or MC\ :sup:`2`-3.
-        There are also instance methods for generating an AAAZZZS ID and an ENDF
-        MAT number.
+        Instances of this class provide a data structure for information about a single nuclide,
+        including the atom number, atomic weight, element, isomeric state, half-life, and name. The
+        class contains static methods for creating an internal ARMI name or label for a nuclide.
+        There are instance methods for generating the nuclide ID for external codes, e.g. MCNP or
+        Serpent, and retrieving the nuclide ID for MC\ :sup:`2`-2 or MC\ :sup:`2`-3. There are also
+        instance methods for generating an AAAZZZS ID and an ENDF MAT number.
     """
 
     def __init__(self, element, a, weight, abundance, state, halflife):
@@ -554,10 +540,10 @@ class NuclideBase(INuclide, IMcnpNuclide):
         """
         Make label for nuclide base.
 
-        The logic causes labels for things with A<10 to be zero padded like H03 or tritium
-        instead of H3. This avoids the metastable tritium collision which would look
-        like elemental HE. It also allows things like MO100 to be held within 4 characters,
-        which is a constraint of the ISOTXS format if we append 2 characters for XS type.
+        The logic causes labels for things with A<10 to be zero padded like H03 or tritium instead
+        of H3. This avoids the metastable tritium collision which would look like elemental HE. It
+        also allows things like MO100 to be held within 4 characters, which is a constraint of the
+        ISOTXS format if we append 2 characters for XS type.
         """
         # len(e.symbol) is 1 or 2 => a % (either 1000 or 100)
         #                         => gives exact a, or last two digits.
@@ -591,9 +577,9 @@ class NuclideBase(INuclide, IMcnpNuclide):
             :implements: R_ARMI_ND_ISOTOPES
 
             This method returns the ``mcc2id`` attribute of a
-            :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>`
-            instance.  This attribute is initially populated by reading from the
-            mcc-nuclides.yaml file in the ARMI resources folder.
+            :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>` instance. This
+            attribute is initially populated by reading from the mcc-nuclides.yaml file in the ARMI
+            resources folder.
         """
         return self.mcc2id
 
@@ -609,9 +595,9 @@ class NuclideBase(INuclide, IMcnpNuclide):
             :implements: R_ARMI_ND_ISOTOPES
 
             This method returns the ``mcc3idEndfbVII0`` attribute of a
-            :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>`
-            instance.  This attribute is initially populated by reading from the
-            mcc-nuclides.yaml file in the ARMI resources folder.
+            :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>` instance. This
+            attribute is initially populated by reading from the mcc-nuclides.yaml file in the ARMI
+            resources folder.
         """
         return self.mcc3idEndfbVII0
 
@@ -623,9 +609,9 @@ class NuclideBase(INuclide, IMcnpNuclide):
             :implements: R_ARMI_ND_ISOTOPES
 
             This method returns the ``mcc3idEndfbVII1`` attribute of a
-            :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>`
-            instance.  This attribute is initially populated by reading from the
-            mcc-nuclides.yaml file in the ARMI resources folder.
+            :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>` instance. This
+            attribute is initially populated by reading from the mcc-nuclides.yaml file in the ARMI
+            resources folder.
         """
         return self.mcc3idEndfbVII1
 
@@ -637,10 +623,9 @@ class NuclideBase(INuclide, IMcnpNuclide):
             :id: I_ARMI_ND_ISOTOPES4
             :implements: R_ARMI_ND_ISOTOPES
 
-            This method generates the MCNP ID for an isotope using the standard
-            MCNP format based on the atomic number A, number of protons Z, and
-            excited state. The implementation includes the special rule for
-            Am-242m, which is 95242. 95642 is used for the less common ground
+            This method generates the MCNP ID for an isotope using the standard MCNP format based on
+            the atomic number A, number of protons Z, and excited state. The implementation includes
+            the special rule for Am-242m, which is 95242. 95642 is used for the less common ground
             state Am-242.
 
         Returns
@@ -672,10 +657,9 @@ class NuclideBase(INuclide, IMcnpNuclide):
             :id: I_ARMI_ND_ISOTOPES5
             :implements: R_ARMI_ND_ISOTOPES
 
-            This method generates the AAAZZZS format ID for an isotope. Where
-            AAA is the mass number, ZZZ is the atomic number, and S is the
-            isomeric state. This is a general format independent of any code that
-            precisely defines an isotope or isomer.
+            This method generates the AAAZZZS format ID for an isotope. Where AAA is the mass
+            number, ZZZ is the atomic number, and S is the isomeric state. This is a general format
+            independent of any code that precisely defines an isotope or isomer.
 
         Notes
         -----
@@ -690,8 +674,8 @@ class NuclideBase(INuclide, IMcnpNuclide):
         Returns
         -------
         id: str
-            The ID of this nuclide based on it's elemental name, weight,
-            and state, eg ``U-235``, ``Te-129m``,
+            The ID of this nuclide based on it's elemental name, weight, and state, eg ``U-235``,
+            ``Te-129m``.
         """
         symbol = self.element.symbol.capitalize()
         return "{}-{}{}".format(symbol, self.a, "m" if self.state else "")
@@ -700,12 +684,11 @@ class NuclideBase(INuclide, IMcnpNuclide):
         """
         Gets the ENDF MAT number.
 
-        MAT numbers are defined as described in section 0.4.1 of the NJOY manual.
-        Basically, it's Z * 100 + I where I is an isotope number. I=25 is defined
-        as the lightest known stable isotope of element Z, so for Uranium,
-        Z=92 and I=25 refers to U234. The values of I go up by 3 for each
-        mass number, so U235 is 9228. This leaves room for three isomeric
-        states of each nuclide.
+        MAT numbers are defined as described in section 0.4.1 of the NJOY manual. Basically, it's
+        Z * 100 + I where I is an isotope number. I=25 is defined as the lightest known stable
+        isotope of element Z, so for Uranium, Z=92 and I=25 refers to U234. The values of I go up by
+        3 for each mass number, so U235 is 9228. This leaves room for three isomeric states of each
+        nuclide.
 
         Returns
         -------
@@ -734,9 +717,9 @@ class NaturalNuclideBase(INuclide, IMcnpNuclide):
 
     Notes
     -----
-    This is meant to represent the combination of all naturally occurring nuclides
-    within an element. The abundance is forced to zero here so that it does not
-    have any interactions with the NuclideBase objects.
+    This is meant to represent the combination of all naturally occurring nuclides within an
+    element. The abundance is forced to zero here so that it does not have any interactions with the
+    NuclideBase objects.
     """
 
     def __init__(self, name, element):
@@ -824,8 +807,8 @@ class DummyNuclideBase(INuclide):
 
     Notes
     -----
-    This may be used to store mass from a depletion calculation, specifically
-    in the instances where the burn chain is truncated.
+    This may be used to store mass from a depletion calculation, specifically in the instances where
+    the burn chain is truncated.
     """
 
     def __init__(self, name, weight):
@@ -966,7 +949,7 @@ class LumpNuclideBase(INuclide):
         return self.mcc3idEndfbVII1
 
 
-def initReachableActiveNuclidesThroughBurnChain(numberDensityDict, activeNuclides):
+def initReachableActiveNuclidesThroughBurnChain(nuclides, numberDensities, activeNuclides):
     """
     March through the depletion chain and find all nuclides that can be reached by depleting nuclides passed in.
 
@@ -974,19 +957,22 @@ def initReachableActiveNuclidesThroughBurnChain(numberDensityDict, activeNuclide
 
     Parameters
     ----------
-    numberDensityDict : dict
-        Starting number densities.
-
+    nuclides : np.array, dtype="S6"
+        Starting array of nuclide names
+    numberDensities : np.array, dtype=np.float64
+        Starting array of number densities
     activeNuclides : OrderedSet
         Active nuclides defined on the reactor blueprints object. See: armi.reactor.blueprints.py
     """
     if not burnChainImposed:
-        return
+        return nuclides, numberDensities
 
     missingActiveNuclides = set()
     memo = set()
-    difference = set(numberDensityDict).difference(memo)
+    nucNames = [nucName.decode() for nucName in nuclides]
+    difference = set(nucNames).difference(memo)
     while any(difference):
+        newNucs = set()
         nuclide = difference.pop()
         memo.add(nuclide)
         # Skip the nuclide if it is not `active` in the burn-chain
@@ -997,20 +983,28 @@ def initReachableActiveNuclidesThroughBurnChain(numberDensityDict, activeNuclide
 
         for interaction in nuclideObj.trans + nuclideObj.decays:
             try:
-                # Interaction nuclides can only be added to the number density
-                # dictionary if they are a part of the user-defined active nuclides
+                # Interaction nuclides can only be added to the number density arrays if they
+                # are a part of the user-defined active nuclides
                 productNuclide = interaction.getPreferredProduct(activeNuclides)
-                if productNuclide not in numberDensityDict:
-                    numberDensityDict[productNuclide] = 0.0
+                if productNuclide not in nucNames:
+                    newNucs.add(productNuclide.encode())
             except KeyError:
                 # Keep track of the first production nuclide
                 missingActiveNuclides.add(interaction.productNuclides)
 
-        difference = set(numberDensityDict).difference(memo)
+        # add the new nuclides to the number density arrays
+        newNDens = np.zeros(len(newNucs), dtype=np.float64)
+        nuclides = np.append(nuclides, list(newNucs))
+        numberDensities = np.append(numberDensities, newNDens)
+
+        nucNames = [nucName.decode() for nucName in nuclides]
+        difference = set(nucNames).difference(memo)
 
     if burnChainImposed and missingActiveNuclides:
         _failOnMissingActiveNuclides(missingActiveNuclides)
 
+    return nuclides, numberDensities
+
 
 def _failOnMissingActiveNuclides(missingActiveNuclides):
     """Raise ValueError with notification of which nuclides to include in the burn-chain."""
@@ -1089,9 +1083,9 @@ def single(predicate):
     """
     Return a single :py:class:`INuclide` object meeting the specified condition.
 
-    Similar to :py:func:`where`, this function uses a lambda input to filter
-    the :py:attr:`INuclide instances <instances>`. If there is not 1 and only
-    1 match for the specified condition, an exception is raised.
+    Similar to :py:func:`where`, this function uses a lambda input to filter the
+    :py:attr:`INuclide instances <instances>`. If there is not 1 and only 1 match for the specified
+    condition, an exception is raised.
 
     Examples
     --------
@@ -1118,8 +1112,8 @@ def changeLabel(nuclideBase, newLabel):
 
     Notes
     -----
-    Since nuclide objects are defined and stored globally, any change to the
-    attributes will be maintained.
+    Since nuclide objects are defined and stored globally, any change to the attributes will be
+    maintained.
     """
     nuclideBase.label = newLabel
     byLabel[newLabel] = nuclideBase
@@ -1137,16 +1131,16 @@ def imposeBurnChain(burnChainStream):
     Notes
     -----
     You cannot impose a burn chain twice. Doing so would require that you clean out the
-    transmutations and decays from all the module-level nuclide bases, which generally
-    requires that you rebuild them. But rebuilding those is not an option because some
-    of them get set as class-level attributes and would be orphaned. If a need to change
-    burn chains mid-run re-arises, then a better nuclideBase-level burnchain cleanup
-    should be implemented so the objects don't have to change identity.
+    transmutations and decays from all the module-level nuclide bases, which generally requires that
+    you rebuild them. But rebuilding those is not an option because some of them get set as class-
+    level attributes and would be orphaned. If a need to change burn chains mid-run re-arises, then
+    a better nuclideBase-level burnchain cleanup should be implemented so the objects don't have to
+    change identity.
 
     Notes
     -----
-    We believe the transmutation information would probably be better stored on a
-    less fundamental place (e.g. not on the NuclideBase).
+    We believe the transmutation information would probably be better stored on a less fundamental
+    place (e.g. not on the NuclideBase).
 
     See Also
     --------
@@ -1174,19 +1168,17 @@ def factory():
     Reads data files to instantiate the :py:class:`INuclides <INuclide>`.
 
     Reads NIST, MC**2 and burn chain data files to instantiate the :py:class:`INuclides <INuclide>`.
-    Also clears and fills in the
-    :py:data:`~armi.nucDirectory.nuclideBases.instances`,
+    Also clears and fills in the :py:data:`~armi.nucDirectory.nuclideBases.instances`,
     :py:data:`byName`, :py:attr:`byLabel`, :py:data:`byMcc3IdEndfbVII0`, and
-    :py:data:`byMcc3IdEndfbVII1` module attributes. This method is automatically run upon
-    loading the module, hence it is not usually necessary to re-run it unless there is a
-    change to the data files, which should not happen during run time, or a *bad*
-    :py:class`INuclide` is created.
+    :py:data:`byMcc3IdEndfbVII1` module attributes. This method is automatically run upon loading
+    the module, hence it is not usually necessary to re-run it unless there is a change to the data
+    files, which should not happen during run time, or a *bad* :py:class`INuclide` is created.
 
     Notes
     -----
-    This cannot be run more than once. NuclideBase instances are used throughout the ARMI
-    ecosystem and are even class attributes in some cases. Re-instantiating them would orphan
-    any existing ones and break everything.
+    This cannot be run more than once. NuclideBase instances are used throughout the ARMI ecosystem
+    and are even class attributes in some cases. Re-instantiating them would orphan any existing
+    ones and break everything.
     """
     if len(instances) != 0:
         raise RuntimeError(
@@ -1218,13 +1210,11 @@ def addNuclideBases():
         :id: I_ARMI_ND_DATA0
         :implements: R_ARMI_ND_DATA
 
-        This function reads the ``nuclides.dat`` file from the ARMI resources
-        folder. This file contains metadata for 4,614 nuclides, including
-        number of protons, number of neutrons, atomic number, excited
-        state, element symbol, atomic mass, natural abundance, half-life,
-        and spontaneous fission yield. The data in ``nuclides.dat`` have been
-        collected from multiple different sources; the references are given
-        in comments at the top of that file.
+        This function reads the ``nuclides.dat`` file from the ARMI resources folder. This file
+        contains metadata for 4,614 nuclides, including number of protons, number of neutrons,
+        atomic number, excited state, element symbol, atomic mass, natural abundance, half-life, and
+        spontaneous fission yield. The data in ``nuclides.dat`` have been collected from multiple
+        different sources; the references are given in comments at the top of that file.
     """
     with open(os.path.join(context.RES, "nuclides.dat")) as f:
         for line in f:
@@ -1260,13 +1250,7 @@ def __addNaturalNuclideBases():
 
 
 def __addDummyNuclideBases():
-    """
-    Generates a set of dummy nuclides.
-
-    Notes
-    -----
-    These nuclides can be used to truncate a depletion / burn-up chain within the
-    """
+    """Generates a set of dummy nuclides."""
     DummyNuclideBase(name="DUMP1", weight=10.0)
     DummyNuclideBase(name="DUMP2", weight=240.0)
 
@@ -1287,14 +1271,13 @@ def readMCCNuclideData():
         :id: I_ARMI_ND_DATA1
         :implements: R_ARMI_ND_DATA
 
-        This function reads the mcc-nuclides.yaml file from the ARMI resources
-        folder. This file contains the MC\ :sup:`2`-2 ID (from ENDF/B-V.2) and MC\ :sup:`2`-3 ID
-        (from ENDF/B-VII.0) for all nuclides in MC\ :sup:`2`. The ``mcc2id``,
-        ``mcc3idEndfVII0``, and  ``mcc3idEndfVII1`` attributes of each :py:class:`NuclideBase
-        <armi.nucDirectory.nuclideBases.NuclideBase>` instance are updated as
-        the data is read, and the global dictionaries ``byMcc2Id``
-        ``byMcc3IdEndfVII0`` and ``byMcc3IdEndfVII1`` are populated with the nuclide bases
-        keyed by their corresponding ID for each code.
+        This function reads the mcc-nuclides.yaml file from the ARMI resources folder. This file
+        contains the MC\ :sup:`2`-2 ID (from ENDF/B-V.2) and MC\ :sup:`2`-3 ID (from ENDF/B-VII.0)
+        for all nuclides in MC\ :sup:`2`. The ``mcc2id``, ``mcc3idEndfVII0``, and ``mcc3idEndfVII1``
+        attributes of each :py:class:`NuclideBase <armi.nucDirectory.nuclideBases.NuclideBase>`
+        instance are updated as the data is read, and the global dictionaries ``byMcc2Id``
+        ``byMcc3IdEndfVII0`` and ``byMcc3IdEndfVII1`` are populated with the nuclide bases keyed by
+        their corresponding ID for each code.
     """
     global byMcc2Id
     global byMcc3Id
@@ -1333,18 +1316,16 @@ def updateNuclideBasesForSpecialCases():
         :implements: R_ARMI_ND_ISOTOPES
 
         This function updates the keys for the :py:class:`NuclideBase
-        <armi.nucDirectory.nuclideBases.NuclideBase>` instances for Am-242m and
-        Am-242 in the ``byName`` and ``byDBName`` global dictionaries.  This
-        function associates the more common isomer Am-242m with the name
-        "AM242", and uses "AM242G" to denote the ground state.
+        <armi.nucDirectory.nuclideBases.NuclideBase>` instances for Am-242m and Am-242 in the
+        ``byName`` and ``byDBName`` global dictionaries. This function associates the more common
+        isomer Am-242m with the name "AM242", and uses "AM242G" to denote the ground state.
 
     Notes
     -----
-    This function is specifically added to change the definition of
-    `AM242` to refer to its metastable isomer, `AM242M` by default. `AM242M`
-    is most common isomer of `AM242` and is typically the desired isomer
-    when being requested rather than than the ground state (i.e., S=0) of
-    `AM242`.
+    This function is specifically added to change the definition of `AM242` to refer to its
+    metastable isomer, `AM242M` by default. `AM242M` is most common isomer of `AM242` and is
+    typically the desired isomer when being requested rather than than the ground state (i.e., S=0)
+    of `AM242`.
     """
     # Change the name of `AM242` to specific represent its ground state.
     am242g = byName["AM242"]
diff --git a/armi/physics/neutronics/crossSectionGroupManager.py b/armi/physics/neutronics/crossSectionGroupManager.py
index f1e3d93b..05dbd592 100644
--- a/armi/physics/neutronics/crossSectionGroupManager.py
+++ b/armi/physics/neutronics/crossSectionGroupManager.py
@@ -482,12 +482,21 @@ def getBlockNuclideTemperatureAvgTerms(block, allNucNames):
 
     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
-        ]
+        if component.p.nuclides is None:
+            return [0.0 for _nuc in allNucNames]
+
+        allByteNucs = [nucName.encode() for nucName in allNucNames]
+        ndens = []
+        nucCopy = np.array(component.p.nuclides)
+        nDensCopy = np.array(component.p.numberDensities)
+        reverseIndex = {nuc: i for i, nuc in enumerate(nucCopy)}
+        for nuc in allByteNucs:
+            i = reverseIndex.get(nuc, -1)
+            if i >= 0:
+                ndens.append(max(nDensCopy[i], TRACE_NUMBER_DENSITY))
+            else:
+                ndens.append(0.0)
+        return ndens
 
     vol = block.getVolume()
     components, volFracs = zip(*block.getVolumeFractions())
diff --git a/armi/reactor/blueprints/componentBlueprint.py b/armi/reactor/blueprints/componentBlueprint.py
index 027da8ca..ba1ad394 100644
--- a/armi/reactor/blueprints/componentBlueprint.py
+++ b/armi/reactor/blueprints/componentBlueprint.py
@@ -415,7 +415,7 @@ def insertDepletableNuclideKeys(c, blueprint):
         :py:func:`~armi.nucDirectory.nuclideBases.initReachableActiveNuclidesThroughBurnChain` to
         apply the user-specifications in the "nuclide flags" section of the blueprints to the
         Component such that all active isotopes and derivatives of those isotopes in the burn chain
-        are initialized to have an entry in the component's ``numberDensities`` dictionary.
+        are initialized to have an entry in the component's ``nuclides`` array.
 
         Note that certain case settings, including ``fpModel`` and ``fpModelLibrary``, may trigger
         modifications to the active nuclides specified by the user in the "nuclide flags" section of
@@ -433,7 +433,11 @@ def insertDepletableNuclideKeys(c, blueprint):
     """
     if c.hasFlags(Flags.DEPLETABLE):
         # depletable components, whether auto-derived or explicitly flagged need expanded nucs
-        nuclideBases.initReachableActiveNuclidesThroughBurnChain(
+        (
+            c.p.nuclides,
+            c.p.numberDensities,
+        ) = nuclideBases.initReachableActiveNuclidesThroughBurnChain(
+            c.p.nuclides,
             c.p.numberDensities,
             blueprint.activeNuclides,
         )
diff --git a/armi/reactor/components/component.py b/armi/reactor/components/component.py
index f2b61cc8..ae79691b 100644
--- a/armi/reactor/components/component.py
+++ b/armi/reactor/components/component.py
@@ -13,7 +13,8 @@
 # limitations under the License.
 
 """
-Components represent geometric objects within an assembly such as fuel, bond, coolant, ducts, wires, etc.
+Components represent geometric objects within an assembly such as fuel, bond, coolant, ducts, wires,
+etc.
 
 This module contains the abstract definition of a Component.
 """
@@ -82,9 +83,8 @@ class _DimensionLink(tuple):
     def __str__(self):
         """Return a string representation of a dimension link.
 
-        These look like ``otherComponentName.otherDimensionName``.
-        For example, if a link were to a ``fuel`` component's
-        ``od`` param, the link would render as ``fuel.od``.
+        These look like ``otherComponentName.otherDimensionName``. For example, if a link were to a
+        ``fuel`` component's ``od`` param, the link would render as ``fuel.od``.
         """
         return f"{self[0].name}.{self[1]}"
 
@@ -93,10 +93,10 @@ class ComponentType(composites.CompositeModelType):
     """
     ComponetType is a metaclass for storing and initializing Component subclass types.
 
-    The construction of Component subclasses is being done through factories for ease of
-    user input.  As a consequence, the ``__init__`` methods' arguments need to be known
-    in order to conform them to the correct format. Additionally, the constructors
-    arguments can be used to determine the Component subclasses dimensions.
+    The construction of Component subclasses is being done through factories for ease of user input.
+    As a consequence, the ``__init__`` methods' arguments need to be known in order to conform them
+    to the correct format. Additionally, the constructors arguments can be used to determine the
+    Component subclasses dimensions.
 
     Warning
     -------
@@ -217,7 +217,7 @@ class Component(composites.Composite, metaclass=ComponentType):
         self.temperatureInC = Thot
         self.material = None
         self.setProperties(material)
-        self.applyMaterialMassFracsToNumberDensities()  # not necessary when duplicating...
+        self.applyMaterialMassFracsToNumberDensities()  # not necessary when duplicating
         self.setType(name)
         self.p.mergeWith = mergeWith
         self.p.customIsotopicsName = isotopics
@@ -239,10 +239,11 @@ class Component(composites.Composite, metaclass=ComponentType):
 
     def __lt__(self, other):
         """
-        True if a circle encompassing this object has a smaller diameter than one encompassing another component.
+        True if a circle encompassing this object has a smaller diameter than one encompassing
+        another component.
 
-        If the bounding circles for both components have identical size, then revert to checking the inner
-        diameter of each component for sorting.
+        If the bounding circles for both components have identical size, then revert to checking the
+        inner diameter of each component for sorting.
 
         This allows sorting because the Python sort functions only use this method.
         """
@@ -286,9 +287,7 @@ class Component(composites.Composite, metaclass=ComponentType):
             Some Components are fluids and are thus defined by the shapes surrounding
             them. This method cycles through each dimension defining the border of this
             Component and converts the name of that Component to a link to the object
-            itself. This series of links is then used downstream to resolve
-            dimensional information.
-
+            itself. This series of links is then used downstream to resolve dimensional information.
         """
         for dimName in self.DIMENSION_NAMES:
             value = self.p[dimName]
@@ -346,7 +345,7 @@ class Component(composites.Composite, metaclass=ComponentType):
         # setting number densities can sometimes change thermal expansion behavior. Call again so
         # the material has access to its parent's comp when providing the reference initial density.
         densityBasedOnParentComposition = self.material.getProperty("pseudoDensity", Tc=self.temperatureInC)
-        self.p.numberDensities = densityTools.getNDensFromMasses(
+        self.p.nuclides, self.p.numberDensities = densityTools.getNDensFromMasses(
             densityBasedOnParentComposition, self.material.massFrac
         )
 
@@ -364,12 +363,10 @@ class Component(composites.Composite, metaclass=ComponentType):
 
         Notes
         -----
-        Call before setTemperature since we need old hot temp.
-        This works well if there is only 1 solid component.
-        If there are multiple components expanding at different rates during thermal
-        expansion this becomes more complicated and, and axial expansion should be used.
+        Call before setTemperature since we need old hot temp. This works well if there is only 1
+        solid component. If there are multiple components expanding at different rates during
+        thermal expansion this becomes more complicated and, and axial expansion should be used.
         Multiple expansion rates cannot trivially be accommodated.
-        See AxialExpansionChanger.
         """
         self.changeNDensByFactor(1.0 / self.getHeightFactor(newHot))
 
@@ -472,11 +469,9 @@ class Component(composites.Composite, metaclass=ComponentType):
 
         Notes
         -----
-        ``self.p.volume`` is not set until this method is called,
-        so under most circumstances it is probably not safe to
-        access ``self.p.volume`` directly. This is because not
-        all components (e.g., ``DerivedShape``) can compute
-        their volume during initialization.
+        ``self.p.volume`` is not set until this method is called, so under most circumstances it is
+        probably not safe to access ``self.p.volume`` directly. This is because not all components
+        (e.g., ``DerivedShape``) can compute their volume during initialization.
         """
         if self.p.volume is None:
             self._updateVolume()
@@ -518,9 +513,9 @@ class Component(composites.Composite, metaclass=ComponentType):
         """
         Check for negative area and warn/error when appropriate.
 
-        Negative component area is allowed for Void materials (such as gaps)
-        which may be placed between components that will overlap during thermal expansion
-        (such as liners and cladding and annular fuel).
+        Negative component area is allowed for Void materials (such as gaps) which may be placed
+        between components that will overlap during thermal expansion (such as liners and cladding
+        and annular fuel).
 
         Overlapping is allowed to maintain conservation of atoms while sticking close to the
         as-built geometry. Modules that need true geometries will have to handle this themselves.
@@ -623,7 +618,9 @@ class Component(composites.Composite, metaclass=ComponentType):
 
         This includes anything that has been specified in here, including trace nuclides.
         """
-        return list(self.p.numberDensities.keys())
+        if self.p.nuclides is None:
+            return []
+        return [nucName.decode() for nucName in self.p.nuclides]
 
     def getNumberDensity(self, nucName):
         """
@@ -639,14 +636,59 @@ class Component(composites.Composite, metaclass=ComponentType):
         number density : float
             number density in atoms/bn-cm.
         """
-        return self.p.numberDensities.get(nucName, 0.0)
+        i = np.where(self.p.nuclides == nucName.encode())[0]
+        if i.size > 0:
+            return self.p.numberDensities[i[0]]
+        else:
+            return 0.0
 
     def getNuclideNumberDensities(self, nucNames):
         """Return a list of number densities for the nuc names requested."""
-        return [self.p.numberDensities.get(nucName, 0.0) for nucName in nucNames]
+        if isinstance(nucNames, (list, tuple, np.ndarray)):
+            byteNucs = np.asanyarray(nucNames, dtype="S6")
+        else:
+            byteNucs = [nucName.encode() for nucName in nucNames]
+
+        if self.p.numberDensities is None:
+            return np.zeros(len(byteNucs), dtype=np.float64)
+
+        # trivial case where nucNames is the full set of nuclides in the same order
+        if np.array_equal(byteNucs, self.p.nuclides):
+            return np.array(self.p.numberDensities)
+
+        if len(byteNucs) < len(self.p.nuclides) / 10:
+            return self._getNumberDensitiesArray(byteNucs)
+
+        nDensDict = dict(zip(self.p.nuclides, self.p.numberDensities))
+        return [nDensDict.get(nuc, 0.0) for nuc in byteNucs]
+
+    def _getNumberDensitiesArray(self, byteNucs):
+        """
+        Get number densities using direct array lookup.
+
+        When only a small subset of nuclide number densities are requested, it is
+        likely faster to lookup the index for each nuclide than to recreate the
+        entire dictionary for a lookup.
+
+        Parameters
+        ----------
+        byteNucs : np.ndarray, dtype="S6"
+            List of nuclides for which to retrieve number densities, as encoded byte strings
+        """
+        ndens = np.zeros(len(byteNucs), dtype=np.float64)
+        nuclides = self.p.nuclides
+        numberDensities = self.p.numberDensities
+
+        # if it's just a small subset of nuclides, use np.where for direct index lookup
+        for i, nuc in enumerate(byteNucs):
+            j = np.where(nuclides == nuc)[0]
+            if j.size > 0:
+                ndens[i] = numberDensities[j[0]]
+        return ndens
 
     def _getNdensHelper(self):
-        return dict(self.p.numberDensities)
+        nucs = self.getNuclides()
+        return dict(zip(nucs, self.p.numberDensities)) if len(nucs) > 0 else {}
 
     def setName(self, name):
         """Components use name for type and name."""
@@ -661,9 +703,9 @@ class Component(composites.Composite, metaclass=ComponentType):
             :id: I_ARMI_COMP_NUCLIDE_FRACS0
             :implements: R_ARMI_COMP_NUCLIDE_FRACS
 
-            The method allows a user or plugin to set the number density of a Component.
-            It also indicates to other processes that may depend on a Component's
-            status about this change via the ``assigned`` attribute.
+            The method allows a user or plugin to set the number density of a Component. It also
+            indicates to other processes that may depend on a Component's status about this change
+            via the ``assigned`` attribute.
 
         Parameters
         ----------
@@ -706,8 +748,8 @@ class Component(composites.Composite, metaclass=ComponentType):
         numberDensities : dict
             nucName: ndens pairs.
         wipe : bool, optional
-            Controls whether the old number densities are wiped. Any nuclide densities not
-            provided in numberDensities will be effectively set to 0.0.
+            Controls whether the old number densities are wiped. Any nuclide densities not provided
+            in numberDensities will be effectively set to 0.0.
 
         Notes
         -----
@@ -721,11 +763,12 @@ class Component(composites.Composite, metaclass=ComponentType):
         the pre-perturbed component's volume/dimensions.
 
         This has no effect if the material thermal expansion has no dependence on component
-        composition. If this is not desired, `self.p.numberDensities` can be set directly.
+        composition. If this is not desired, `self.p.numberDensities` and `self.p.nuclides` can be
+        set directly.
         """
-        # prepare to change the densities with knowledge that dims could change due to
-        # material thermal expansion dependence on composition
-        if len(self.p.numberDensities) > 0:
+        # prepare to change the densities with knowledge that dims could change due to material
+        # thermal expansion dependence on composition
+        if self.p.numberDensities.size > 0:
             dLLprev = self.material.linearExpansionPercent(Tc=self.temperatureInC) / 100.0
             materialExpansion = True
         else:
@@ -735,15 +778,31 @@ class Component(composites.Composite, metaclass=ComponentType):
         try:
             vol = self.getVolume()
         except (AttributeError, TypeError):
-            # either no parent to get height or parent's height is None
-            # which would be AttributeError and TypeError respectively, but other errors could be possible
+            # Either no parent to get height or parent's height is None. Which would be
+            # AttributeError and TypeError respectively, but other errors could be possible.
             vol = None
             area = self.getArea()
 
         # change the densities
         if wipe:
-            self.p.numberDensities = {}  # clear things not passed
-        self.p.numberDensities.update(numberDensities)
+            self.p.numberDensities = None  # clear things not passed
+            self.p.nuclides = None  # clear things not passed
+            self.p.nuclides = np.asanyarray(list(numberDensities.keys()), dtype="S6")
+            self.p.numberDensities = np.array(list(numberDensities.values()))
+        else:
+            newNucs = []
+            newNumDens = []
+            nucs = self.p.nuclides
+            ndens = self.p.numberDensities
+            for nucName, dens in numberDensities.items():
+                i = np.where(nucs == nucName.encode())[0]
+                if i.size > 0:
+                    ndens[i[0]] = dens
+                else:
+                    newNucs.append(nucName.encode())
+                    newNumDens.append(dens)
+            self.p.nuclides = np.append(nucs, newNucs)
+            self.p.numberDensities = np.append(ndens, newNumDens)
 
         # check if thermal expansion changed
         dLLnew = self.material.linearExpansionPercent(Tc=self.temperatureInC) / 100.0
@@ -759,15 +818,15 @@ class Component(composites.Composite, metaclass=ComponentType):
                 factor = area / self.getArea()
             self.changeNDensByFactor(factor)
 
-        # since we're updating the object the param points to but not the param itself, we have to inform
-        # the param system to flag it as modified so it properly syncs during ``syncMpiState``.
+        # since we are updating the object the param points to but not the param itself, we have to
+        # inform the param system to flag it as modified so it syncs during ``syncMpiState``.
         self.p.assigned = parameters.SINCE_ANYTHING
         self.p.paramDefs["numberDensities"].assigned = parameters.SINCE_ANYTHING
 
     def changeNDensByFactor(self, factor):
         """Change the number density of all nuclides within the object by a multiplicative factor."""
-        newDensities = {nuc: dens * factor for nuc, dens in self.p.numberDensities.items()}
-        self.p.numberDensities = newDensities
+        if self.p.numberDensities is not None:
+            self.p.numberDensities *= factor
         self._changeOtherDensParamsByFactor(factor)
 
     def _changeOtherDensParamsByFactor(self, factor):
@@ -789,8 +848,7 @@ class Component(composites.Composite, metaclass=ComponentType):
         Notes
         -----
         Getting mass enrichment on any level higher than this is ambiguous because you may have
-        enriched boron in one pin and enriched uranium in another and blending those doesn't make
-        sense.
+        enriched boron in one pin and uranium in another and blending those doesn't make sense.
         """
         if self.material.enrichedNuclide is None:
             raise ValueError(f"Cannot get enrichment of {self.material} because `enrichedNuclide` is not defined.")
@@ -834,10 +892,16 @@ class Component(composites.Composite, metaclass=ComponentType):
             The mass in grams.
         """
         volume = self.getVolume() / (self.parent.getSymmetryFactor() if self.parent else 1.0)
-        nuclideNames = self._getNuclidesFromSpecifier(nuclideNames)
-        # densities comes from self.p.numberDensities
-        densities = self.getNuclideNumberDensities(nuclideNames)
-        nDens = {nuc: dens for nuc, dens in zip(nuclideNames, densities)}
+        if nuclideNames is None:
+            nDens = self._getNdensHelper()
+        else:
+            nuclideNames = self._getNuclidesFromSpecifier(nuclideNames)
+            # densities comes from self.p.numberDensities
+            if len(nuclideNames) > 0:
+                densities = self.getNuclideNumberDensities(nuclideNames)
+                nDens = dict(zip(nuclideNames, densities))
+            else:
+                nDens = {}
         return densityTools.calculateMassDensity(nDens) * volume
 
     def setDimension(self, key, val, retainLink=False, cold=True):
@@ -853,11 +917,11 @@ class Component(composites.Composite, metaclass=ComponentType):
             dimension is for a cold configuration or not. If it is not for a cold
             configuration, the thermal expansion factor is considered when setting the dimension.
 
-            If the ``retainLink`` argument is ``True``, any Components linked to this
-            one will also have its dimensions changed consistently. After a dimension
-            is updated, the ``clearLinkedCache`` method is called which sets the
-            volume of this Component to ``None``. This ensures that when the volume is
-            next accessed it is recomputed using the updated dimensions.
+            If the ``retainLink`` argument is ``True``, any Components linked to this one will also
+            have its dimensions changed consistently. After a dimension is updated, the
+            ``clearLinkedCache`` method is called which sets the volume of this Component to
+            ``None``. This ensures that when the volume is next accessed it is recomputed using the
+            updated dimensions.
 
         Parameters
         ----------
@@ -866,8 +930,8 @@ class Component(composites.Composite, metaclass=ComponentType):
         val : float
             The value to set on the dimension
         retainLink : bool, optional
-            If True, the val will be applied to the dimension of linked
-            component which indirectly changes this component's dimensions.
+            If True, the val will be applied to the dimension of linked component which indirectly
+            changes this component's dimensions.
         cold : bool, optional
             If True sets the component cold dimension to the specified value.
         """
@@ -892,12 +956,11 @@ class Component(composites.Composite, metaclass=ComponentType):
             :id: I_ARMI_COMP_DIMS
             :implements: R_ARMI_COMP_DIMS
 
-            Due to thermal expansion, Component dimensions depend on their temperature.
-            This method retrieves a dimension from the Component at a particular
-            temperature, if provided. If the Component is a LinkedComponent then the
-            dimensions are resolved to ensure that any thermal expansion that has
-            occurred to the Components that the LinkedComponent depends on is reflected
-            in the returned dimension.
+            Due to thermal expansion, Component dimensions depend on their temperature. This method
+            retrieves a dimension from the Component at a particular temperature, if provided. If
+            the Component is a LinkedComponent then the dimensions are resolved to ensure that any
+            thermal expansion that has occurred to the Components that the LinkedComponent depends
+            on is reflected in the returned dimension.
 
         Parameters
         ----------
@@ -927,10 +990,9 @@ class Component(composites.Composite, metaclass=ComponentType):
 
         Notes
         -----
-        The inner circle is meaningful for annular shapes, i.e., circle with non-zero ID,
-        hexagon with non-zero IP, etc. For shapes with corners (e.g., hexagon, rectangle, etc)
-        the inner circle intersects the corners of the inner bound, opposed to intersecting
-        the "flats".
+        The inner circle is meaningful for annular shapes, i.e., circle with non-zero ID, hexagon
+        with non-zero IP, etc. For shapes with corners (e.g., hexagon, rectangle, etc) the inner
+        circle intersects the corners of the inner bound, opposed to intersecting the "flats".
         """
         raise NotImplementedError
 
@@ -975,13 +1037,11 @@ class Component(composites.Composite, metaclass=ComponentType):
             :id: I_ARMI_COMP_EXPANSION0
             :implements: R_ARMI_COMP_EXPANSION
 
-            This method enables the calculation of the thermal expansion factor
-            for a given material. If the material is solid, the difference
-            between ``T0`` and ``Tc`` is used to calculate the thermal expansion
-            factor. If a solid material does not have a linear expansion factor
-            defined and the temperature difference is greater than
-            a predetermined tolerance, an
-            error is raised. Thermal expansion of fluids or custom materials is
+            This method enables the calculation of the thermal expansion factor for a given
+            material. If the material is solid, the difference between ``T0`` and ``Tc`` is used to
+            calculate the thermal expansion factor. If a solid material does not have a linear
+            expansion factor defined and the temperature difference is greater than a predetermined
+            tolerance, an error is raised. Thermal expansion of fluids or custom materials is
             neglected, currently.
 
         Parameters
@@ -1094,8 +1154,7 @@ class Component(composites.Composite, metaclass=ComponentType):
         aMerge = compToMergeWith.getArea()
         meNDens = {nucName: aMe / aMerge * self.getNumberDensity(nucName) for nucName in self.getNuclides()}
         mergeNDens = {nucName: compToMergeWith.getNumberDensity(nucName) for nucName in compToMergeWith.getNuclides()}
-        # set the new homogenized number densities from both. Allow
-        # overlapping nuclides.
+        # set the new homogenized number densities from both. Allow overlapping nuclides.
         for nucName in set(meNDens) | set(mergeNDens):
             compToMergeWith.setNumberDensity(nucName, (meNDens.get(nucName, 0.0) + mergeNDens.get(nucName, 0.0)))
 
@@ -1107,8 +1166,8 @@ class Component(composites.Composite, metaclass=ComponentType):
         """
         Create and store a backup of the state.
 
-        This needed to be overridden due to linked components which actually have a parameter value of another
-        ARMI component.
+        This needed to be overridden due to linked components which actually have a parameter value
+        of another ARMI component.
         """
         linkedDims = self._getLinkedDimsAndValues()
         composites.Composite.backUp(self)
@@ -1193,20 +1252,19 @@ class Component(composites.Composite, metaclass=ComponentType):
         """
         Return the multigroup neutron flux in [n/cm^2/s].
 
-        The first entry is the first energy group (fastest neutrons). Each additional
-        group is the next energy group, as set in the ISOTXS library.
+        The first entry is the first energy group (fastest neutrons). Each additional group is the
+        next energy group, as set in the ISOTXS library.
 
         Parameters
         ----------
         adjoint : bool, optional
             Return adjoint flux instead of real
         average : bool, optional
-            If True, will return average flux between latest and previous. Doesn't work
-            for pin detailed.
+            If True, will return average flux between latest and previous. Does not work for pin
+            detailed.
         volume: float, optional
-            The volume-integrated flux is divided by volume before
-            being returned. The user may specify a volume here, or the function
-            will obtain the block volume directly.
+            The volume-integrated flux is divided by volume before being returned. The user may
+            specify a volume here, or the function will obtain the block volume directly.
         gamma : bool, optional
             Whether to return the neutron flux or the gamma flux.
 
@@ -1274,16 +1332,15 @@ class Component(composites.Composite, metaclass=ComponentType):
         Returns
         -------
         np.ndarray
-            A ``(N, nGroup)`` array of pin multigroup fluxes, where ``N`` is the
-            equivalent to the multiplicity of the component (``self.p.mult``)
-            and ``nGroup`` is the number of energy groups of the flux.
+            A ``(N, nGroup)`` array of pin multigroup fluxes, where ``N`` is the equivalent to the
+            multiplicity of the component (``self.p.mult``) and ``nGroup`` is the number of energy
+            groups of the flux.
 
         Raises
         ------
         ValueError
-            If the location(s) of the component are not aligned with pin indices
-            from the block. This would happen if this component is not actually
-            a pin.
+            If the location(s) of the component are not aligned with pin indices from the block.
+            This would happen if this component is not actually a pin.
         """
         # Get the (i, j, k) location of all pins from the parent block
         indicesAll = {(loc.i, loc.j): i for i, loc in enumerate(self.parent.getPinLocations())}
@@ -1345,7 +1402,7 @@ class Component(composites.Composite, metaclass=ComponentType):
 
         See Also
         --------
-        armi.physics.neutronics.fissionProductModel.lumpedFissionProduct.LumpedFissionProduct : LFP object
+        armi.physics.neutronics.fissionProductModel.lumpedFissionProduct.LumpedFissionProduct
         """
         if self.parent:
             return self.parent.getLumpedFissionProductCollection()
@@ -1376,15 +1433,14 @@ class Component(composites.Composite, metaclass=ComponentType):
     def finalizeLoadingFromDB(self):
         """Apply any final actions after creating the component from database.
 
-        This should **only** be called internally by the database loader. Otherwise
-        some properties could be doubly applied.
+        This should **only** be called internally by the database loader. Otherwise some properties
+        could be doubly applied.
 
-        This exists because the theoretical density is initially defined as a material
-        modification, and then stored as a Material attribute. When reading from blueprints,
-        the blueprint loader sets the theoretical density parameter from the Material
-        attribute. Component parameters are also set when reading from the database.
-        But, we need to set the Material attribute so routines that fetch a material's
-        density property account for the theoretical density.
+        This exists because the theoretical density is initially defined as a material modification,
+        and then stored as a Material attribute. When reading from blueprints, the blueprint loader
+        sets the theoretical density parameter from the Material attribute. Component parameters are
+        also set when reading from the database. But, we need to set the Material attribute so
+        routines that fetch a material's density property account for the theoretical density.
         """
         self.material.adjustTD(self.p.theoreticalDensityFrac)
 
diff --git a/armi/reactor/components/componentParameters.py b/armi/reactor/components/componentParameters.py
index 945cfb89..14b2c0f2 100644
--- a/armi/reactor/components/componentParameters.py
+++ b/armi/reactor/components/componentParameters.py
@@ -59,10 +59,18 @@ def getComponentParameterDefinitions():
 
         pb.defParam(
             "numberDensities",
+            setter=isNumpyArray("numberDensities"),
             units=f"#/(bn*{units.CM})",
             description="Number densities of each nuclide.",
         )
 
+        pb.defParam(
+            "nuclides",
+            setter=isNumpyArray("nuclides"),
+            units=units.UNITLESS,
+            description="Nuclide names corresponding to numberDensities array.",
+        )
+
         pb.defParam(
             "detailedNDens",
             setter=isNumpyArray("detailedNDens"),
diff --git a/armi/reactor/composites.py b/armi/reactor/composites.py
index eccbea84..1357e149 100644
--- a/armi/reactor/composites.py
+++ b/armi/reactor/composites.py
@@ -430,11 +430,11 @@ class ArmiObject(metaclass=CompositeModelType):
         """
         Make a clean copy of this object.
 
-        .. warning:: Be careful with inter-object dependencies. If one object contains a
-            reference to another object which contains links to the entire hierarchical
-            tree, memory can fill up rather rapidly. Weak references are designed to help
-            with this problem.
-
+        Warning
+        -------
+        Be careful with inter-object dependencies. If one object contains a reference to another
+        object which contains links to the entire hierarchical tree, memory can fill up rather
+        rapidly. Weak references are designed to help with this problem.
         """
         raise NotImplementedError
 
@@ -448,11 +448,9 @@ class ArmiObject(metaclass=CompositeModelType):
         """
         Obtain a value from the cache.
 
-        Cached values can be used to temporarily store frequently read but
-        long-to-compute values.  The practice is generally discouraged because it's
-        challenging to make sure to properly invalidate the cache when the state
-        changes.
-
+        Cached values can be used to temporarily store frequently read but long-to-compute values.
+        The practice is generally discouraged because it's challenging to make sure to properly
+        invalidate the cache when the state changes.
         """
         return self.cached.get(name, None)
 
@@ -935,9 +933,7 @@ class ArmiObject(metaclass=CompositeModelType):
         elif isinstance(nucSpec, (str)):
             nuclideNames = [nucSpec]
         elif isinstance(nucSpec, list):
-            nuclideNames = []
-            for ns in nucSpec:
-                nuclideNames.extend(self._getNuclidesFromSpecifier(ns))
+            nuclideNames = nucSpec
         else:
             raise TypeError(f"nucSpec={nucSpec} is an invalid specifier. It is a {type(nucSpec)}")
 
@@ -1240,12 +1236,7 @@ class ArmiObject(metaclass=CompositeModelType):
             # there are no children so no volume or number density
             return [0.0] * len(nucNames)
 
-        densListForEachComp = []
-        for c in self:
-            numberDensityDict = c.getNumberDensities()
-            densListForEachComp.append([numberDensityDict.get(nuc, 0.0) for nuc in nucNames])
-        nucDensForEachComp = np.array(densListForEachComp)  # c x n
-
+        nucDensForEachComp = np.array([c.getNuclideNumberDensities(nucNames) for c in self])  # c x n
         return volumes.dot(nucDensForEachComp) / totalVol
 
     def _getNdensHelper(self):
@@ -1880,7 +1871,7 @@ class ArmiObject(metaclass=CompositeModelType):
             The total heavy metal number (atom) density in atoms/bn-cm.
         """
         hmNuclides = [nuclide for nuclide in self.getNuclides() if nucDir.isHeavyMetal(nuclide)]
-        hmDens = sum(self.getNuclideNumberDensities(hmNuclides))
+        hmDens = np.sum(self.getNuclideNumberDensities(hmNuclides))
         return hmDens
 
     def getFPMass(self):
@@ -2273,16 +2264,20 @@ class ArmiObject(metaclass=CompositeModelType):
 
         Parameters
         ----------
-        elementalNuclide : :class:`armi.nucDirectory.nuclideBases.NaturalNuclide`
-            natural nuclide to replace.
+        elementalNuclide : :class:`armi.nucDirectory.nuclideBases.NaturalNuclide` natural nuclide to
+            replace.
         """
         natName = elementalNuclide.name
         for component in self.iterComponents():
             elementalDensity = component.getNumberDensity(natName)
             if elementalDensity == 0.0:
                 continue
-            component.setNumberDensity(natName, 0.0)  # clear the elemental
-            del component.p.numberDensities[natName]
+
+            keepIndex = np.where(component.p.nuclides != natName.encode())[0]
+            newNuclides = [nuc.decode() for nuc in component.p.nuclides[keepIndex]]
+            newNDens = component.p.numberDensities[keepIndex]
+            component.updateNumberDensities(dict(zip(newNuclides, newNDens)), wipe=True)
+
             # add in isotopics
             for natNuc in elementalNuclide.getNaturalIsotopics():
                 component.setNumberDensity(natNuc.name, elementalDensity * natNuc.abundance)
@@ -2343,8 +2338,7 @@ class Composite(ArmiObject):
     ``spatialLocator`` and ``spatialGrid`` parameters. The ``spatialLocator`` is a numpy
     triple representing either:
 
-    1. Indices in the parent's ``spatialGrid`` (for lattices, etc.), used when the dtype
-    is int.
+    1. Indices in the parent's ``spatialGrid`` (for lattices, etc.), used when the dtype is int.
 
     2. Coordinates in the parent's universe in cm, used when the dtype is float.
 
@@ -2581,17 +2575,16 @@ class Composite(ArmiObject):
             :id: I_ARMI_CMP1
             :implements: R_ARMI_CMP
 
-            This method retrieves all children within a given Composite object. Children
-            of any generation can be retrieved. This is achieved by visiting all
-            children and calling this method recursively for each generation requested.
+            This method retrieves all children within a given Composite object. Children of any
+            generation can be retrieved. This is achieved by visiting all children and calling this
+            method recursively for each generation requested.
 
-            If the method is called with ``includeMaterials``, it will additionally
-            include information about the material for each child. If a function is
-            supplied as the ``predicate`` argument, then this method will be used
-            to evaluate all children as a filter to include or not. For example, if the
-            caller of this method only desires children with a certain flag, or children
-            which only contain a certain material, then the ``predicate`` function
-            can be used to perform this filtering.
+            If the method is called with ``includeMaterials``, it will additionally include
+            information about the material for each child. If a function is supplied as the
+            ``predicate`` argument, then this method will be used to evaluate all children as a
+            filter to include or not. For example, if the caller of this method only desires
+            children with a certain flag, or children which only contain a certain material, then
+            the ``predicate`` function can be used to perform this filtering.
 
         Parameters
         ----------
@@ -2599,21 +2592,20 @@ class Composite(ArmiObject):
             Return all children of all levels.
 
         generationNum : int, optional
-            Which generation to return. 1 means direct children, 2 means children of
-            children. Setting this parameter will only return children of this
-            generation, not their parents. Default: Just return direct children.
+            Which generation to return. 1 means direct children, 2 means children of children.
+            Setting this parameter will only return children of this generation, not their parents.
+            Default: Just return direct children.
 
         includeMaterials : bool, optional
             Include the material properties
 
         predicate : callable, optional
-            An optional unary predicate to use for filtering results. This can be used
-            to request children of specific types, or with desired attributes. Not all
-            ArmiObjects have the same methods and members, so care should be taken to
-            make sure that the predicate executes gracefully in all cases (e.g., use
-            ``getattr(obj, "attribute", None)`` to access instance attributes). Failure
-            to meet the predicate only affects the object in question; children will
-            still be considered.
+            An optional unary predicate to use for filtering results. This can be used to request
+            children of specific types, or with desired attributes. Not all ArmiObjects have the
+            same methods and members, so care should be taken to make sure that the predicate
+            executes gracefully in all cases (e.g., use ``getattr(obj, "attribute", None)`` to
+            access instance attributes). Failure to meet the predicate only affects the object in
+            question; children will still be considered.
 
         See Also
         --------
@@ -2667,22 +2659,21 @@ class Composite(ArmiObject):
 
     def syncMpiState(self):
         """
-        Synchronize all parameters of this object and all children to all worker nodes
-        over the network using MPI.
+        Synchronize all parameters of this object and all children to all worker nodes over the
+        network using MPI.
 
-        In parallelized runs, if each process has its own copy of the entire reactor
-        hierarchy, this method synchronizes the state of all parameters on all objects.
+        In parallelized runs, if each process has its own copy of the entire reactor hierarchy, this
+        method synchronizes the state of all parameters on all objects.
 
         .. impl:: Composites can be synchronized across MPI threads.
             :id: I_ARMI_CMP_MPI
             :implements: R_ARMI_CMP_MPI
 
-            Parameters need to be handled properly during parallel code execution.This
-            method synchronizes all parameters of the composite object across all
-            processes by cycling through all the children of the Composite and ensuring
-            that their parameters are properly synchronized. If it fails to synchronize,
-            an error message is displayed which alerts the user to which Composite has
-            inconsistent data across the processes.
+            Parameters need to be handled properly during parallel code execution. This method
+            synchronizes all parameters of the composite object across all processes by cycling
+            through all the children of the Composite and ensuring that their parameters are
+            properly synchronized. If it fails to synchronize, an error message is displayed which
+            alerts the user to which Composite has inconsistent data across the processes.
 
         Returns
         -------
@@ -2786,8 +2777,8 @@ class Composite(ArmiObject):
         """
         Mark the composite and child parameters as synchronized across MPI.
 
-        We clear SINCE_LAST_DISTRIBUTE_STATE so that anything after this point will set
-        the SINCE_LAST_DISTRIBUTE_STATE flag, indicating it has been modified
+        We clear SINCE_LAST_DISTRIBUTE_STATE so that anything after this point will set the
+        SINCE_LAST_DISTRIBUTE_STATE flag, indicating it has been modified
         SINCE_LAST_DISTRIBUTE_STATE.
         """
         paramDefs = set()
@@ -2812,8 +2803,8 @@ class Composite(ArmiObject):
         Parameters
         ----------
         paramsToApply : iterable
-            Parameters that should be applied to the state after existing the state
-            retainer. All others will be reverted to their values upon entering.
+            Parameters that should be applied to the state after existing the state retainer. All
+            others will be reverted to their values upon entering.
 
         Notes
         -----
@@ -2825,8 +2816,8 @@ class Composite(ArmiObject):
         """
         Create and store a backup of the state.
 
-        This needed to be overridden due to linked components which actually have a
-        parameter value of another ARMI component.
+        This needed to be overridden due to linked components which actually have a parameter value
+        of another ARMI component.
         """
         self._backupCache = (self.cached, self._backupCache)
         self.cached = {}  # don't .clear(), using reference above!
@@ -2872,7 +2863,7 @@ class Composite(ArmiObject):
 
         See Also
         --------
-        armi.physics.neutronics.fissionProductModel.lumpedFissionProduct.LumpedFissionProduct : LFP object
+        armi.physics.neutronics.fissionProductModel.lumpedFissionProduct.LumpedFissionProduct
         """
         lfps = ArmiObject.getLumpedFissionProductCollection(self)
         if lfps is None:
@@ -2899,8 +2890,8 @@ class Composite(ArmiObject):
         """
         Returns the multigroup neutron tracklength in [n-cm/s].
 
-        The first entry is the first energy group (fastest neutrons). Each additional
-        group is the next energy group, as set in the ISOTXS library.
+        The first entry is the first energy group (fastest neutrons). Each additional group is the
+        next energy group, as set in the ISOTXS library.
 
         Parameters
         ----------
@@ -3034,7 +3025,7 @@ class Composite(ArmiObject):
     def getPuMoles(self):
         """Returns total number of moles of Pu isotopes."""
         nucNames = [nuc.name for nuc in elements.byZ[94].nuclides]
-        puN = sum(self.getNuclideNumberDensities(nucNames))
+        puN = np.sum(self.getNuclideNumberDensities(nucNames))
 
         return puN / units.MOLES_PER_CC_TO_ATOMS_PER_BARN_CM * self.getVolume()
 
diff --git a/armi/reactor/cores.py b/armi/reactor/cores.py
index b2919e66..312698c9 100644
--- a/armi/reactor/cores.py
+++ b/armi/reactor/cores.py
@@ -1173,15 +1173,20 @@ class Core(composites.Composite):
             fuelNuclides = set()
             structureNuclides = set()
             for c in self.iterComponents():
+                compNuclides = []
                 # get only nuclides with non-zero number density
                 # nuclides could be present at 0.0 density just for XS generation
-                nuclides = [nuc for nuc, dens in c.getNumberDensities().items() if dens > 0.0]
+                if c.p.numberDensities is None:
+                    continue
+                for nuc, dens in zip(c.p.nuclides, c.p.numberDensities):
+                    if dens > 0.0:
+                        compNuclides.append(nuc.decode())
                 if c.getName() == "coolant":
-                    coolantNuclides.update(nuclides)
+                    coolantNuclides.update(compNuclides)
                 elif "fuel" in c.getName():
-                    fuelNuclides.update(nuclides)
+                    fuelNuclides.update(compNuclides)
                 else:
-                    structureNuclides.update(nuclides)
+                    structureNuclides.update(compNuclides)
             structureNuclides -= coolantNuclides
             structureNuclides -= fuelNuclides
             remainingNuclides = set(self.parent.blueprints.allNuclidesInProblem) - structureNuclides - coolantNuclides
diff --git a/armi/utils/densityTools.py b/armi/utils/densityTools.py
index 38dd5b1e..197aec13 100644
--- a/armi/utils/densityTools.py
+++ b/armi/utils/densityTools.py
@@ -16,6 +16,8 @@
 
 from typing import Dict, List, Tuple
 
+import numpy as np
+
 from armi import runLog
 from armi.nucDirectory import elements, nucDir, nuclideBases
 from armi.utils import units
@@ -43,18 +45,22 @@ def getNDensFromMasses(rho, massFracs, normalize=False):
 
     Returns
     -------
-    numberDensities : dict
-        vector of number densities (#/bn-cm) keyed by their nuclide name
+    nuclides : np.ndarray[np.bytes_]
+        vector of nuclide names as byte strings
+    numberDensities : np.ndarray[np.float64]
+        vector of number densities (#/bn-cm) for each nuclide in nuclides
     """
     if normalize:
         massFracs = normalizeNuclideList(massFracs, normalization=normalize)
 
-    numberDensities = {}
+    nuclides = []
+    numberDensities = []
     rho = rho * units.MOLES_PER_CC_TO_ATOMS_PER_BARN_CM
     for nucName, massFrac in massFracs.items():
         atomicWeight = nuclideBases.byName[nucName].weight
-        numberDensities[nucName] = massFrac * rho / atomicWeight
-    return numberDensities
+        nuclides.append(nucName.encode())
+        numberDensities.append(massFrac * rho / atomicWeight)
+    return np.array(nuclides), np.array(numberDensities)
 
 
 def getMassFractions(numberDensities):
diff --git a/armi/utils/flags.py b/armi/utils/flags.py
index a09af8aa..0acaa0f7 100644
--- a/armi/utils/flags.py
+++ b/armi/utils/flags.py
@@ -15,42 +15,39 @@
 """
 A Flag class, similar to ``enum.Flag``.
 
-This is an alternate implementation of the standard-library ``enum.Flag`` class. We use
-this to implement :py:class:`armi.reactor.flags.Flags`. We used to use the
-standard-library implementation, but that became limiting when we wanted to make it
-possible for plugins to define their own flags; the standard implementation does not
-support extension. We also considered the ``aenum`` package, which permits extension of
-``Enum`` classes, but unfortunately does not support extension of ``Flags``. So, we had to
-make our own. This is a much simplified version of what comes with ``aenum``, but still
-provides most of the safety and functionality.
+This is an alternate implementation of the standard-library ``enum.Flag`` class. We use this to implement
+:py:class:`armi.reactor.flags.Flags`. We used to use the standard-library implementation, but that became limiting when
+we wanted to make it possible for plugins to define their own flags; the standard implementation does not support
+extension. We also considered the ``aenum`` package, which permits extension of ``Enum`` classes, but unfortunately does
+not support extension of ``Flags``. So, we had to make our own. This is a much simplified version of what comes with
+``aenum``, but still provides most of the safety and functionality.
 """
 
 import math
 from typing import Dict, List, Sequence, Tuple, Union
 
+from armi import runLog
+
 
 class auto:  # noqa: N801
     """
     Empty class for requesting a lazily-evaluated automatic field value.
 
-    This can be used to automatically provision a value for a field, when the specific
-    value does not matter.
+    This can be used to automatically provision a value for a field, when the specific value does not matter.
 
-    In the future, it would be nice to support some arithmetic for these so that
-    automatically-derived combinations of other automatically defined fields can be
-    specified as well.
+    In the future, it would be nice to support some arithmetic for these so that automatically-derived combinations of
+    other automatically defined fields can be specified as well.
     """
 
     def __iter__(self):
         """
         Dummy __iter__ implementation.
 
-        This is only needed to make mypy happy when it type checks things that have
-        FlagTypes in them, since these can normally be iterated over, but mypy doesn't
-        know that the metaclass consumes the autos.
+        This is only needed to make mypy happy when it type checks things that have FlagTypes in them, since these can
+        normally be iterated over, but mypy doesn't know that the metaclass consumes the autos.
         """
         raise NotImplementedError(
-            "__iter__() is not actually implemented on {}; it is only defined to appease mypy.".format(type(self))
+            f"__iter__() is not actually implemented on {type(self)}; it is only defined to appease mypy."
         )
 
 
@@ -58,14 +55,12 @@ class _FlagMeta(type):
     """
     Metaclass for defining new Flag classes.
 
-    This attempts to do the minimum required to make the Flag class and its subclasses
-    function properly. It mostly digests the class attributes, resolves automatic values
-    and creates instances of the class as it's own class attributes for each field. The
-    rest of the functionality lives in the base ``Flag`` class as plain-old code.
+    This attempts to do the minimum required to make the Flag class and its subclasses function properly. It mostly
+    digests the class attributes, resolves automatic values and creates instances of the class as it's own class
+    attributes for each field. The rest of the functionality lives in the base ``Flag`` class as plain-old code.
 
-    .. tip:: Because individual flags are defined as *class* attributes (as opposed to
-        instance attributes), we have to customize the way a Flag subclass itself is built,
-        which requires a metaclass.
+    .. tip:: Because individual flags are defined as *class* attributes (as opposed to instance attributes), we have to
+        customize the way a Flag subclass itself is built, which requires a metaclass.
     """
 
     def __new__(cls, name, bases, attrs):
@@ -110,9 +105,8 @@ class _FlagMeta(type):
         """
         Implement indexing at the class level.
 
-        This has to be done at the metaclass level, since the python interpreter looks
-        to ``type(klass).__getitem__(klass, key)``, which for an implementation of Flag
-        is this metaclass.
+        This has to be done at the metaclass level, since the python interpreter looks to
+        ``type(klass).__getitem__(klass, key)``, which for an implementation of Flag is this metaclass.
         """
         return cls(cls._nameToValue[key])
 
@@ -121,31 +115,26 @@ class Flag(metaclass=_FlagMeta):
     """
     A collection of bitwise flags.
 
-    This is intended to emulate ``enum.Flag``, except with the possibility of extension after the
-    class has been defined. Most docs for ``enum.Flag`` should be relevant here, but there are sure
-    to be occasional differences.
+    This is intended to emulate ``enum.Flag``, except with the possibility of extension after the class has been
+    defined. Most docs for ``enum.Flag`` should be relevant here, but there are sure to be occasional differences.
 
     .. impl:: No two flags have equivalence.
         :id: I_ARMI_FLAG_DEFINE
         :implements: R_ARMI_FLAG_DEFINE
 
-        A bitwise flag class intended to emulate the standard library's ``enum.Flag``, with the
-        added functionality that it allows for extension after the class has been defined. Each Flag
-        is unique; no two Flags are equivalent.
+        A bitwise flag class intended to emulate the standard library's ``enum.Flag``, with the added functionality that
+        it allows for extension after the class has been defined. Each Flag is unique; no two Flags are equivalent.
 
-        Note that while Python allows for arbitrary-width integers, exceeding the system-native
-        integer size can lead to challenges in storing data, e.g. in an HDF5 file. In this case, the
-        ``from_bytes()`` and ``to_bytes()`` methods are provided to represent a Flag's values in
-        smaller chunks so that writeability can be maintained.
+        Note that while Python allows for arbitrary-width integers, exceeding the system-native integer size can lead to
+        challenges in storing data, e.g. in an HDF5 file. In this case, the ``from_bytes()`` and ``to_bytes()`` methods
+        are provided to represent a Flag's values in smaller chunks so that writeability can be maintained.
 
     .. warning::
-        Python features arbitrary-width integers, allowing one to represent an
-        practically unlimited number of fields. *However*, including more flags than can
-        be represented in the system-native integer types may lead to strange behavior
-        when interfacing with non-pure Python code. For instance, exceeding 64 fields
-        makes the underlying value not trivially-storable in an HDF5 file. In such
-        circumstances, the ``from_bytes()`` and ``to_bytes()`` methods are available to
-        represent a Flag's values in smaller chunks.
+        Python features arbitrary-width integers, allowing one to represent an practically unlimited number of fields.
+        *However*, including more flags than can be represented in the system-native integer types may lead to strange
+        behavior when interfacing with non-pure Python code. For instance, exceeding 64 fields makes the underlying
+        value not trivially-storable in an HDF5 file. In such circumstances, the ``from_bytes()`` and ``to_bytes()``
+        methods are available to represent a Flag's values in smaller chunks.
     """
 
     _autoAt = None
@@ -161,13 +150,14 @@ class Flag(metaclass=_FlagMeta):
         for k, v in self._nameToValue.items():
             if self._value & v:
                 flagsOn.add(k)
+
         return flagsOn
 
     def __repr__(self):
-        return "<{}.{}: {}>".format(type(self).__name__, "|".join(self._flagsOn()), self._value)
+        return f"<{type(self).__name__}.{'|'.join(self._flagsOn())}: {self._value}>"
 
     def __str__(self):
-        return "{}.{}".format(type(self).__name__, "|".join(self._flagsOn()))
+        return f"{type(self).__name__}.{'|'.join(self._flagsOn())}"
 
     def __getstate__(self):
         return self._value
@@ -180,10 +170,13 @@ class Flag(metaclass=_FlagMeta):
         """
         Plug a new field into the Flags.
 
-        This makes sure everything is consistent and does error/collision checks. Mostly
-        useful for extending an existing class with more fields.
+        This makes sure everything is consistent and does error/collision checks. Mostly useful for extending an
+        existing class with more fields.
         """
-        assert value not in cls._nameToValue
+        if name in cls._nameToValue:
+            runLog.debug(f"The flag {name} already exists and does not need to be recreated.")
+            return
+
         cls._valuesTaken.add(value)
         cls._nameToValue[name] = value
         cls._width = math.ceil(len(cls._nameToValue) / 8)
@@ -224,22 +217,21 @@ class Flag(metaclass=_FlagMeta):
         Extend the Flags object with new fields.
 
         .. warning::
-            This alters the class that it is called upon! Existing instances should see
-            the new data, since classes are mutable.
+            This alters the class that it is called upon! Existing instances should see the new data, since classes are
+            mutable.
 
         .. impl:: Set of flags are extensible without loss of uniqueness.
             :id: I_ARMI_FLAG_EXTEND0
             :implements: R_ARMI_FLAG_EXTEND
 
-            A class method to extend a ``Flag`` with a vector of provided additional ``fields``,
-            with field names as keys, without loss of uniqueness. Values for the additional
-            ``fields`` can be explicitly specified, or an instance of ``auto`` can be supplied.
+            A class method to extend a ``Flag`` with a vector of provided additional ``fields``, with field names as
+            keys, without loss of uniqueness. Values for the additional ``fields`` can be explicitly specified, or an
+            instance of ``auto`` can be supplied.
 
         Parameters
         ----------
         fields : dict
-            A dictionary containing field names as keys, and their desired values, or an instance of
-            ``auto`` as values.
+            A dictionary containing field names as keys, and their desired values, or an instance of ``auto`` as values.
 
         Example
         -------
@@ -256,7 +248,10 @@ class Flag(metaclass=_FlagMeta):
         # add explicit values first, so that autos know about them
         for field, value in ((f, v) for f, v in fields.items() if isinstance(v, int)):
             cls._registerField(field, value)
+
+        # find auto values (ignore if they already exist)
         toResolve = [field for field, val in fields.items() if isinstance(val, auto)]
+        toResolve = [field for field in toResolve if field not in cls._nameToValue]
         resolved = cls._resolveAutos(toResolve)
         for field, value in resolved:
             cls._registerField(field, value)
@@ -265,13 +260,11 @@ class Flag(metaclass=_FlagMeta):
         """
         Return a byte stream representing the flag.
 
-        This is useful when storing Flags in a data type of limited size. Python ints
-        can be of arbitrary size, while most other systems can only represent integers
-        of 32 or 64 bits. For compatibility, this function allows to convert the flags
-        to a sequence of single-byte elements.
+        This is useful when storing Flags in a data type of limited size. Python ints can be of arbitrary size, while
+        most other systems can only represent integers of 32 or 64 bits. For compatibility, this function allows to
+        convert the flags to a sequence of single-byte elements.
 
-        Note that this uses snake_case to mimic the method on the Python-native int
-        type.
+        Note that this uses snake_case to mimic the method on the Python-native int type.
         """
         return self._value.to_bytes(self.width(), byteorder=byteorder)
 
@@ -298,15 +291,12 @@ class Flag(metaclass=_FlagMeta):
 
         Note
         ----
-        This is avoiding just ~ on the ``_value`` because it might not be safe.  Using
-        the int directly is slightly dangerous in that python ints are not of fixed
-        width, so the result of inverting one Flag might not be as wide as the result of
-        inverting another Flag. Typically, one would want to invert a Flag to create a
-        mask for unsetting a bit on another Flag, like ``f1 &= ~f2``. If ``f2`` is
-        narrower than ``f1`` the field of ones that you need to keep ``f1`` bits on
-        might not cover the width of ``f1``, erroneously turning off its upper bits. Not
-        sure if this was an issue before or not. Once things are working, might make
-        sense to play with this more.
+        This is avoiding just ~ on the ``_value`` because it might not be safe. Using the int directly is slightly
+        dangerous in that python ints are not of fixed width, so the result of inverting one Flag might not be as wide
+        as the result of inverting another Flag. Typically, one would want to invert a Flag to create a mask for
+        unsetting a bit on another Flag, like ``f1 &= ~f2``. If ``f2`` is narrower than ``f1`` the field of ones that
+        you need to keep ``f1`` bits on might not cover the width of ``f1``, erroneously turning off its upper bits. Not
+        sure if this was an issue before or not. Once things are working, might makes sense to play with this more.
         """
         new = self._value
         for _, val in self._nameToValue.items():
@@ -334,7 +324,6 @@ class Flag(metaclass=_FlagMeta):
         return hash(self._value)
 
 
-# Type alias to reliably check for a proper Flag type. This cannot just be `Flag`, since
-# mypy gets confused by `auto` because it doesn't go to the trouble of resolving them in
-# the metaclass.
+# Type alias to reliably check for a proper Flag type. This cannot just be `Flag`, since mypy gets confused by `auto`
+# because it doesn't go to the trouble of resolving them in the metaclass.
 FlagType = Union[Flag, auto]


Test patch:
diff --git a/armi/bookkeeping/db/tests/test_databaseInterface.py b/armi/bookkeeping/db/tests/test_databaseInterface.py
index 0d320810..98d9ac9d 100644
--- a/armi/bookkeeping/db/tests/test_databaseInterface.py
+++ b/armi/bookkeeping/db/tests/test_databaseInterface.py
@@ -400,9 +400,8 @@ class TestDatabaseReading(unittest.TestCase):
         cls.td = directoryChangers.TemporaryDirectoryChanger()
         cls.td.__enter__()
 
-        # The database writes the settings object to the DB rather
-        # than the original input file. This allows settings to be
-        # changed in memory like this and survive for testing.
+        # The database writes the settings object to the DB rather than the original input file.
+        # This allows settings to be changed in memory like this and survive for testing.
         newSettings = {"verbosity": "extra"}
         cls.nCycles = 2
         newSettings["nCycles"] = cls.nCycles
@@ -565,18 +564,8 @@ class TestDatabaseReading(unittest.TestCase):
         c1 = b1.getComponent(Flags.FUEL)
         c2 = b2.getComponent(Flags.FUEL)
 
-        self.assertIsInstance(c1.p.numberDensities, dict)
-        self.assertIsInstance(c2.p.numberDensities, dict)
-        keys1 = set(c1.p.numberDensities.keys())
-        keys2 = set(c2.p.numberDensities.keys())
-        self.assertEqual(keys1, keys2)
-
-        numDensVec1, numDensVec2 = [], []
-        for k in keys1:
-            numDensVec1.append(c1.p.numberDensities[k])
-            numDensVec2.append(c2.p.numberDensities[k])
-
-        assert_allclose(numDensVec1, numDensVec2)
+        for i, v1 in enumerate(c1.p.numberDensities):
+            self.assertAlmostEqual(v1, c2.p.numberDensities[i])
 
     def test_timesteps(self):
         with Database(self.dbName, "r") as db:
diff --git a/armi/materials/tests/test_air.py b/armi/materials/tests/test_air.py
index 90e9b6f1..ec2970a2 100644
--- a/armi/materials/tests/test_air.py
+++ b/armi/materials/tests/test_air.py
@@ -21,8 +21,8 @@ from armi.materials.air import Air
 from armi.utils import densityTools
 
 """
-Reference thermal physical properties from Table A.4 in Incropera, Frank P.,
-et al. Fundamentals of heat and mass transfer. Vol. 5. New York: Wiley, 2002.
+Reference thermal physical properties from Table A.4 in Incropera, Frank P., et al. Fundamentals of
+heat and mass transfer. Vol. 5. New York: Wiley, 2002.
 """
 
 REFERENCE_Tk = [
@@ -232,16 +232,22 @@ class Test_Air(unittest.TestCase):
         refO = 0.210748
         refAR = 0.004671
 
-        nDens = densityTools.getNDensFromMasses(air.pseudoDensity(Tk=300), air.massFrac)
+        nuclides, nDens = densityTools.getNDensFromMasses(air.pseudoDensity(Tk=300), air.massFrac)
 
-        error = math.fabs(nDens["C"] / sum(nDens.values()) - refC)
-        self.assertLess(error, 1e-4)
-        error = math.fabs(nDens["N"] / sum(nDens.values()) - refN)
-        self.assertLess(error, 1e-4)
-        error = math.fabs(nDens["O"] / sum(nDens.values()) - refO)
-        self.assertLess(error, 1e-4)
-        error = math.fabs(nDens["AR"] / sum(nDens.values()) - refAR)
-        self.assertLess(error, 1e-4)
+        diff = 1e-4
+        error = abs(nDens[0] / sum(nDens) - refC)
+        self.assertLess(error, diff)
+        error = abs(nDens[1] / sum(nDens) - refN)
+        self.assertLess(error, diff)
+        error = abs(nDens[2] / sum(nDens) - refO)
+        self.assertLess(error, diff)
+        error = abs(nDens[3] / sum(nDens) - refAR)
+        self.assertLess(error, diff)
+
+        self.assertEqual(nuclides[0].decode(), "C")
+        self.assertEqual(nuclides[1].decode(), "N")
+        self.assertEqual(nuclides[2].decode(), "O")
+        self.assertEqual(nuclides[3].decode(), "AR")
 
     def test_validRanges(self):
         air = Air()
diff --git a/armi/reactor/blueprints/tests/test_customIsotopics.py b/armi/reactor/blueprints/tests/test_customIsotopics.py
index 23661818..ce951dfa 100644
--- a/armi/reactor/blueprints/tests/test_customIsotopics.py
+++ b/armi/reactor/blueprints/tests/test_customIsotopics.py
@@ -17,6 +17,7 @@
 import unittest
 from logging import DEBUG
 
+import numpy as np
 import yamlize
 
 from armi import runLog, settings
@@ -352,12 +353,8 @@ assemblies:
     def test_unmodified(self):
         """Ensure that unmodified components have the correct isotopics."""
         fuel = self.a[0].getComponent(Flags.FUEL)
-        self.assertEqual(
-            self.numUZrNuclides,
-            len(fuel.p.numberDensities),
-            msg=fuel.p.numberDensities.keys(),
-        )
-        # Note this density does not come from the material but is based on number densities
+        self.assertEqual(self.numUZrNuclides, len(fuel.p.numberDensities))
+        # NOTE: This density does not come from the material but is based on number densities.
         self.assertAlmostEqual(15.5, fuel.density(), 0)  # i.e. it is not 19.1
 
     def test_massFractionsAreApplied(self):
@@ -372,7 +369,10 @@ assemblies:
         self.assertEqual(self.numCustomNuclides, len(fuel1.p.numberDensities))
         self.assertAlmostEqual(19.1, fuel1.density())
 
-        self.assertEqual(set(fuel2.p.numberDensities.keys()), set(fuel1.p.numberDensities.keys()))  # keys are same
+        # keys are same
+        keys1 = set([i for i, v in enumerate(fuel1.p.numberDensities) if v == 0.0])
+        keys2 = set([i for i, v in enumerate(fuel2.p.numberDensities) if v == 0.0])
+        self.assertEqual(keys1, keys2)
 
     def test_densitiesAppliedToNonCustomMaterials(self):
         """Ensure that a density can be set in custom isotopics for components using library materials."""
@@ -489,13 +489,16 @@ assemblies:
             :id: T_ARMI_MAT_USER_INPUT4
             :tests: R_ARMI_MAT_USER_INPUT
         """
-        # fuel blocks 2 and 4 should be the same, one is defined as mass fractions, and the other as number fractions
+        # fuel blocks 2 and 4 should be the same, one is defined as mass fractions, and the other as
+        # number fractions
         fuel2 = self.a[1].getComponent(Flags.FUEL)
         fuel4 = self.a[3].getComponent(Flags.FUEL)
         self.assertAlmostEqual(fuel2.density(), fuel4.density())
 
-        for nuc in fuel2.p.numberDensities.keys():
-            self.assertAlmostEqual(fuel2.p.numberDensities[nuc], fuel4.p.numberDensities[nuc])
+        keys2 = set([i for i, v in enumerate(fuel2.p.numberDensities) if v == 0.0])
+        keys4 = set([i for i, v in enumerate(fuel4.p.numberDensities) if v == 0.0])
+        self.assertEqual(keys2, keys4)
+        np.testing.assert_almost_equal(fuel2.p.numberDensities, fuel4.p.numberDensities)
 
     def test_numberDensities(self):
         """Ensure that the custom isotopics can be specified via number densities.
@@ -504,21 +507,22 @@ assemblies:
             :id: T_ARMI_MAT_USER_INPUT5
             :tests: R_ARMI_MAT_USER_INPUT
         """
-        # fuel blocks 2 and 5 should be the same, one is defined as mass fractions, and the other as number densities
+        # fuel blocks 2 and 5 should be the same, one is defined as mass fractions, and the other as
+        # number densities
         fuel2 = self.a[1].getComponent(Flags.FUEL)
         fuel5 = self.a[4].getComponent(Flags.FUEL)
         self.assertAlmostEqual(fuel2.density(), fuel5.density())
 
-        for nuc in fuel2.p.numberDensities.keys():
-            self.assertAlmostEqual(fuel2.p.numberDensities[nuc], fuel5.p.numberDensities[nuc])
+        for i, nuc in enumerate(fuel2.p.nuclides):
+            self.assertIn(nuc, fuel5.p.nuclides)
+            j = np.where(fuel5.p.nuclides == nuc)[0][0]
+            self.assertAlmostEqual(fuel2.p.numberDensities[i], fuel5.p.numberDensities[j])
 
     def test_numberDensitiesAnchor(self):
         fuel4 = self.a[4].getComponent(Flags.FUEL)
         fuel5 = self.a[5].getComponent(Flags.FUEL)
         self.assertAlmostEqual(fuel4.density(), fuel5.density())
-
-        for nuc in fuel4.p.numberDensities.keys():
-            self.assertAlmostEqual(fuel4.p.numberDensities[nuc], fuel5.p.numberDensities[nuc])
+        np.testing.assert_almost_equal(fuel4.p.numberDensities, fuel5.p.numberDensities)
 
     def test_expandedNatural(self):
         cs = settings.Settings()
diff --git a/armi/reactor/tests/test_assemblies.py b/armi/reactor/tests/test_assemblies.py
index f3ab1893..c2db7f38 100644
--- a/armi/reactor/tests/test_assemblies.py
+++ b/armi/reactor/tests/test_assemblies.py
@@ -461,6 +461,9 @@ class Assembly_TestCase(unittest.TestCase):
         self.assertEqual(reactiveHeight, 20.0)
 
     def test_getFissileMass(self):
+        for b in self.assembly:
+            b.p.massHmBOL = b.getHMMass()
+            b.p.enrichmentBOL = b.getFissileMassEnrich()
         cur = self.assembly.getFissileMass()
         ref = sum(bi.getMass(["U235", "PU239"]) for bi in self.assembly)
         self.assertAlmostEqual(cur, ref)
diff --git a/armi/reactor/tests/test_components.py b/armi/reactor/tests/test_components.py
index 33de821b..d8d1eacb 100644
--- a/armi/reactor/tests/test_components.py
+++ b/armi/reactor/tests/test_components.py
@@ -150,7 +150,10 @@ class TestComponentFactory(unittest.TestCase):
             component = components.factory(name, [], thisAttrs)
             duped = copy.deepcopy(component)
             for key, val in component.p.items():
-                if key not in ["area", "volume", "serialNum"]:
+                if key in ["numberDensities", "nuclides"]:
+                    for i in range(len(val)):
+                        self.assertEqual(val[i], duped.p[key][i])
+                elif key not in ["area", "volume", "serialNum"]:
                     # they get recomputed
                     self.assertEqual(
                         val,
@@ -173,9 +176,9 @@ class TestGeneralComponents(unittest.TestCase):
 
     def setUp(self, component=None):
         """
-        Most of the time nothing will be passed as `component` and the result will
-        be stored in self, but you can also pass a component object as `component`,
-        in which case the object will be returned with the `parent` attribute assigned.
+        Most of the time nothing will be passed as `component` and the result will be stored in
+        self, but you can also pass a component object as `component`, in which case the object will
+        be returned with the `parent` attribute assigned.
         """
 
         class _Parent:
@@ -811,12 +814,14 @@ class TestCircle(TestShapedComponent):
 
     def test_getNumberDensities(self):
         """Test that demonstrates that number densities can be retrieved on from component."""
-        self.component.p.numberDensities = {"NA23": 1.0}
+        self.component.p.numberDensities = np.ones(1, dtype=np.float64)
+        self.component.p.nuclides = np.array(["NA23"], dtype="S6")
         self.assertEqual(self.component.getNumberDensity("NA23"), 1.0)
 
     def test_changeNumberDensities(self):
         """Test that demonstrates that the number densities on a component can be modified."""
-        self.component.p.numberDensities = {"NA23": 1.0}
+        self.component.p.numberDensities = np.ones(1, dtype=np.float64)
+        self.component.p.nuclides = np.array(["NA23"], dtype="S6")
         self.component.p.detailedNDens = [1.0]
         self.component.p.pinNDens = [1.0]
         self.assertEqual(self.component.getNumberDensity("NA23"), 1.0)
@@ -895,12 +900,12 @@ class TestComponentExpansion(unittest.TestCase):
         # mass density is proportional to Fe number density and derived from
         # all the number densities and atomic masses
         self.assertAlmostEqual(
-            circle1.p.numberDensities[isotope] / circle2.p.numberDensities[isotope],
+            circle1.getNumberDensity(isotope) / circle2.getNumberDensity(isotope),
             circle1.density() / circle2.density(),
         )
 
-        # the colder one has more because it is the same cold outer diameter
-        # but it would be taller at the same temperature
+        # the colder one has more because it is the same cold outer diameter but it would be taller
+        # at the same temperature
         mass1 = circle1.density() * circle1.getArea() * hotHeight
         mass2 = circle2.density() * circle2.getArea() * hotHeight
         self.assertGreater(mass1, mass2)
@@ -954,10 +959,7 @@ class TestComponentExpansion(unittest.TestCase):
         circle1.setTemperature(self.tHot)
 
         # now its density is same as hot component
-        self.assertAlmostEqual(
-            circle1.density(),
-            circle2.density(),
-        )
+        self.assertAlmostEqual(circle1.density(), circle2.density())
 
         # show that mass is conserved after expansion
         circle1NewHotHeight = hotHeight * heightFactor
@@ -1003,8 +1005,7 @@ class TestComponentExpansion(unittest.TestCase):
                 circle.density(),
                 circle.material.density(Tc=circle.temperatureInC),
             )
-            # total mass consistent between hot and cold
-            # Hot height will be taller
+            # total mass consistent between hot and cold. Hot height will be taller
             hotHeight = coldHeight * circle.getThermalExpansionFactor()
             self.assertAlmostEqual(
                 coldHeight * circle.getArea(cold=True) * circle.material.density(Tc=circle.inputTemperatureInC),
diff --git a/armi/reactor/tests/test_composites.py b/armi/reactor/tests/test_composites.py
index 76dbec6e..0206e995 100644
--- a/armi/reactor/tests/test_composites.py
+++ b/armi/reactor/tests/test_composites.py
@@ -244,7 +244,7 @@ class TestCompositePattern(unittest.TestCase):
         c1 = [c.name for c in self.container]
         self.assertNotEqual(c0, c1)
 
-        # verify repeated sortings behave
+        # verify repeated sorting behave
         for _ in range(3):
             self.container.sort()
             ci = [c.name for c in self.container]
diff --git a/armi/utils/tests/test_densityTools.py b/armi/utils/tests/test_densityTools.py
index e4ca4916..1a85d955 100644
--- a/armi/utils/tests/test_densityTools.py
+++ b/armi/utils/tests/test_densityTools.py
@@ -15,6 +15,8 @@
 
 import unittest
 
+import numpy as np
+
 from armi.materials.material import Material
 from armi.materials.uraniumOxide import UO2
 from armi.nucDirectory import elements, nuclideBases
@@ -123,10 +125,12 @@ class TestDensityTools(unittest.TestCase):
             :id: T_ARMI_UTIL_MASS2N_DENS
             :tests: R_ARMI_UTIL_MASS2N_DENS
         """
-        nDens = densityTools.getNDensFromMasses(1, {"O": 1, "H": 2})
+        nucs, nDens = densityTools.getNDensFromMasses(1, {"O": 1, "H": 2})
+        O = np.where(nucs == "O".encode())[0]
+        H = np.where(nucs == "H".encode())[0]
 
-        self.assertAlmostEqual(nDens["O"], 0.03764, 5)
-        self.assertAlmostEqual(nDens["H"], 1.19490, 5)
+        self.assertAlmostEqual(nDens[O][0], 0.03764, 5)
+        self.assertAlmostEqual(nDens[H][0], 1.19490, 5)
 
     def test_getMassFractions(self):
         """Number densities to mass fraction."""
diff --git a/armi/utils/tests/test_flags.py b/armi/utils/tests/test_flags.py
index 95a73174..abcc26e1 100644
--- a/armi/utils/tests/test_flags.py
+++ b/armi/utils/tests/test_flags.py
@@ -15,6 +15,7 @@
 
 import unittest
 
+from armi.reactor.composites import FlagSerializer
 from armi.utils.flags import Flag, auto
 
 
@@ -152,3 +153,47 @@ class TestFlag(unittest.TestCase):
 
     def test_getitem(self):
         self.assertEqual(ExampleFlag["FOO"], ExampleFlag.FOO)
+
+    def test_duplicateFlags(self):
+        """Show that duplicate flags can be added and silently ignored."""
+
+        class F(Flag):
+            @classmethod
+            def len(cls):
+                return len(cls._nameToValue)
+
+        F.extend({"FLAG0": auto()})
+        for i in range(1, 12):
+            F.extend({f"FLAG{i}": auto()})
+            num = F.len()
+            F.extend({f"FLAG{i - 1}": auto()})
+            self.assertEqual(F.len(), num)
+
+            # While the next two lines do not assert anything, these lines used to raise an error.
+            # So these lines remain as proof against that error in the future.
+            ff = getattr(F, f"FLAG{i}")
+            FlagSerializer._packImpl(
+                [
+                    ff,
+                ],
+                F,
+            )
+            self.assertEqual(F.len(), num)
+
+    def test_soManyFlags(self):
+        """Show that many flags can be added without issue."""
+
+        class F(Flag):
+            @classmethod
+            def len(cls):
+                return len(cls._nameToValue)
+
+        for i in range(1, 100):
+            num = F.len()
+            flagName = f"FLAG{i}"
+            F.extend({flagName: auto()})
+            self.assertEqual(F.len(), num + 1)
+
+            flag = getattr(F, flagName)
+            flag.to_bytes()
+            self.assertEqual(F.len(), num + 1)