`expandTo` in `nuclide flags` produces different number densities than using elemental nuclide

[keckler]

I noticed that the total atom density summed over all isotopes of a given element will be different in the ARMI reactor model depending on if one uses the expandTo attribute of the nuclide flags in the blueprints. I have confirmed this at least for the case of oxygen, though I have not tried any other elements.

I have attached two very simple examples that contain only a single block with only oxygen in it. The only difference is that the first case (named "combined") uses elemental oxygen, while the second case (named "expanded") uses the expandTo option to break the oxygen into O16 and O17. I would expect that both cases would produce the exact same number density when summed over all nuclides, but they do not.

import armi
armi.configure()

db_combined = databaseFactory("anl-afci-177-combined.h5", 'r')

db_expanded = databaseFactory("anl-afci-177-expanded.h5", 'r')

with db_combined:
     r_combined = db_combined.load(0,0)

with db_expanded:
     r_expanded = db_expanded.load(0,0)

sum(r_combined.core[0].getNumberDensities().values())
1.0002744781201818

sum(r_expanded.core[0].getNumberDensities().values())
1.0002505665856218

While the difference is small, it is still not expected, and it seems like a bug to me.

expandTo bug.zip

[john-science]

I can't prove it yet, but my guess is that (as in your example above) the normal/combined case will always yield a higher number than the "expandTo" case.

What I expect to see somewhere in the code is...

1. User selects expandTo: [O16, O17]
2. The other Oxygen nuclides are dropped.
3. The total mass of each of O16 and O17 are normalized to the total mass of all types of Oxygen.

I can't find (3) anywhere in the code. So, IF that doesn't exist, that would be the problem.

armi/armi/reactor/blueprints/isotopicOptions.py

Line 112 in 9dc6760

expandTo = yamlize.Attribute(type=yamlize.StrList, default=None)

armi/armi/reactor/blueprints/__init__.py

Lines 416 to 423 in 9dc6760

	if (
	elemental.name in nuclideFlags
	and nuclideFlags[elemental.element.symbol].expandTo
	):
	# user-input expandTo has precedence
	newNuclides = [
	nuclideBases.byName[nn]
	for nn in nuclideFlags[elemental.element.symbol].expandTo

armi/armi/reactor/blueprints/componentBlueprint.py

Lines 381 to 385 in 9dc6760

	nuclidesToBecome = (
	[nuclideBases.byName[nn] for nn in nucFlags.expandTo]
	if (nucFlags and nucFlags.expandTo)
	else None
	)

I also do not see a unit test checking the math on this. I would expect it to be here.

[john-science]

@keckler Okay, this seems like the easiest/shortest test I could devise to prove this bug exists: 84001de

armi/armi/reactor/blueprints/tests/test_customIsotopics.py

Lines 317 to 339 in 84001de

	def test_expandTo(self):
	"""A quick, end-to-end test that using "expandTo" in blueprints works as intended."""
	# build a new Assembly, with only C12 instead of C
	cs = settings.Settings()
	cs = cs.modified(newSettings={CONF_XS_KERNEL: "MC2v2"})
	newYaml = self.yamlString.replace(
	"C: {burn: false, xs: true}",
	'C: {burn: false, xs: true, expandTo: ["C12"]}',
	)
	self.assertNotEqual(self.yamlString, newYaml)
	bp = blueprints.Blueprints.load(newYaml)
	newAssem = bp.constructAssem(cs, name="fuel a")
	# prove the correct nuclides are in each assembly
	self.assertIn("C", self.a.getNumberDensities())
	self.assertNotIn("C12", self.a.getNumberDensities())
	self.assertIn("C12", newAssem.getNumberDensities())
	self.assertNotIn("C", newAssem.getNumberDensities())
	# prove that the total mass was not changed
	oldSum = sum(self.a.getNumberDensities().values())
	newSum = sum(newAssem.getNumberDensities().values())
	self.assertEqual(oldSum, newSum)

This should pass, but it fails with:

AssertionError: 0.06289884851817294 != 0.06290986051086071

[keckler]

Yeesh, it fails even for Carbon, that is damning!

[john-science]

Okay, @keckler , I still don't have a solution. But I used your "Oxygen-only block" example to inform 3 unit tests: ff543c3

• expandTo O16 - passes
• expandTo O17 - fails, underpredicts
• expandTo O16 and O17 - fails, underpredicts

This should probably be telling, but so far, no joy.

[john-science]

@keckler I talked to @jakehader , and I am not sure there should be a warning here. I mean, the feature exists, so it's not "wrong". But what do you think?

My current feeling is that we should beef up some docstrings, to explain to people how potentially strange and non-physical using this nuclide flag can be, but then we close this ticket and focus on the redesign ticket Jake opened over here: #320

I think the re-design is our endgame here.

I'm just thinking aloud though. Thoughts?

[keckler]

Can somebody give me a good example/reason why it would be a good idea for a user to just zero out a real isotope from all materials in a run? At this point, I cannot understand how that is either (1) reasonable nor (2) safe.

I honestly don't know the use case for this. My only thought is someone might ditch a particular isotope from an element because it isn't present in their XS library. But I have multiple confusions about that:

1. Why shouldn't that be handled by the specific plugin that has the problem?
2. What happens when your depletion solver gives you some of that isotope from the depletion process? Do you just throw it away? Probably not? But then you end up with the problem-isotope anyways, so why did you bother tossing it from the BOL reactor in the first place?
3. Even if you don't have a cross-section for a particular isotope of an element (let's say O18, which makes up 0.2% of naturally abundant oxygen), how do you justify including MORE of a different isotope in its place? A different isotope has different characteristics. So you maybe replaced an isotope that was neutronically unimportant (hence its lack of inclusion in the XS lib) with one that has more neutronic impact. That does not sound reasonable to me without some kind of justification.

[john-science]

@keckler My notion was not that I find this feature safe, or comforting. But that I know people are using this in the real world. And Jake already has the "redesign this feature" ticket open: #320

But to answer your above questions more directly, I can only imagine the problem is common in nuclear engineering that you lack good cross section data sometimes, and dealing with that problem will always be awkward and semi-non-physical. Yuck.

1. Yeah, this feels like the kind of hack that was built for a single downstream project, let's find that project and move this there.
2. How does this affect the burn chain? I can't imagine an answer that would make either of us happy.
3. Agreed. This will be non-physical.