njoy / NJOY2016

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

Charged particle processing fails in acer #87

Closed whaeck closed 6 years ago

whaeck commented 6 years ago

I created an ENDF file for incident alpha particles on Be9 based on the TENDL 2017 evaluation for this incident particle and nuclide. The original ENDF file had an upper energy of 200 MeV. For the new ENDF file, this upper energy was cut to 15 MeV (several reactions with thresholds above 15 MeV were removed from the file and outgoing particles with a multiplicity of zero below 15 MeV were removed from the file). This file was created to graft JENDL4 data on it as the JENDL4 ENDF file is incomplete and not capable of being used for transport. Proper formatting was tested using several ENDF parsers, including the one inside JANIS. This cut down ENDF file is attached to this issue.

When running the original TENDL file with the attached NJOY input file, NJOY runs to completion. It does not however using the cut down file. In this case, NJOY crashes with the following error: Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:

0 0x7feff54f027f in ???

1 0x7feff5d07100 in ???

2 0x7feff63bccf7 in acelcp

at /home/wim/Projects/NJOY2016/njoy2016/src/acefc.f90:10631

3 0x7feff63fcf31 in acelod

at /home/wim/Projects/NJOY2016/njoy2016/src/acefc.f90:6175

4 0x7feff64330ea in __acefc_MOD_acetop

at /home/wim/Projects/NJOY2016/njoy2016/src/acefc.f90:178

5 0x7feff64783a9 in __acem_MOD_acer

at /home/wim/Projects/NJOY2016/njoy2016/src/acer.f90:424

6 0x401c58 in njoy

at /home/wim/Projects/NJOY2016/njoy2016/src/main.f90:238

7 0x401e42 in main

at /home/wim/Projects/NJOY2016/njoy2016/src/main.f90:112

Floating point exception

input.txt endf.txt

whaeck commented 6 years ago

The problem is traced back to the acelcp subroutine which handles charged particle data. This routine handles the energy and angle distributions but also the recoil energy associated with charged particles. For this purpose, it will go over all of the MF6 data and retrieve the appropriate position in the ace file for the reaction under consideration.

With the slimmed down TENDL file, things seem to go wrong when the subroutine tries to retrieve the MT2 data from the ACE file using a do loop (line 10624 of acefc.f90):

         do while (mtt.ne.mt)
            ir=ir+1
            mtt=nint(xss(mtr+ir-1))
            k=nint(xss(lsig+ir-1))+sig-1
            n=nint(xss(k+1))
            iaa=nint(xss(k))
            q=xss(lqr+ir-1)
         enddo

While going over the xss array, acer tries to find MT2, which it cannot find as it is never added. the xss(mtr+ir) values are set inside an if statement explicitly excluding MT2 starting at line 5345:

         if (mth.ne.2) then
            !-- code removed between line 5346 up to and including 5361 - irrelevant for this illustration
            !--store reaction parameters
            j=j+1
            xss(mtr+ir)=mth
            xss(lqr+ir)=sigfig(scr(2)/emev,7,0)
            xss(tyr+ir)=0
            xss(lsig+ir)=next-sig+1

While trying to find MT2, we most likely go over the array size into unallocated memory because the loop can never exit until it finds a 2 in xss(mtr+ir-1). In most cases this actually appears to happen, except for the slimmed down ENDF file.

Now, acelcp actually knows it does not have to treat MT2 at this point, because on line 10647, it explicitly excludes MT2:

            !--compute the heating from this recoil nuclide
            if (izap.gt.2004.and.mt.gt.2) then

but it only checks this after the do loop.

A fix for this issue is to change where we check for MT2 by changing line 10620 to:

      if (mfh.eq.6.and.mth.gt.2) then

and line 10647 to:

            if (izap.gt.2004) then

The tosend call on line 10786 should be moved to after the endif statement to ensure that we skip over the MT2 section properly.

whaeck commented 6 years ago

Testing has shown that for the original TENDL2017 file, the produced ACE file does not change when applying these changes. The NJOY run using the cut down ENDF file now runs to completion without crashing.

whaeck commented 6 years ago

Changes added to pull request #88.