search

njoy / NJOY2016

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

Consistency checks introduces spikes in the edges of the secondary energy distributions #130

Closed jchsublet closed 4 years ago

jchsublet commented 5 years ago

When processing the Fe56 from TENDL 2019, we found a potential problem in the way ACER corrects some data issues in the file while running the consistency checks.

The first issue is related to the epmax value determined by NJOY. Whenever the upper energy of the secondary energy distribution is larger than this value, consis will reset the energy values that are larger to be below this value and "renormalise" the pdf function associated to it. The normalisation does not appear to be correct as this leads to the introduction of artificial spikes in the energy distributions.

Whenever we change to secondary energy distributions to have an upper energy limit equal to the NJOY calculated epmax value, the distributions look as they should, as intended by the evaluator.

A second issue is also related to consis but cannot be corrected in the actual evaluation. ENDF rules (apparently) require MT5 crossections and other related data to begin at 1e-5 eV because the reaction Q value is zero and it is therefore not a threshold reaction. consis also changes the epmax (MF6 MT5 neutron production) given in the COM frame as it should to a LAB frame, but because the data is below the MT5 cutoff, the evaluator only put a dummy distribution here. In this case, it consists of a distribution at an incident energy of 1e-5 eV with secondary energies between 0 and 1e-5 eV. consis wants to change this last value in 9.648869e-6 eV. Suffice to say, this makes little sense.

This is the original TENDL 2019 file: n-Fe056.tendl-19.txt This is the corrected file using the epmax value from NJOY: n-Fe056.tendl-19-c5.txt

This are the input files we used to illustrate the problem: inFe056.txt inFe056c5.txt

These are the pdf graphs produced by acer for the original and corrected evaluations (after consis made its corrections), see page 57: Fe056chk.pdf Fe056c5chk.pdf

jlconlin commented 5 years ago

Thanks @jchsublet I've assigned this to @whaeck to fix it since he thinks he knows how to do it already.

jchsublet commented 4 years ago

I'll be most grateful for a move on the first issue: the "spike" creation, wrong re-normalisation as I have many new files coming up. TEFAL has been corrected to better follow NJOY epmax rules but others file exits and the re-normalisation sequences contains a bug

whaeck commented 4 years ago

@jchsublet Sorry about the delay, I have been getting some backlog out of the way on NJOY2016 before I got to this one.

I have been looking at the code to see what is going on. It essentially happens in these lines in the consis subroutine (starting at line 14717 of acefc.f90) :

                  if (n2big.gt.0) then
                     write(nsyso,'(''   consis:'',&
                       &'' shifting eprimes greater than epmax'',&
                       &'' and renorming the distribution'')')
                     do j=nn-n2big+1,nn
                        ishift=j-nn-1
                        ep=xss(j+loci)
                        xss(j+loci)=sigfig(epmax,7,ishift)
                        if (intt.eq.1) then
                           p=(xss(j+2*nn+loci)-xss(j-1+2*nn+loci))&
                            /(xss(j+loci)-xss(j-1+loci))
                           xss(j-1+nn+loci)=p
                           xss(j+nn+loci)=p
                        else
                           p=2*(xss(j+2*nn+loci)-xss(j-1+2*nn+loci))&
                            /(xss(j+loci)-xss(j-1+loci))&
                             -xss(j-1+nn+loci)
                           xss(j+nn+loci)=p
                        endif
                     enddo
                  endif
               enddo

When the consis routine determines the maximum outgoing energy is too high compared to epmax, it sets the last value of the outgoing energy array to this value, and all values before it to monotonically increasing values that are lower than epmax. It then adjusts the pdf values to maintain the associated cdf. In total the last n2big values of the outgoing energy values and pdf values are impacted.

As an example, for the iron example you provided, these are the values from the first acer run for two incident energy values (1e-11 MeV and 30 MeV):

E_i = 1E-11    n2big = 22    epmax = 9.6488690000009646E-012    intt = 1
           1   0.0000000000000000        110.02290000000001        0.0000000000000000
           2   3.8722100000000000E-005   172.61089999999999        4.2603192700000000E-003
           3   4.6048399999999998E-005   188.23290000000000        5.5249185499999997E-003
           4   5.4760899999999999E-005   205.26990000000001        7.1648976699999999E-003
           5   6.5121800000000004E-005   223.84690000000001        9.2916784399999994E-003
           6   7.7442999999999997E-005   244.10589999999999        1.2049740600000001E-002
           7   9.2095400000000006E-005   266.19990000000001        1.5626477400000002E-002
           8   1.0951999999999999E-004   290.29280000000000        2.0264903399999999E-002
           9   1.3024200000000000E-004   316.56479999999999        2.6280351600000001E-002
          10   1.5488399999999999E-004   345.21580000000000        3.4081142100000003E-002
          11   1.8418800000000001E-004   376.46080000000001        4.4197346200000001E-002
          12   2.1903700000000000E-004   410.53280000000001        5.7316628400000000E-002
          13   2.6047899999999998E-004   447.68880000000001        7.4329927700000006E-002
          14   3.0976199999999999E-004   488.20670000000001        9.6393372599999999E-002
          15   3.6837000000000003E-004   532.39170000000001       0.12500619299999999
          16   4.3806599999999998E-004   580.57669999999996       0.16211176499999999
          17   5.2094900000000004E-004   633.12170000000003       0.21023170199999999
          18   6.1951400000000005E-004   690.42359999999996       0.27263533800000000
          19   7.3672699999999995E-004   752.90959999999995       0.35356196200000001
          20   8.7611700000000004E-004   821.05359999999996       0.45851003000000001
          21   1.0418800000000000E-003   895.36249999999995       0.59461032899999999
          22   1.2390100000000000E-003   976.39649999999995       0.77111314099999995
          23   1.4734300000000000E-003   0.0000000000000000        1.0000000000000000

E_i = 30.0    n2big = 1    epmax = 28.946600000002896    intt = 1
         108   28.731880000000000        2.7093920000000001E-003  0.99866930899999995
         109   29.223020000000002        0.0000000000000000        1.0000000000000000

After running consis, these values have become the following:

E_i = 1E-11    n2big = 22    epmax = 9.6488690000009646E-012    intt = 1
           1   0.0000000000000000        441536617.79480743        0.0000000000000000
           2   9.6488470000009648E-012   1264599279269990.5        4.2603192700000000E-003
           3   9.6488480000009654E-012   1639979119053297.2        5.5249185499999997E-003
           4   9.6488490000009660E-012   2126780772208282.5        7.1648976699999999E-003
           5   9.6488500000009649E-012   2758062158407868.0        9.2916784399999994E-003
           6   9.6488510000009655E-012   3576736803713805.0        1.2049740600000001E-002
           7   9.6488520000009645E-012   4638425997322396.0        1.5626477400000002E-002
           8   9.6488530000009651E-012   6015448196527493.0        2.0264903399999999E-002
           9   9.6488540000009656E-012   7800790508099733.0        2.6280351600000001E-002
          10   9.6488550000009646E-012   10116204094160266.        3.4081142100000003E-002
          11   9.6488560000009652E-012   13119282192426694.        4.4197346200000001E-002
          12   9.6488570000009657E-012   17013299317665288.        5.7316628400000000E-002
          13   9.6488580000009647E-012   22063444887263532.        7.4329927700000006E-002
          14   9.6488590000009653E-012   28612820383482804.        9.6393372599999999E-002
          15   9.6488600000009659E-012   37105572038527536.       0.12500619299999999
          16   9.6488610000009648E-012   48119936972222032.       0.16211176499999999
          17   9.6488620000009654E-012   62403635963976560.       0.21023170199999999
          18   9.6488630000009660E-012   80926624084027888.       0.27263533800000000
          19   9.6488640000009649E-012   1.0494806793941714E+017  0.35356196200000001
          20   9.6488650000009655E-012   1.3610029914131590E+017  0.45851003000000001
          21   9.6488660000009645E-012   1.7650281189811101E+017  0.59461032899999999
          22   9.6488670000009651E-012   2.2888685886787162E+017  0.77111314099999995
          23   9.6488680000009656E-012   2.2888685886787162E+017   1.0000000000000000

E_i = 30.0    n2big = 1    epmax = 28.946600000002896    intt = 1
         108   28.731880000000000        6.1976200455596125E-003  0.99866930899999995
         109   28.946590000002896        6.1976200455596125E-003   1.0000000000000000

The reason why the pdf values become so big is due to the use of sigfig to modify the last significant digits to obtain a monotonically increasing outgoing energy grid. This is evident in the example above for the first incident energy value.

After analysing the algorithm, I still think there may be issues with the way this is done. Instead of using sigfig to produce outgoing energy values that are one significant digit apart, we could use equal distant energy values or renormalise the outgoing energy grid to maintain the same relative distance in each of the affected bins.

Or, we could keep the actual energy grid and change the first energy value that is larger than epmax to be equal to epmax and change the pdf value so that the associated cdf value becomes 1. All other pdf and cdf values of the next energy values is then to be set to 0.0 and 1.0 respectively.

I'll try out a few things.

whaeck commented 4 years ago

I tried the scheme indicated at the end of my previous comment, and acer now produces plots that are more appropriate. However, calling acer multiple times on the corrected ace file degrades the file even further. In addition, I do not know if MCNP will actually accept a pdf/cdf like this where multiple points at the end have a cdf value of 1.0.

In addition, I have an issue with the way the probabilities are calculated for pdf/cdf with linear interpolation. For every energy above epmax, consis sets the following probabilities:

                           p=(xss(j+2*nn+loci)-xss(j-1+2*nn+loci))&
                            /(xss(j+loci)-xss(j-1+loci))
                           xss(j-1+nn+loci)=p
                           xss(j+nn+loci)=p

It uses the cdf values and then sets the previous energ's probability and the current one to the same value. This becomes problematic when multiple energies are affected. The resulting pdf and cdf will no longer be in sync. This can be seen in teh example above for the 1e-11 MeV incident energy point: all but the last two pdf values are different because of this.

jchsublet commented 4 years ago

Thanks Wim to have make a move on this one, it is appreciated, however let us keep in mind the physics behind all this, not only the numerics. It does make sense to output outgoing particles spectra below the energy of the incident particle (not always though) but doing so below 1e-11 MeV (1e-5 eV) does not..for NJOY or anybody else make much sense

Regarding when epmax in the file is well above the minimum energy but still too high when compared with the value calculated by NJOY one would like to take a simple approach: remove, skim all data (energy point and KM parameters) when ep-file is > ep-mjoy (sigfig 4 is more then enough there) then terminate with the njoy-epmax pdf=0.0 cdf=1.0, the last KM values in the file.

You certainly noticed the with or without the missing pdf, the cdf will be one, or even has been for a while (does not make sense either). I am always amazed by the arrogance of the numerical precision given or forced on the values of the energy grid and Kalbach-Mann parameters. We are after all taking about Kalbach-Mann systematics representation.

whaeck commented 4 years ago

I take the opportunity of getting these issues to delve deeper into the source code, this understanding will pay off later when I have to redevelop some of these things in NJOY21.

In my opinion, these things should be done when creating the actual ACE file (modifying the XSS array is easy at that point).

I will try the following next (assuming that the original secondary energy distribution has N points and that the first secondary energy value that is larger than epmax is point Nnew):

jchsublet commented 4 years ago

Indeed, this seems a good path to follow, having said that in light of what you've said above and the importance-relevance of the KM parameters in the tails: f0 and r from which a is derived I would argue that the simple following steps are acceptable, would respect the evaluator original wish:

reset epmax to epmax-njoy with f0, r at 0.000000+0 remove all ep > epmax-njoy, e.g. x energy points reset epmax-1 f0 to the sum of the f0 (of the x points removed) with r equal to 1.000000

In doing so one shift the epmax of the energy grid to the NJOY calculated value but also keep (sum-up) the pdf (f0) of the removed energy points (no renormalisation needed). Lets us not forgot that the parameters are valid between ep-1 and ep.

whaeck commented 4 years ago

I have implemented the changes proposed above into acefc (see the branch fix/acer-consis). The current implementation has some temporary print statements in it to verify that the table is modified correctly. No locator changes have to be made, only the number of outgoing energy values are modified. In the output file for the second acer run which performs the consis corrections, the data is being represented correctly: e.g. the incident energy point at 1e-11 MeV no longer has 23 outgoing energy points but only two, as it should be after the correction. The entire ace data print out in the output file appears to be correct. The graphs also appear to be correct, without the large spikes.

However, this is where things go wrong again. While the xss array is output correctly in the acer output file (along with the plots being correct), it is not correct in the actual corrected ace file (tape34 in your example). This is due to one rather annoying fact: while the interpretation of the xss array in acer for the output file takes into account the locators in the file, the routine writing the xss array to the actual ace file does NOT. aceout and the underlying change routine (which prints either an int or double depending on context) assumes the xss array is CONTIGUOUS and ignores locators.

As a side note: I was wrong about the interpolation. It is being done correctly (I believed 1 was linear interpolation but it is actual the histogram representation).

jchsublet commented 4 years ago

Am I to understand that the 'correction' works only for the printed output and the graphs but not the actual output tape that serve MCNP? does this occurs on every occasions of consis interferences of just that one? What can I do to help?

whaeck commented 4 years ago

When consis corrects the data, we do this on the xss array read from the first ace file. Once the changes are made, acer goes over the xss array to print the data to the output file. Once that is done, the aceout subroutine is used to create the actual new ace file.

The problem is that aceout does not look at the locators in the xss array. Taking the example of the secondary energy distribution for the 1e-11 MeV incident energy: aceout sees that the number of energies is 2 so it writes 5*2 values (Eprime, pdf, cdf, r and a for every outgoing energy). At that point, aceout thinks it is done with the secondary distribution for this incident energy and it then moves on to the next one. However, because we zeroed out 21 energies, we still have to write 105 values in the xss array before we actually get to the next energy distribution. The reason why the ace data printing in the output file works, is because that actually looks at the locators, contrary to aceout. As a result, aceout transforms data into integers or doubles, depending on what context aceout thinks it is in. If we output the entire xss array using doubles only, we shouldn't have a problem.

I already have a partial solution by modifying aceout for law=44 using the secondary energy locator values (except at the end of the last energy distribution because I do not have a locator, it is supposed to be the next block). Once I get this to work for our law 44 specific issue, I'll expand this aceout correction over most blocks.

whaeck commented 4 years ago

@jchsublet Ain't this pretty? :-) tape42.pdf

consis now corrects the xss array by resetting a single energy point, shifting the pdf, cdf, r and a values and "erasing" the remainder of the array for the secondary energy distribution where the change was made. The locators are left intact and now NJOY prints the ACE file using the actual locators so that these "gaps" in the xss array are actually taken into account.

Now that this is done, I have to go over all the other consistency checks and laws to add locator advancing (just in case).

whaeck commented 4 years ago

I forgot to say: the changes are available in the fix/acer-consis branch (there still are print statements in the code but those are temporary, I use them during development to check if the gaps appear in the proper places).

jlconlin commented 4 years ago

@whaeck That plot does look good. Can you please show the same plot from before the change was made? I'm not sure I fully comprehend what the changes are.

whaeck commented 4 years ago

@jlconlin JC added a pdf of his own, see the top of the page: Fe056chk.pdf

jlconlin commented 4 years ago

Ah yes, thanks for reminding me. The new plots do look much better!

jchsublet commented 4 years ago

Absolutely beautiful, professionalism at its best. Now please teach me how to access/get the corrected branch

whaeck commented 4 years ago

To get this branch, you have to do the following:

cd < njoy2016 source location >
git pull
git checkout fix/acer-consis

and you'll be ready to build it. njoy ;-)

I have also started on symplifying the routine that writes the xss array to a file by integrating repetitive code into a few subroutines (write_real_list, write_integer_list, etc.). That'll make things easier to read later on.

I'll also try to collect a few files with consistency problems and add them as test files once all this work is done.

Keep in mind that we still need to test these new ace files with MCNP to see if MCNP does indeed use only the locators (I have been told that it SHOULD be the case).

whaeck commented 4 years ago

@jchsublet I now applied the fix for law=44 to law=4 (it's practically the same, except that there is no array for r and a as in law=44). I also changed the messages from consis a bit: it will now tell you how many values we "remove", what MT and LAW the change was applied to.

Do you know of a case where consis does this for a law=4? I would like to test the changes on that case as well. I added your example of Fe56 as test 55 in the fix/acer-consis branch. Once I have a case for LAW=4, I'll add that as test 56.

Once all that is done, I'll finalise this and perform a pull request to put this up as NJOY2016.54

jchsublet commented 4 years ago

@whaeck excuse the delay in answering, out of my control. No I do not know or spotted such behaviour for a law4, this one was pretty difficult to spot already

whaeck commented 4 years ago

@jchsublet I went ahead and made the pull request to incorporate these changes into master. I have not found any LAW=4 case with the same problem so I have not added a test for that particular case.

Would you be willing to verify my changes independently?

jchsublet commented 4 years ago

@whaeck I just downloaded and tested 55 but sorry I do not even see the changed you've made in 53 regarding acer consis eprime they did not make it to this acefc.f90. Or I did no tested what you wanted me to test so please specify what you want me to do

whaeck commented 4 years ago

The NJOY version with those acefc fixes was in a separate branch (granted, the version number on that one was 53 because it was branched off from NJOY2016.53) and we're pulling it into the master branch now. The official version of NJOY containing these changes will most likely be 2016.56.

See the latest pull request: #147

jchsublet commented 4 years ago

So I already tested it as my .53 had this branched in, learning the hard way github handling instead of upd, poor me

jlconlin commented 4 years ago

This was resolved in #147