udpate density to account for material change

Author: onufer

armi/reactor/blocks.py

@@ -2441,10 +2441,18 @@ def getWettedPerimeter(self):
         )
 
         # flags pertaining to circular pin components where the exterior of the circle is wetted
-        wettedPinComponentFlags = (Flags.CLAD, Flags.WIRE, Flags.CLAD | Flags.DEPLETABLE, Flags.WIRE | Flags.DEPLETABLE)
+        wettedPinComponentFlags = (
+            Flags.CLAD,
+            Flags.WIRE,
+            Flags.CLAD | Flags.DEPLETABLE,
+            Flags.WIRE | Flags.DEPLETABLE,
+        )
 
         # flags pertaining to components where both the interior and exterior are wetted
-        wettedHollowComponentFlags = (Flags.DUCT | Flags.INNER, Flags.DUCT | Flags.INNER | Flags.DEPLETABLE,)
+        wettedHollowComponentFlags = (
+            Flags.DUCT | Flags.INNER,
+            Flags.DUCT | Flags.INNER | Flags.DEPLETABLE,
+        )
 
         # obtain all wetted components based on type
         wettedHollowHexagonComponents = []

armi/reactor/components/component.py

@@ -375,11 +375,20 @@ def applyMaterialMassFracsToNumberDensities(self):
         # `density` is 3D density
         # call getProperty to cache and improve speed
         density = self.material.getProperty("pseudoDensity", Tc=self.temperatureInC)
-
         self.p.numberDensities = densityTools.getNDensFromMasses(
             density, self.material.massFrac
         )
 
+        # Sometimes the material thermal expansion depends on its parents composition (eg Pu frac) so
+        # setting number densities can sometimes change thermal expansion behavior.
+        # so call again so that the material has access to its parents comp when providing the reference initial density.
+        densityBasedOnParentComposition = self.material.getProperty(
+            "pseudoDensity", Tc=self.temperatureInC
+        )
+        self.p.numberDensities = densityTools.getNDensFromMasses(
+            densityBasedOnParentComposition, self.material.massFrac
+        )
+
         # material needs to be expanded from the material's cold temp to hot,
         # not components cold temp, so we don't use mat.linearExpansionFactor or
         # component.getThermalExpansionFactor.
@@ -721,12 +730,7 @@ def setNumberDensity(self, nucName, val):
         val : float
             Number density to set in atoms/bn-cm (heterogeneous)
         """
-        self.p.numberDensities[nucName] = val
-        self.p.assigned = parameters.SINCE_ANYTHING
-        # necessary for syncMpiState
-        parameters.ALL_DEFINITIONS[
-            "numberDensities"
-        ].assigned = parameters.SINCE_ANYTHING
+        self.updateNumberDensities({nucName: val})
 
     def setNumberDensities(self, numberDensities):
         """
@@ -747,12 +751,19 @@ def setNumberDensities(self, numberDensities):
 
         Notes
         -----
-        We don't just call setNumberDensity for each nuclide because we don't want to call ``getVolumeFractions``
-        for each nuclide (it's inefficient).
+        Note that sometimes volume/dimensions can change due to the number density change when the material thermal
+        expansion depends on the component's composition (eg its plutonium fraction). In this case, changing the
+        density will implicitly change the area/volume. Since it its very difficult to predict the new dims ahead of time,
+        and perturbation/depletion calculations are almost exclusively done assuming constant volume,
+        the densities sent are automatically perturbed to conserve mass with the original dimensions.
+        That is, the components densities are not exactly as passed, but whatever they would need to be to preserve volume
+        integrated number densities (moles) from the pre-perturbed components volume/dims.
+        Note this has no effect if the material thermal expansion has no dependence on component composition fracs.
+        If this is not desired, `self.p.numberDensities` can be set directly.
         """
-        self.p.numberDensities = numberDensities
+        self.updateNumberDensities(numberDensities, wipe=True)
 
-    def updateNumberDensities(self, numberDensities):
+    def updateNumberDensities(self, numberDensities, wipe=False):
         """
         Set one or more multiple number densities. Leaves unlisted number densities alone.
 
@@ -761,12 +772,61 @@ def updateNumberDensities(self, numberDensities):
         numberDensities : dict
             nucName: ndens pairs.
 
+        Notes
+        -----
+        Note that sometimes volume/dimensions can change due to the number density change when the material thermal
+        expansion depends on the component's composition (eg its plutonium fraction). In this case, changing the
+        density will implicitly change the area/volume. Since it its very difficult to predict the new dims ahead of time,
+        and perturbation/depletion calculations are almost exclusively done assuming constant volume,
+        the densities sent are automatically perturbed to conserve mass with the original dimensions.
+        That is, the components densities are not exactly as passed, but whatever they would need to be to preserve volume
+        integrated number densities (moles) from the pre-perturbed components volume/dims.
+        Note this has no effect if the material thermal expansion has no dependence on component composition fracs.
+        If this is not desired, `self.p.numberDensities` 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:
+            dLLprev = (
+                self.material.linearExpansionPercent(Tc=self.temperatureInC) / 100.0
+            )
+        else:
+            dLLprev = 0.0
+        try:
+            vol = self.getVolume()
+        except:
+            # 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 = {}
         self.p.numberDensities.update(numberDensities)
-        # 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``.
-        self.p.assigned = parameters.SINCE_ANYTHING
-        self.p.paramDefs["numberDensities"].assigned = parameters.SINCE_ANYTHING
+
+        # check if thermal expansion changed
+        dLLnew = self.material.linearExpansionPercent(Tc=self.temperatureInC) / 100.0
+        if dLLprev != dLLnew and dLLprev != 0.0:
+            # the thermal expansion changed so the volume change is happening at same time as
+            # density change was requested. Attempt to make mass consistent with old dims
+            # (since the density change was for the old volume and otherwise mass wouldn't be conserved)
+
+            # enable recalculation of volume, otherwise it uses stored.
+            self.clearLinkedCache()
+            if vol is not None:
+                factor = vol / self.getVolume()
+            else:
+                factor = area / self.getArea()
+            self.changeNDensByFactor(factor)
+
+    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
+        self._changeOtherDensParamsByFactor(factor)
 
     def getEnrichment(self):
         """Get the mass enrichment of this component, as defined by the material."""
@@ -1247,7 +1307,9 @@ def getIntegratedMgFlux(self, adjoint=False, gamma=False):
             if not self.parent:
                 return np.zeros(1)
 
-            volumeFraction = (self.getVolume() / self.parent.getSymmetryFactor() ) / self.parent.getVolume()
+            volumeFraction = (
+                self.getVolume() / self.parent.getSymmetryFactor()
+            ) / self.parent.getVolume()
             return volumeFraction * self.parent.getIntegratedMgFlux(adjoint, gamma)
 
         # pin-level flux is available. Note that it is NOT integrated on the param level.
@@ -1262,7 +1324,11 @@ def getIntegratedMgFlux(self, adjoint=False, gamma=False):
             else:
                 pinFluxes = self.parent.p.pinMgFluxes
 
-        return pinFluxes[self.p.pinNum - 1] * self.getVolume() / self.parent.getSymmetryFactor()
+        return (
+            pinFluxes[self.p.pinNum - 1]
+            * self.getVolume()
+            / self.parent.getSymmetryFactor()
+        )
 
     def getPinMgFluxes(
         self, adjoint: Optional[bool] = False, gamma: Optional[bool] = False

armi/reactor/composites.py

@@ -1575,6 +1575,10 @@ def changeNDensByFactor(self, factor):
             nuc: val * factor for nuc, val in self.getNumberDensities().items()
         }
         self.setNumberDensities(densitiesScaled)
+        self._changeOtherDensParamsByFactor(factor)
+
+    def _changeOtherDensParamsByFactor(self, factor):
+        """Change the other density storage params by factor."""
         # Update detailedNDens
         if self.p.detailedNDens is not None:
             self.p.detailedNDens *= factor

armi/reactor/converters/geometryConverters.py

@@ -1568,6 +1568,7 @@ def scaleParamsRelatedToSymmetry(core, paramsToScaleSubset=None):
             for b, bSymmetric in zip(a, aSymmetric):
                 _scaleParamsInBlock(b, bSymmetric, completeListOfParamsToScale)
 
+
 def _generateListOfParamsToScale(core, paramsToScaleSubset):
     fluxParamsToScale = (
         core.getFirstBlock()