njoy / NJOY2016

Nuclear data processing with legacy NJOY
https://www.njoy21.io/NJOY2016
Other
96 stars 85 forks source link

THERMR failing when iform=1 #27

Closed mattooca closed 6 years ago

mattooca commented 7 years ago

Trying to process some thermal neutron scattering law data from the latest ENDF-VIII ln -s endfSVN/neutrons/neutrons/n-040_Zr_090.endf tape20 ln -s endfSVN/thermal_scatt/tsl-ZrinZrH.endf tape21

Then run njoy2016 with the following input file:

-- Process TSL data for Zr in ZrH
 'reconr'
 20 22 /
 'pendf tape for MAT=4025 evaluation'/
 4025 1 0 /
 .001/
 'HinZrH from ENDF'/
 0/
 'broadr'
 20 22 23
 4025 1/
 .001/
 296/
 0/
 'thermr'
 21 23 24 /
 58 4025 8 1 2 1 1 1 229 2 /
 296 /
 0.005 5 /
 'stop'

Results:

njoy 2016.0   21Dec16                                       05/31/17 17:20:56
 *****************************************************************************

 reconr...                                                                0.0s

 broadr...                                                                0.9s

 thermr...                                                                1.9s

 ***warning***maximum value of beta limits the allowed energy transfer
 the sct approx. will be used for transfers larger than  1.000 ev.
At line 56 of file /g/g16/mattoon1/apps/NJOY2016/src/endf.f90 (unit = 12, file = '/var/tmp/mattoon1/gfortrantmpIyuTTK')
Fortran runtime error: Bad value during floating point read

The same input file works with njoy2012.50 (although the resulting tape24 contains lots of "-Inf+***" for the MT=229 incoherent inelastic cross section)

Changing iform to 0 works better, gives back reasonable results for MT=229 cross section. However, iform=1 is very useful for comparing to ENDL results (with iform=1 the distribution is given as a pointwise function of mu rather than as Legendre expansions).

jlconlin commented 7 years ago

I have repeated the calculation and can confirm your results. We'll take a look at this.

whaeck commented 6 years ago

I was looking through the backlog of issues and I saw this one so I tried it. It still fails on NJOY2016.30.

I do not know if this is still an issue for you Caleb, but this is due to you using iform=1. If you use iform=0, the problem goes away (all LANL processing of Sab files uses iform=0). I think that there is a division by zero when you use iform=1. I will have to look in more detail to see what exactly is going on.

mattooca commented 6 years ago

Thanks for taking a look! iform=1 is handy for us since it is similar to how outgoing neutron distributions for thermal scattering law data are stored in the ENDL library: P(mu | incident energy) * P(outgoing energy | incident energy, mu). Last year I did some comparisons between NJOY2012 and FUDGE for tsl data, and will probably do more of those comparisons eventually. Not high priority for me right now, though. Cheers, Caleb

On Thu, 22 Mar 2018 at 12:50 Wim Haeck notifications@github.com wrote:

I was looking through the backlog of issues and I saw this one so I tried it. It still fails on NJOY2016.30.

I do not know if this is still an issue for you Caleb, but this is due to you using iform=1. If you use iform=0, the problem goes away (all LANL processing of Sab files uses iform=0). I think that there is a division by zero when you use iform=1. I will have to look in more detail to see what exactly is going on.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/njoy/NJOY2016/issues/27#issuecomment-375435684, or mute the thread https://github.com/notifications/unsubscribe-auth/AFYQL8JFy64DjSaAzC9xWj_5_j1afnFOks5thACZgaJpZM4NsW7G .

whaeck commented 6 years ago

I found the problem. It is due to the use of the incident energy array not being initialised properly when iform=1. That goes to show that nobody actually uses this option, as it was also in NJOY2012 so it has been in there for a very long time. Non regression tests, people!

Anyway, just add the following line at line 2234 of therm.f90 (4 lines under label 520): esi(ie)=enow

This will solve your issue. However, as apparantly nobody has ever tested this, I need to be sure that this part of the code actually does its job correctly. Soooo ..., seeing as you brought it up, would you be willing to look at the numbers coming out to see if it is OK?

I will use your input to create a non-regression test for iform=0 and 1 so that we can detect this in the future.

njoy, Wim

whaeck commented 6 years ago

The fix is now in pull request #74, I added a non regression test as well for these format options.

If you have the time Caleb, checking the results on your cases would be appreciated to see if there are other problems.

This issue will now be closed.