Repository: terrapower/armi
Pull request #1822 (2024-08-14 15:20:05)

Fixing the affect of custom isotopics on material density
## What is the change?

This fixes some problems that were found with #1745. Those issues are described in #1818 .

## Why is the change being made?

We are fixing a bug that made it into the codebase.

---

## Checklist

<!--
    You (the pull requester) should put an `x` in the boxes below you have completed.

    (If a checkbox doesn't apply to your PR, check it anyway.)
-->

- [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.

<!-- Check the code quality -->

- [x] The code style follows [good practices](https://terrapower.github.io/armi/developer/standards_and_practices.html).
- [x] The commit message(s) follow [good practices](https://terrapower.github.io/armi/developer/tooling.html).

<!-- Check the project-level cruft -->

- [ ] The [release notes](https://terrapower.github.io/armi/developer/tooling.html#add-release-notes) have been updated if necessary.
- [x] The [documentation](https://terrapower.github.io/armi/developer/tooling.html#document-it) is still up-to-date in the `doc` folder.
- [x] The dependencies are still up-to-date in `pyproject.toml`.

john-science: @jakehader   Just a reminder of the ARMI PR policy:

https://github.com/terrapower/armi/blob/bd0b43adadf6e9c0063e27e490353f4a16e1b623/doc/developer/tooling.rst?plain=1#L45-L49

We are all extremely busy, and if a PR gets opened before it is ready I will get an email:

* every single time you commit, and the CI passes
* for each and every time a single CI job fails

I love the help, and the interesting in fixing the feature. But in the future, please wait until the feature is ready for final review before opening a PR.
jakehader: > @jakehader Just a reminder of the ARMI PR policy:
> 
> https://github.com/terrapower/armi/blob/bd0b43adadf6e9c0063e27e490353f4a16e1b623/doc/developer/tooling.rst?plain=1#L45-L49
> 
> We are all extremely busy, and if a PR gets opened before it is ready I will get an email:
> 
> * every single time you commit, and the CI passes
> * for each and every time a single CI job fails
> 
> I love the help, and the interesting in fixing the feature. But in the future, please wait until the feature is ready for final review before opening a PR.

Thanks for the reminder @john-science - I put the PR here but forgot to make it a draft. This was just intended to provide some ideas for how to resolve the ticket and to put something out there to be stress tested by @keckler and others. 

There are failing unit tests because I changed the construct arguments so that needs to be resolved still (probably a poor design to require case settings to build a component). Also, we will need to add more unit tests around custom Isotopics (same comment on the ticket).
keckler: I ran some tests on my FFTF model to see if this branch works. A couple points:
1. The keff in my model increased by nearly 5000 pcm to a value that is very incorrect.
2. The predicted isothermal temperature reactivity coefficient looks correct.
3. The fuel density when the model is initialized to Thot is way too high (~10%).
4. The _change_ in density between different Thot values appears to be roughly correct, and the number of atoms is being correctly conserved.

So it seems like this change is half working. Somehow the initialized model is getting the densities very incorrect, though.
jakehader: Was this fixed elsewhere @keckler?
keckler: Nope, still broken
keckler: Running my FFTF tests, I now believe that this PR is working correctly.
keckler: @john-science I believe that this is ready for your eyes. 

In addition to the updated unit tests, I did quite a bit of offline testing with my FFTF model to convince myself that this is working properly.
john-science: Random thought: I think this PR is big enough for a release note in `doc/release/0.5.rst`.

Patch:
diff --git a/armi/reactor/blueprints/blockBlueprint.py b/armi/reactor/blueprints/blockBlueprint.py
index 692ae84a..a1e7a554 100644
--- a/armi/reactor/blueprints/blockBlueprint.py
+++ b/armi/reactor/blueprints/blockBlueprint.py
@@ -28,6 +28,7 @@ from armi.reactor.components.component import Component
 from armi.reactor.composites import Composite
 from armi.reactor.converters import blockConverters
 from armi.reactor.flags import Flags
+from armi.settings.fwSettings.globalSettings import CONF_INPUT_HEIGHTS_HOT
 
 
 def _configureGeomOptions():
@@ -144,7 +145,11 @@ class BlockBlueprint(yamlize.KeyedList):
             filteredMaterialInput, byComponentMatModKeys = self._filterMaterialInput(
                 materialInput, componentDesign
             )
-            c = componentDesign.construct(blueprint, filteredMaterialInput)
+            c = componentDesign.construct(
+                blueprint,
+                filteredMaterialInput,
+                cs[CONF_INPUT_HEIGHTS_HOT],
+            )
             components[c.name] = c
 
             # check that the mat mods for this component are valid options
diff --git a/armi/reactor/blueprints/componentBlueprint.py b/armi/reactor/blueprints/componentBlueprint.py
index 41dcc0dc..289a0de8 100644
--- a/armi/reactor/blueprints/componentBlueprint.py
+++ b/armi/reactor/blueprints/componentBlueprint.py
@@ -186,7 +186,7 @@ class ComponentBlueprint(yamlize.Object):
     mergeWith = yamlize.Attribute(type=str, default=None)
     area = yamlize.Attribute(type=float, default=None)
 
-    def construct(self, blueprint, matMods):
+    def construct(self, blueprint, matMods, inputHeightsConsideredHot):
         """Construct a component or group.
 
         .. impl:: User-defined on material alterations are applied here.
@@ -205,6 +205,17 @@ class ComponentBlueprint(yamlize.Object):
             The ``applyInputParams()`` method of that material class is then called,
             passing in the associated material modifications data, which the material
             class can then use to modify the isotopics as necessary.
+
+        Parameters
+        ----------
+        blueprint : Blueprints
+            Blueprints object containing various detailed information, such as nuclides to model
+
+        matMods : dict
+            Material modifications to apply to the component.
+
+        inputHeightsConsideredHot : bool
+            See the case setting of the same name.
         """
         runLog.debug("Constructing component {}".format(self.name))
         kwargs = self._conformKwargs(blueprint, matMods)
@@ -214,7 +225,9 @@ class ComponentBlueprint(yamlize.Object):
             constructedObject = composites.Composite(self.name)
             for groupedComponent in group:
                 componentDesign = blueprint.componentDesigns[groupedComponent.name]
-                component = componentDesign.construct(blueprint, matMods=dict())
+                component = componentDesign.construct(
+                    blueprint, {}, inputHeightsConsideredHot
+                )
                 # override free component multiplicity if it's set based on the group definition
                 component.setDimension("mult", groupedComponent.mult)
                 _setComponentFlags(component, self.flags, blueprint)
@@ -229,29 +242,33 @@ class ComponentBlueprint(yamlize.Object):
                 constructedObject.material.getTD()
             )
 
-        # set the custom density for non-custom material components after construction
-        self.setCustomDensity(constructedObject, blueprint, matMods)
+        self._setComponentCustomDensity(
+            constructedObject,
+            blueprint,
+            matMods,
+            inputHeightsConsideredHot,
+        )
 
         return constructedObject
 
-    def setCustomDensity(self, constructedComponent, blueprint, matMods):
+    def _setComponentCustomDensity(
+        self, comp, blueprint, matMods, inputHeightsConsideredHot
+    ):
         """Apply a custom density to a material with custom isotopics but not a 'custom material'."""
         if self.isotopics is None:
             # No custom isotopics specified
             return
 
-        density = blueprint.customIsotopics[self.isotopics].density
-        if density is None:
+        densityFromCustomIsotopic = blueprint.customIsotopics[self.isotopics].density
+        if densityFromCustomIsotopic is None:
             # Nothing to do
             return
 
-        if density <= 0:
+        if densityFromCustomIsotopic <= 0:
             runLog.error(
                 "A zero or negative density was specified in a custom isotopics input. "
                 "This is not permitted, if a 0 density material is needed, use 'Void'. "
-                "The component is {} and the isotopics entry is {}.".format(
-                    constructedComponent, self.isotopics
-                )
+                f"The component is {comp} and the isotopics entry is {self.isotopics}."
             )
             raise ValueError(
                 "A zero or negative density was specified in the custom isotopics for a component"
@@ -261,30 +278,40 @@ class ComponentBlueprint(yamlize.Object):
         if not isinstance(mat, materials.Custom):
             # check for some problem cases
             if "TD_frac" in matMods.keys():
-                runLog.warning(
-                    "Both TD_frac and a custom density (custom isotopics) has been specified for "
-                    "material {}. The custom density will override the density calculated using "
-                    "TD_frac.".format(self.material)
+                runLog.error(
+                    f"Both TD_frac and a custom isotopic with density {blueprint.customIsotopics[self.isotopics]} "
+                    f"has been specified for material {self.material}. This is an overspecification."
                 )
             if not mat.density(Tc=self.Tinput) > 0:
                 runLog.error(
-                    "A custom density has been assigned to material '{}', which has no baseline "
+                    f"A custom density has been assigned to material '{self.material}', which has no baseline "
                     "density. Only materials with a starting density may be assigned a density. "
-                    "This comes up e.g. if isotopics are assigned to 'Void'.".format(
-                        self.material
-                    )
+                    "This comes up e.g. if isotopics are assigned to 'Void'."
                 )
                 raise ValueError(
                     "Cannot apply custom densities to materials without density."
                 )
 
-            densityRatio = density / constructedComponent.density()
-            constructedComponent.changeNDensByFactor(densityRatio)
+            # Apply a density scaling to account for the temperature change between Tinput
+            # Thot. There may be a better place in the initialization to determine
+            # if the block height will be interpreted as hot dimensions, which would
+            # allow us to not have to pass the case settings down this far
+            dLL = comp.material.linearExpansionFactor(
+                Tc=comp.temperatureInC, T0=comp.inputTemperatureInC
+            )
+            if inputHeightsConsideredHot:
+                f = 1.0 / (1 + dLL) ** 2
+            else:
+                f = 1.0 / (1 + dLL) ** 3
+
+            scaledDensity = comp.density() / f
+            densityRatio = densityFromCustomIsotopic / scaledDensity
+            comp.changeNDensByFactor(densityRatio)
 
             runLog.important(
                 "A custom material density was specified in the custom isotopics for non-custom "
-                "material {}. The component density has been altered to "
-                "{}.".format(mat, constructedComponent.density()),
+                f"material {mat}. The component density has been altered to "
+                f"{comp.density()} at temperature {comp.temperatureInC} C",
                 single=True,
             )
 
diff --git a/armi/reactor/blueprints/isotopicOptions.py b/armi/reactor/blueprints/isotopicOptions.py
index 37ded762..e8b2b354 100644
--- a/armi/reactor/blueprints/isotopicOptions.py
+++ b/armi/reactor/blueprints/isotopicOptions.py
@@ -395,10 +395,10 @@ class CustomIsotopic(yamlize.Map):
         if self.density is not None:
             if not isinstance(material, materials.Custom):
                 runLog.important(
-                    "A custom density or number densities has been specified for non-custom "
-                    "material {}. The material object's density will not be updated to prevent unintentional "
-                    "density changes across the model. Only custom materials may have a density "
-                    "specified.".format(material),
+                    "A custom isotopic with associated density has been specified for non-`Custom` "
+                    f"material {material}. The reference density of materials in the materials library "
+                    "will not be changed, but the associated components will use the density "
+                    "implied by the custom isotopics.",
                     single=True,
                 )
                 # specifically, non-Custom materials only use refDensity and dLL, mat.customDensity has no effect
diff --git a/doc/release/0.5.rst b/doc/release/0.5.rst
index 0b92c954..f730979d 100644
--- a/doc/release/0.5.rst
+++ b/doc/release/0.5.rst
@@ -98,6 +98,7 @@ Bug Fixes
 #. Fixed hex block rotation in ``plotBlockDiagram``. (`PR#1926 <https://github.com/terrapower/armi/pull/1926>`_)
 #. Fixed edge case in ``assemblyBlueprint._checkParamConsistency()``. (`PR#1928 <https://github.com/terrapower/armi/pull/1928>`_)
 #. Fixed wetted perimeter for hex inner ducts. (`PR#1985 <https://github.com/terrapower/armi/pull/1985>`_)
+#. Fixing number densities when custom isotopics and material properties are combined. (`PR#1822 <https://github.com/terrapower/armi/pull/1822>`_)
 
 Quality Work
 ------------


Test patch:
diff --git a/armi/reactor/blueprints/tests/test_blueprints.py b/armi/reactor/blueprints/tests/test_blueprints.py
index d63c1246..7827dd5e 100644
--- a/armi/reactor/blueprints/tests/test_blueprints.py
+++ b/armi/reactor/blueprints/tests/test_blueprints.py
@@ -32,6 +32,7 @@ from armi.tests import TEST_ROOT
 from armi.utils import directoryChangers
 from armi.utils import textProcessors
 from armi.reactor.blueprints.gridBlueprint import saveToStream
+from armi.settings.fwSettings.globalSettings import CONF_INPUT_HEIGHTS_HOT
 
 
 class TestBlueprints(unittest.TestCase):
@@ -612,7 +613,7 @@ assemblies:
         # which is required during construction of a component
         design._resolveNuclides(cs)
         componentDesign = design.componentDesigns["freefuel"]
-        topComponent = componentDesign.construct(design, matMods=dict())
+        topComponent = componentDesign.construct(design, {}, cs[CONF_INPUT_HEIGHTS_HOT])
         self.assertEqual(topComponent.getDimension("od", cold=True), 4.0)
         self.assertGreater(topComponent.getVolume(), 0.0)
         self.assertGreater(topComponent.getMass("U235"), 0.0)
diff --git a/armi/reactor/blueprints/tests/test_customIsotopics.py b/armi/reactor/blueprints/tests/test_customIsotopics.py
index f20cec81..d971cb65 100644
--- a/armi/reactor/blueprints/tests/test_customIsotopics.py
+++ b/armi/reactor/blueprints/tests/test_customIsotopics.py
@@ -309,7 +309,12 @@ assemblies:
     @classmethod
     def setUpClass(cls):
         cs = settings.Settings()
-        cs = cs.modified(newSettings={CONF_XS_KERNEL: "MC2v2"})
+        cs = cs.modified(
+            newSettings={
+                CONF_XS_KERNEL: "MC2v2",
+                "inputHeightsConsideredHot": False,
+            }
+        )
 
         cls.bp = blueprints.Blueprints.load(cls.yamlString)
         cls.a = cls.bp.constructAssem(cs, name="fuel a")
@@ -356,14 +361,19 @@ assemblies:
         # A block with custom density set via number density
         fuel8 = self.a[8].getComponent(Flags.FUEL)
 
+        dLL = fuel2.material.linearExpansionFactor(Tc=600, T0=25)
+        # the exponent here is 3 because inputHeightsConsideredHot = False.
+        # if inputHeightsConsideredHot were True, then we would use a factor of 2 instead
+        f = 1 / ((1 + dLL) ** 3)
+
         # Check that the density is set correctly on the custom density block,
         # and that it is not the same as the original
-        self.assertAlmostEqual(19.1, fuel2.density())
+        self.assertAlmostEqual(19.1 * f, fuel2.density())
         self.assertNotAlmostEqual(fuel0.density(), fuel2.density(), places=2)
         # Check that the custom density block has the correct material
         self.assertEqual("UZr", fuel2.material.name)
         # Check that the block with only number densities set has a new density
-        self.assertAlmostEqual(19.1, fuel8.density())
+        self.assertAlmostEqual(19.1 * f, fuel8.density())
         # original material density should not be changed after setting a custom density component,
         # so a new block without custom isotopics and density should have the same density as the original
         self.assertAlmostEqual(fuel6.density(), fuel0.density())
@@ -387,12 +397,16 @@ assemblies:
 
             # Check for log messages
             streamVal = mockLog.getStdout()
-            self.assertIn("Both TD_frac and a custom density", streamVal, msg=streamVal)
+            self.assertIn(
+                "Both TD_frac and a custom isotopic with density",
+                streamVal,
+                msg=streamVal,
+            )
             self.assertIn(
                 "A custom material density was specified", streamVal, msg=streamVal
             )
             self.assertIn(
-                "A custom density or number densities has been specified",
+                "A custom isotopic with associated density has been specified for non-`Custom`",
                 streamVal,
                 msg=streamVal,
             )