Allowing XS plotting module to handle excitation reactions

[eitan-weinstein]

Closes #303.

This PR implements a number of changes in service of the goal of allowing xs_plotting.extract_continuous_data() to accurately produce cross-sections for specific excitation reactions. Similarly to the approach used #229, where we had to reference the explicit excitation pathways from MF9/10, here instead of just using the 'LFS' tag, to extract the cross-sections, each pathway's subsection data is needed. For excitation reactions specified in MF10, the update is fairly simple, as MF10 already contains cross-section data. For MF9, however, the energy dependent variable contained within is the reaction multiplicity, a fraction ranging from 0-1 corresponding the proportion of the total cross-section from MF3 (which contains all excitation pathways) going to each specific pathway. As such, the formula to calculate an energy-dependent cross-section for an MF9 pathway-specific reaction is:

$\sigma_{\mathrm{excitation}}(E) = \sigma_{\mathrm{MF3,total}}(E)*multiplicity(E)$

In order for these multiplicities to be useful in reconstructing pathway-specific cross-sections, however, they must first be interpolated according to the interpolation scheme(s) designated in each subsection header by the 'INT' flag. Moreover, the specific energies/number of energies contained in MF9 vs MF3 for a given reaction may not match, so the MF3 cross-sections also need to be interpolated with the MF9 energies.

Built into NJOY's endf utility module is a function terpa(), which is used to interpolate ENDF TAB1 (records for "one-dimensional tabulated functions such as y(x)", as defined in Section 0.6.4.7 of the ENDF-6 manual) according the TAB1-specified interpolation scheme(s). In order to make use of this existing functionality and avoid re-inventing the wheel on this issue, this PR also includes a small Fortran wrapper of terpa(), njoy_endf_wrapper.f90. I've incorporated numpy.f2py to interface with this Fortran code to make it callable within the MF9 handling block of xs_plotting.extract_continuous_data() so that an array of pathway-specific cross-sections can be appropriately produced at each excitation.

[gonuke]

I'd prefer to avoid adding Fortran to this

[gonuke]

Has no-one written a python version of this yet??

[gonuke]

FWIW, I think it would be simple to write a python version for TAB1 interpolation

[eitan-weinstein]

Has no-one written a python version of this yet??

It looks like I may actually be able to utilize openmc.data to accomplish this, so long as we're fine further depending on OpenMC

[eitan-weinstein]

It looks like I may actually be able to utilize openmc.data to accomplish this, so long as we're fine further depending on OpenMC

As a follow up, I know that I just migrated everything to endf_parserpy (though that was not so onerous), but now I wonder whether all of our ENDF interfacing could just be done through OpenMC anyways. I'm not sure if OpenMC's ENDF capabilities include the ability to parse PENDF files, though, so that could be the limiting factor for total migration.

Either way, for the specific application of the continuous-energy parsing needed for this PR, I think OpenMC should be able to accomplish what we need.

[gonuke]

A few suggestions here...

[gonuke]

This suggests that if pathways is at least as long as isomeric_state that it will be filled with data that is in the order of the isomeric states...

[gonuke]

Also, it seems that the only way you index pathways is by an ordinal value, so couldn't/shouldn't this be a vector instead of a dictionary?

[eitan-weinstein]

To the first point, I believe that this is the correct implementation, based on my previous work in handling how TENDL formats isomeric states, in that a nuclide's "m" and "n" isomers may correspond to excitations not equal to 1 and 2 respectively. For example, coming across a reaction with daughters in [0,4,29], these would correspond to the ground state (0 will always be present if multiple pathways are given), the first ("m") excited state and the second ("n") excited state.

To your point on vectorization, however, that's well-taken; there's no justifiable reason why this would need to be a dictionary based on my current usage.

[gonuke]

Thanks for that clarification on the difference between excitation state & isomeric flag - that makes sense

[gonuke]

More thoughts as I understand/infer some of the data patterns

[gonuke]

Wait... these excited states are not in the same order as the isomeric flags?

Then maybe the dictionary is simpler after all...

I think this is a place where a comment stating the assumptions/expectations would be warranted:

• isomeric states are flagged in the order of their excited state value, even if the excited states don't appear in that order

[eitan-weinstein]

rxn.products maintains the ordering of the original file, so in principle they should always be in ascending order by ENDF conventions, but I put this there purely as a safety mechanism. Mostly, I was concerned about the possibility of mismatched indexing if we're reading a deprecated MF 9/10 section that places them out of order. This is not an issue I've encountered, but given that TENDL already exhibits some weird unconventional formatting within MF 9/10 (i.e. any MT appearing in one should necessarily appear in the other, but exclusively does not), I figured guaranteeing an ascending order would probably be prudent.

[gonuke]

I agree that ascending order makes sense in principle, but since we are relying on this order to be consistent with the labeling of isomers, it makes me nervous that it's not very robust.

[eitan-weinstein]

Given that, if this is not a reliable safeguard, should we just stick to relying on the listed order to avoid further propagating any mislabeling issues? For that, we could simplify the whole block down to something like this:

pathways = [product for product in rxn.products if product not in {'neutron', 'photon', 'electron'}]

If we're not sorting or keying based on excitation level, but rather just preserving the index, your original array suggestion may be best suited after all.

[gonuke]

I guess I don't know what the correct way to do it is, and whether it's documented anywhere?

[eitan-weinstein]

From what I can find in the ENDF-6 manual, there is no explicit ordering of LFS states listed in MF9/10, but for MF8, the subsection structure is described as follows:

The structure of each section always starts with a HEAD record and ends with a SEND record. Subsections contain data for particular final state of the reaction product (LFS). The number of subsections NS is given on the HEAD record for the section. The subsections are ordered by increasing value of LFS... (ENDF6 Manual, Section 8.2.1)

Why this is not stated explicitly in the chapters for MF9/10 is not clear, though I've only ever encountered the subsections to be increasing order of LFS. That said, this is the specific language describing subsection structuring in MF9 is as follows:

The sections are ordered by increasing MT number. Within a section for a given MT are subsections for different final states of the daughter product (LFS). (ENDF6 Manual, Section 9.1)

While there is no mention of the ordering of the subsections, below, in the "Formats" section, the LFS quantity descriptor is:

LFS Indicator to specify the level number of the nuclide (ZAP) (as defined in File 8) produced in the reaction (MT number).
LFS = 0: the final state is the ground state.
LFS = 1: the final state is the first excited state.
LFS = 2: the final state is the second excited state.
———-
———-
LFS = 98: an unspecified range of final states. (ENDF6 Manual, Section 9.2)

Again, this is not an explicit assertion of the ordering, but it seems that the ascending order in the descriptor implies a like-order of the subsections by LFS. Of course, relying on a perceived implication (even if it is consistent with the structures I've encountered in my wading through TENDL files) may not be a strong enough basis by which to make a determination of appropriateness for any ordered subsection unpacking solution.

I'll keep digging around to see if I can find other documentation in support of or countering this notion, though I'm skeptical that it would exist if not in the authoritative ENDF6 manual.