terrapower / armi

An open-source nuclear reactor analysis automation framework that helps design teams increase efficiency and quality
https://terrapower.github.io/armi/
Apache License 2.0
231 stars 90 forks source link

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

Open keckler opened 2 months ago

keckler commented 2 months ago

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 commented 2 months ago

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.

https://github.com/terrapower/armi/blob/9dc6760de7da3800cbe98e83066cf5f343e81a6a/armi/reactor/blueprints/isotopicOptions.py#L112

https://github.com/terrapower/armi/blob/9dc6760de7da3800cbe98e83066cf5f343e81a6a/armi/reactor/blueprints/__init__.py#L416-L423

https://github.com/terrapower/armi/blob/9dc6760de7da3800cbe98e83066cf5f343e81a6a/armi/reactor/blueprints/componentBlueprint.py#L381-L385

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

john-science commented 2 months ago

@keckler Okay, this seems like the easiest/shortest test I could devise to prove this bug exists: https://github.com/terrapower/armi/commit/84001de02f0c5e17b533077a78fde9715813f69d

https://github.com/terrapower/armi/blob/84001de02f0c5e17b533077a78fde9715813f69d/armi/reactor/blueprints/tests/test_customIsotopics.py#L317-L339

This should pass, but it fails with:

AssertionError: 0.06289884851817294 != 0.06290986051086071
keckler commented 2 months ago

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

john-science commented 2 months ago

Okay, @keckler , I still don't have a solution. But I used your "Oxygen-only block" example to inform 3 unit tests: https://github.com/terrapower/armi/commit/ff543c3f3d47b492dd13a007002b97f0ffac1264

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

john-science commented 2 months ago

@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: https://github.com/terrapower/armi/issues/320

I think the re-design is our endgame here.

I'm just thinking aloud though. Thoughts?

keckler commented 2 months ago

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 commented 2 months ago

@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: https://github.com/terrapower/armi/issues/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.