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
236 stars 89 forks source link

`Block.getReactionRates()` doesn't work for certain isotopes #808

Closed keckler closed 2 years ago

keckler commented 2 years ago

The Block.getReactionRates() method supposedly will allow one to specify an isotope and get the rate of a specified reaction on the given block: https://github.com/terrapower/armi/blob/69c778934b9e0cfe3751f2efde836126889a1902/armi/reactor/blocks.py#L1489-L1520

If the user does not specify the nDensity argument, Block.getReactionRates() will attempt to look up the isotope's number density in the block using the Block.getNumberDensity() method, which under the hood iteratively calls the Component.getNumberDensity() method: https://github.com/terrapower/armi/blob/69c778934b9e0cfe3751f2efde836126889a1902/armi/reactor/components/component.py#L596-L610

Of note is that this Component.getNumberDensity() method will return 0.0 if the nucName is not found in the dict.

This is all fine and dandy.

So if you had a block blockInstance, had a properly-loaded ISOTXS file, and called blockInstance.getReactionRates("U235"), you would get the result you're looking for:

In [13]: r.core.getAssemblyWithStringLocation('001-001')[6].getReactionRates('U238')['nF']
Out[13]: 5134734346273870.0

BUT, if you instead called the method for a certain subset of isotopes, you might have problems:

In [14]: r.core.getAssemblyWithStringLocation('001-001')[6].getReactionRates('PU239')['nF']
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-14-0590e0975867> in <module>
----> 1 r.core.getAssemblyWithStringLocation('001-001')[6].getReactionRates('PU239')['nF']

c:\users\ckeckler\codes\nala\framework\armi\reactor\blocks.py in getReactionRates(self, nucName, nDensity)
   1517                 self.p.xsType,
   1518                 self.getIntegratedMgFlux(),
-> 1519                 nDensity,
   1520             )
   1521         except AttributeError:

c:\users\ckeckler\codes\nala\framework\armi\reactor\components\component.py in getReactionRateDict(nucName, lib, xsType, mgFlux, nDens)
   1239     """
   1240     key = "{}{}A".format(nucName, xsType)
-> 1241     libNuc = lib[key]
   1242     rxnRates = {"n3n": 0}
   1243     for rxName, mgXSs in [

c:\users\ckeckler\codes\nala\framework\armi\nuclearDataIO\xsLibraries.py in __getitem__(self, key)
    398
    399     def __getitem__(self, key):
--> 400         return self._nuclides[key]
    401
    402     def get(self, nuclideLabel, default):

KeyError: 'PU239AA'

The underlying issue here ends up being that PU239 in the ISOTXS file is actually shortened to just PU39. And a variety of other isotopes are like this as well, so that the isotope representations are all only 4 characters long. So basically the any element with a 2-character symbol that has more than 99 nucleons.

And then if you think you're clever and you instead supply this shorter isotope representation to Block.getReactionRates(), you get this inaccurate result:

In [16]: r.core.getAssemblyWithStringLocation('001-001')[6].getReactionRates('PU39')['nF']
Out[16]: 0.0

But it is NOT that there is zero PU239 in the block:

In [17]: r.core.getAssemblyWithStringLocation('001-001')[6].getNumberDensity('PU239')
Out[17]: 0.0014610752712455898

Rather, 0.0 is returned because Component.getNumberDensity() doesn't know how to handle the shorter isotope representation, so it just returns 0.0 because it doesn't know what PU39 is:

In [18]: r.core.getAssemblyWithStringLocation('001-001')[6].getNumberDensity('PU39')
Out[18]: 0.0

So basically, as is currently, you CANNOT get the reaction rates for quite a few very important isotopes using this method.

I have no idea how this went undiscovered for so long. Possibly because the method is more typically used for generating cross-sections, as alluded to in it's docstring? In that case, you won't run into the same problem because the number density is never looked up.

keckler commented 2 years ago

Happily, I think this can be fixed rather easily via some utilities in the nuclideBases module: https://terrapower.github.io/armi/.apidocs/armi.nucDirectory.nuclideBases.html#module-armi.nucDirectory.nuclideBases

I can go ahead and get a PR for this.