ENDF support in openmc.data and resonance reconstruction

[paulromano]

This pull request extends the openmc.data module with from_endf() methods for many of the classes. They work the same as the from_ace() methods, so for example, you could do something like:

u235 = openmc.data.IncidentNeutron.from_endf('n-092_U_235.endf')
elastic_xs = u235[2].xs['0K']
nu_prompt = u235[18].products[0].yield_

If you're familiar with the ACE machinery, the extensions should seem natural. Resonance data is stored in a new resonances attribute on IncidentNeutron, and there are a series of classes for dealing with the different resonance formalisms. Notably, the API is capable of reconstructing cross sections from the resonance parameters for SLBW, MLBW, and Reich-Moore formalisms. The most immediate use I can think of for this is generating 0 K elastic scattering cross sections for heavy-nuclide resonance scattering treatments (DBRC and the like). Since reconstructing resonances in Python would be dirt slow, I've implemented it in Cython to get closer to C performance. The use of Cython is optional -- when you setup the package, if Cython is not found, the only consequence is that resonance reconstruction is not available. All other features are still available.

Resonant cross sections are represented using a new ResonancesWithBackground class. This class has a __call__() method just like any other function, so you can evaluate it at one or many energies of your choosing.

A few other notes:

• The FissionEnergyRelease class has been reworked to use the new ENDF functionality.
• TAB2 records in ENDF are not handled super cleanly right now. I'm kind of waiting for the GND folks to settle on how they plan to do low-level data containers for things like 1- and 2-D functions, and at that point, I'll try to move things to be more in line with GND.
• Not all photon energy/angle distributions are handled yet. Not a big deal since we don't use photon production data yet anyway.
• I haven't done anything with thermal scattering yet.

Notes on Resonance Reconstruction

The implementation is pretty slick if you ask me. Each formalism class (SingleLevelBreitWigner, ReichMoore, etc.) has a parameters attribute that is a Pandas dataframe of resonance parameters. The first time you go to reconstruct a cross section, the _prepare_resonances() method gets called which calculates shift and penetration factors and then stores them alongside the resonance parameters in a 2D array, one per (l,J) state. The actual reconstruction function in Cython can then efficiently use that 2D array.

Regarding accuracy of the reconstructed cross sections -- in general, it is really good. As one example, here's the relative difference in cross sections between openmc.data and NJOY2012 for U235 from ENDF/B-VII.1:

I've compared cross sections for all of ENDF/B-VII.1 (except for Cl35 and Pu239 for reasons I'll discuss) and almost everything is in very close agreement (errors < 1e-4) with NJOY2012. The cases where I saw relative differences larger than 1e-4 were:

• Ar40, Nb93, U240, Gd152, Dy160: These nuclides actually have negative elastic scattering cross sections at 0 K. NJOY sets them to 1e-8 b, which of course results in huge relative differences. The following plot shows elastic scattering in Gd152:
  
• Ce136, Sm150: The 1e-8 b minimum cross section also occurs even if the cross section is tiny but positive. This happens for these two nuclides. Here's a plot showing Ce136:
  
• Fe54: There's a huge resonance interference dip in elastic scattering and
  we miss the bottom by 0.1%.
• Pb206, Pb208: These are the only two nuclides with differences that I can't yet explain. The differences are in capture at low energy, which is 1/v. However, our reconstructed cross sections are just slightly jagged if you zoom in closely enough. My guess is we're losing precision somewhere, but I haven't been able to figure out where exactly yet. Both of these use the Reich-Moore approximation. Here's the relative difference in Pb208 capture:
  

For the very curious: a tarball with plots comparing all nuclides with resonances.

I should mention that in order to get the extraordinary agreement on most nuclides, I did make sure I was using the same fundamental physics constants as NJOY (which feed into the calculation of the neutron wave number). If I use updated CODATA values, the comparisons are not quite as perfect, but still very good.

There were two cases I couldn't run:

• Cl35: We can't handle LRF=7 (R-Matirx Limited) yet and this is the only guy in ENDF/B-VII.1 that uses that formalism.
• Pu239: In VII.1, there are actually three resolved resonance ranges. I need to look into whether it's safe to just combine the parameters from all three. For what it's worth, I tested the new B/VIII.0 beta2 file and it compares perfectly with NJOY.

tl;dr resonance reconstruction is not yet quite complete but is pretty damn good already.

Developer note

If, like me, you tend to add the OpenMC directory to your PYTHONPATH, you can still get the reconstruct Cython module built without having to install the entire API. To do this, just run:

python setup.py build_ext --inplace

[nelsonag]

👍 This is awesome

[cjosey]

So, I have been working with Pablo in rewriting OpenW in Julia (and fully eliminating WHOPPER dependencies so that it can be open sourced) over the last week. As such, this is all very fresh in my mind.

As far as I am aware, the level matrix form (I - K) to calculate U is less numerically stable than the X matrix form (SAMMY Revision 7 manual). I have never seen the level matrix form have pretty much any errors, but I wonder if that might help with Pb208. WHOPPER has a specific note for Pb208: "IOVLP > 1 is specially developed to handle the pseudo levels associated with Pb-208 which are strongly overlapping." The X matrix form uses wave function amplitudes (small gammas), which I find makes Reich-Moore easier to implement.

I'm not 100% sure if this is addressed or not (I searched for NLSC and did some spot checks), but check ENDF102 section 2.4.19. In summary, even if a spin group does not exist in the ENDF file, it still contributes to sigmaP. U-238 only has up to l=1 data, but you need l=2 to get the right answer (in GND, it is the attribute LvaluesNeededForConvergence, in ENDF it is NLSC). As far as I am aware, this also means you have to deal with separate channel spins even if they're not present. For U238, this is a 10^-7 error. For nuclides that end at very high E, it can be higher. Fudge says "screw it" and processes all nuclides to l = 3.

The code Pablo and I are working on already has partial LRF=7 support. I haven't implemented a numerically stable coulomb wave function yet (parts of one confluent hypergeometric function has to be analytically continued for l as an integer if written explicitly), but that's the only hurdle left. Unfortunately, I've moved to GND only, so I'm not sure if my code will be any use to you. The Cl-35 nuclide, although is LRF=7, is calculated using the Reich-Moore equations. It just has a different P for the (n, p) channel.

For Pu239, my understanding is that if the resonance regimes are separately defined like that, they're treated as quantum-mechanically separate as well. As such, to get the value at E = 1.1e3 eV, you use table 2 only. If this weren't the case, there would be a 0.16 barn discontinuity at 1.0e3 eV sigT, as Reich-Moore is strictly continuous (and File 3 here is not).

Finally, and I am not surprised by its absence at all as it gives me nightmares to think about, but are you planning on adding the Reich-Moore angular distribution reconstruction?

[paulromano]

@cjosey Thanks for the comments and suggestions. I've heard of the X matrix approach but haven't tried it before -- I'll take a closer look at the SAMMY manual.

I was able to solve the issues with Pb206 and Pb208. The problem was that if the real part of K is small and the phase shift is also small, you can end up losing precision when calculating the total cross section (1 - Re(U)) because Re(U) ends up being very very close to one. The solution (adopted from NJOY) is to expand the phase shift using Euler's formula and a Taylor series, keep the first two terms, and simplify. It very well may be that using the X matrix approach avoids this ad hoc fix in the first place.

For the "l values needed for convergence", you're correct that I don't handle that at present. However, I do handle contributions to the level matrix from channels which don't explicitly have resonances but which have l values that otherwise occur -- if you ignore this you do get very large errors. As far as U238, it looks like the ENDF file says NLSC=2, which I interpret to mean you need two l values (l=0 and l=1). Am I mistaken?

Handling multiple resolved ranges has been fixed now per your suggestion that they are quantum mechanically separate; comparisons on Pu239 look great.

No immediate plans for angular distribution reconstruction. I can name the equation that's used, but beyond that I'm completely ignorant on how it works. Are you guys thinking about adding that in OpenW?

[cjosey]

The X matrix approach might solve that issue, then. It yields a different equation instead of 1-Re[U] that should introduce smaller errors. However, I didn't know about the Taylor series path. I will have to look into that.

You are right in that I did misinterpret NLSC with an off-by-one error. It means l=0, l=1.

Angular distributions are quite painful, and I am unsure if we have long term plans to implement them. The reason being that we do not know if an angular pole representation can exist. If we suddenly find some way to do that, we'll need a reference solution to compare to, though.

Finally, if you are planning to add charged particles, I wrote a toy python script to implement P_l(rho, eta) with coulomb wave functions. It's very numerically stable. I'll send you it if you end up going that route to save some time.

[paulromano]

@cjosey Eventually I'd like to be able to do charged particles from LRF=7 so your toy Python script would be greatly appreciated.

[smharper]

MeV -> eV

[smharper]

Currently, these energies need to be in units of eV. I think this function should request units of MeV to be consistent with our ACE data, and convert to eV before calling _reconstruct

[paulromano]

Unfortunately the inconsistency is deeper than that. Even non-resonant cross sections are stored with energies in eV. Same goes for angle-energy distributions; everything in ENDF is eV. I think for now we should just identify that energies may be in MeV or eV depending on the source data. Longer term, I think we should consider:

• Converting data loaded from ACE files to use eV rather than MeV, and/or
• Have our basic data types know about their units so there is no ambiguity.

What are your thoughts?

[smharper]

Re short term: If the units aren't consistent then we can't use this functionality to generate HDF5 data for our transport solver. In that case, I'm not sure what the use-case is beyond plotting 0K XS. Are you just hoping to merge this chunk in now to prevent a megalithic PR later?

Re long term: I'm not a fan of the second option. I don't want to support both eV and MeV on the Fortran side, and given that constraint, I would also prefer that the PyAPI sticks to the same units as Fortran just for consistency sake. So of the two, I prefer option 1. Why do you feel that it makes more sense to convert ACE to eV rather than ENDF to MeV?

[paulromano]

There is still quite a ways to go before we could actually use this to generate HDF5 data for the transport solver, so I'm not too concerned about that for now. I'm kind of hoping to defer on the units issue because I didn't want this PR to grow too huge. I can envision a future PR just dealing with units issues. As far as the use case, in the short term I would like to use this to generate 0K elastic scattering for DBRC.

As far as MeV versus eV, I find eV to be the more "natural" unit. To me, ACE/MCNP are kind of an oddity being in MeV. I'm sure there are reasons they chose to use MeV but I'd take a guess that it has to do with their application area. And, ultimately, if we really want to go from ENDF/GND -> HDF5, it doesn't make sense to convert to MeV. Of course, if we do decide to eventually adopt eV units, users will need a fair bit of warning.

Regarding the second option (have data types know about their units), I completely agree that on the Fortran side we would not want to support that, but it wouldn't be as terrible on the Python side. I bring this up because the data abstractions in GND do have units listed explicitly. It's conceivable then that a future evaluation might contain data using non-standard units because the format allows it. You might say "no, evaluators wouldn't do something silly like that", but I've learned over time that evaluators do plenty of silly things.

[paulromano]

Tangentially related to this (I should probably just put this in a separate issue), I've long had the idea of incorporating units in our Python API for input generation. For example, having something like:

from openmc.units import inch
import openmc

s = openmc.Sphere(R=3.0*inch)

Rather than expecting a real number, the R argument in this case would expect a physical quantity with units of length. On the backend, the API would automatically convert the user-specified units to canonical units (in this case, centimeters) when generating XML files. I think it would be pretty cool, but probably not a trivial amount of work to implement.

[smharper]

"but I've learned over time that evaluators do plenty of silly things"---lol, yeah. I'm not even sure I'd be surprised if we have to convert to muonVolts in a few years.

I feel like supporting multiple units will be a real pain, even on the Python side. The tally arithmetic logic will grow a lot if we have to do unit checking and conversion. If the transport solver only outputs canonical units, then we'll have to provide functionality for users to convert units en masse that were loaded from statepoint. Third party libraries will have to handle our custom united-numbers rather than just intrinsic floats, and this will probably raise performance issues. I also don't want to contribute to the continued use of imperial units, on principle.

Regarding this particular PR, would it really be that difficult to convert everything to MeV after reading from ENDF? Unless I'm missing something, it'll only be ~20 LoC.

[paulromano]

I never said that it wouldn't be difficult, only that it wouldn't be "as bad" as the Fortran side 😄 All the complexities you've listed are among the reasons I haven't yet taken a stab at it.

For unit conversion, it will definitely take more than 20 LoC (not sure how much more, maybe ~100 LoC?). I'd need to convert units on:

• Reaction cross sections
• Reaction Q-values
• Reaction product yields
• Photon production yields / cross sections
• Prompt, total, and delayed nu (since they appear separately from other product yields)
• Angular distributions
• Energy distributions
• Fission energy release (although this might be net negative LoC?)
• Resonance parameters (and the reconstruction functions, since they're already written)

For distributions, it's often necessary to change both x- and y-values of a function since x is in units of energy and y is in units of inverse energy.

If I'm going to go through the effort of making units consistent, I'd actually rather just convert ACE data to eV rather than ENDF data to MeV, but that would have to be predicated on others' willingness to move to eV units. If you're open to that, I'll probably do it as a separate PR (which can be reviewed/merged prior to this one if you're worried about unit consistency).

[smharper]

Okay, I'm convinced. I'm okay to merge this even with the unit inconsistency. And for the record, I do agree that eV makes more sense for us.

I'll take a look soon to see if I have any other requests before merge.

[smharper]

This is the only place _ENDF_FLOAT_RE is used. Is there some reason you didn't just use it as a local variable?

[paulromano]

Yes, the regex has to be "compiled" before use, so declaring it module-level allows the compilation to happen only once which in theory would reduce the time for repeated calls to this function. After doing a little reading, apparently Python actually caches the compilation anyway, so perhaps there's not such a difference in performance.

[smharper]

Describe the skipC argument

[smharper]

I get an error with U-238 on JEFF 3.2:

ValueError                                Traceback (most recent call last)
<ipython-input-13-219fcbd89c07> in <module>()
----> 1 dat2 = openmc.data.IncidentNeutron.from_endf('/home/smharper/nuclear-data/endf/jeff32/n-92-U-238.jeff32')

/home/smharper/openmc/openmc/data/neutron.py in from_endf(cls, ev_or_filename)
    631             ev = ev_or_filename
    632         else:
--> 633             ev = Evaluation(ev_or_filename)
    634 
    635         atomic_number = ev.target['atomic_number']

/home/smharper/openmc/openmc/data/endf.py in __init__(self, filename)
    322             self.section[MF, MT] = section_data
    323 
--> 324         self._read_header()
    325 
    326     def _read_header(self):

/home/smharper/openmc/openmc/data/endf.py in _read_header(self)
    382         # File numbers, reaction designations, and number of records
    383         for i in range(NXC):
--> 384             items = get_cont_record(file_obj, skipC=True)
    385             MF, MT, NC, MOD = items[2:6]
    386             self.reaction_list.append((MF, MT, NC, MOD))

/home/smharper/openmc/openmc/data/endf.py in get_cont_record(file_obj, skipC)
    122     L2 = int(line[33:44])
    123     N1 = int(line[44:55])
--> 124     N2 = int(line[55:66])
    125     return [C1, C2, L1, L2, N1, N2]
    126 

ValueError: invalid literal for int() with base 10: '           '

[paulromano]

@smharper The JEFF 3.2 U-238 file is apparently ill-formatted (missing many MOD quantities in MF=1, MT=451). The same problem exists in the H-2 evaluation as well. Mn-55 has text encoding that is all f-ed up and Python 3 does not appreciate it. At the end of the day, I'm not sure how much we can or should do to clean up messy work by evaluators.

[smharper]

I'd like to be able MT=458 data from JEFF. Can you work in any other options for just the 458 data?

[paulromano]

Not sure I understand your last comment. Are you asking if I can put in a workaround so that all the fissionable nuclides from JEFF-3.2 work without modifying the files?

[paulromano]

Hmm, apparently there are numerous issues with JEFF 3.2 (beyond just ill-formatted files), so do hold off until I can address them.

[smharper]

Bummer. It's really a shame just how sloppy the evaluations are :(

Regarding the 458 data, I don't expect you to write an API that can identify and correct JEFF errors, but maybe modify _extract_458_data so that it doesn't require an Evaluation object. Here's a rough mock-up of what I'm thinking:

def _extract_458_data(ev):
    ...

    if isinstance(ev, openmc.data.Evaluation):
        if (1, 458) in ev.section:
            ...

    elif isinstance(ev, str):
        with open(ev, 'r') as fh:
            ...

    else:
        raise TypeError

[paulromano]

I fixed a few issues that were coming up with the JEFF 3.2 library:

1. A number of evaluations (Li7, Cl35, Eu152,153, Gd155, Ir193, Th229, Pu238,241,242, Bk250, Cf254) use an obsolete CM-to-lab energy transformation matrix in MF=4. ENDF-102 only makes one weird mention of this in the following sentence from section 4.3: "The ENDF-6 format rules do not allow transformation matrices for conversion between CM and LAB coordinate systems". So, these evaluations are breaking the rules. Nevertheless, it is easy enough to just skip the transformation matrices and continue reading as normal.
2. Pu239 from JEFF 3.2 had an energy-dependent total delayed neutron yield and an energy-dependent precursor group probability (tabulated on a different set of energies) that need to be multiple by one another to get the group yield. Since this was already dealt with for the ACE data, it was easy enough handling it on the ENDF side.
3. When no polynomials are given for fission energy release data, one is to use the Sher-Beck formulas outlined in ENDF-102. However, one of these formulas requires nu. When @smharper first implemented this, there was no support for polynomial nu. I've now implemented support for calculating prompt neutron energy release when nu is given as a polynomial. @smharper Please take a close look at how I've done it to make sure you agree.

The remaining issues that I know of are all what I consider to be broken evaluations, in three categories:

• Some are missing MOD values in MF=1, MT=451 (U235, U236, U237, and U238).
• Some have non-ASCII characters without the proper file encoding (H2, Mn55, Cd111, Gd152, Gd154)
• Some are missing resonance parameters (Cu63, Cr50, Cr54). I'm sure they are supposed to be assumed to be zero, but I don't see anything in ENDF-102 that permits just leaving a field blank and expecting a code to pick it up as a zero.

[smharper]

The Sher-Beck 458 data has some mixed-units issues. All of the Sher-Beck constants are in units of MeV, but the nu table or polynomial uses eV if it originates from ENDF.

[smharper]

Also, the JEFF-3.2 U235 and U238 might have bad formatting errors, but they do have intact 458 data that I'd like to be able to use (ENDF and JEFF differ significantly on their fission Q-values).

[paulromano]

I've implemented a fix for the missing MOD values in the U isotopes. They should work now (except for U237 which has 6 delayed neutron spectra but 8 precursor groups?! Doesn't seem right so I throw a ValueError when this happens).

[paulromano]

@smharper I believe I've handled all your requested changes. Is there anything remaining you would like me to change for this PR?

[smharper]

@paulromano, yeah I want the mixed-units for Sher-Beck 458 to be addressed somehow. If you don't want to sink a lot of time into it, maybe just raise a RuntimeError whenever someone wants to add fission energy release data to an IncidentNeutron built with the from_endf constructor.

Would also be nice if a helpful error message was output when a user tries to run to_hdf5 on data that came from_endf.

[paulromano]

Ok, I've implemented what I think is a fix for the mixed units problem. Per your request, calling IncidentNeutron.export_to_hdf5() will raise an exception if the instance was created with from_endf().

[smharper]

I think we need 'eV' in this call.

[paulromano]

That is the default for the units argument.

[smharper]

Do you mean "units that are used"?

[paulromano]

Yes, will fix that.

[smharper]

I'm running a bunch of simulations right now so I can't test your changes, but I'll try to test them later tonight and merge if it works out.

[smharper]

There's something wrong with ENDF/B-VII.1 Th-232 MF=1, MT=456 (prompt nu). The ENDF file gives a TAB1 with 20 interpolation points, but IncidentNeutron.from_endf is only reading 3.

This error isn't present for all nuclides. Pu-239 works fine.

[paulromano]

Ok, the Th-232 issue is fixed now. Thanks for bringing that up!

[smharper]

Merged! Sorry to drag you through that lengthy review.

[paulromano]

Thanks @smharper! No need to apologize; I appreciate your attention to detail!