Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix/acer consis #147

Merged
merged 19 commits into from
Mar 31, 2020
Merged

Fix/acer consis #147

merged 19 commits into from
Mar 31, 2020

Conversation

whaeck
Copy link
Member

@whaeck whaeck commented Jan 22, 2020

This addresses issue #130

This relates to changes made by the consistency checks when a checker acer run is requested. NJOY only rarely modifies data. It does so for secondary particle distributions that use LAW=4 (isotropic angular distribution and continuous tabulated energy distributions) or LAW=44 (Kalbach-Mann).

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. This often leads to the introduction of artificial spikes in the plots produced by acer.

This behaviour has now been modified. Whenever this problem is detected, the consis routine will now reset the first energy larger to epmax to a value slightly below it and adjust the pdf value so that the cdf value for this energy point is 1.0. All other energy points after this energy point along with their associated data is set to zero. The locators in the xss array are left unchanged.

Due to this new behaviour, and because the routine that writes the xss array does not check locators, the change routine was modified as well to verify consistency of the locators with the actual position in the xss array. Whenever the locator is inconsistent, it will either just advance to the locator (thus taking into account the gaps created by the consis routine) or issue an error when the locator points to a position lower than the current one.

The change routine was also simplified to improve its readability by adding subroutines to write a list of integers, a list of real numbers, a single integer and a single real number.

A large number of compiler warnings were also fixed in this pull request (mostly unused or uninitialised variables).

Test case 55 was added to properly test this for a LAW=44. Up to now, no case for LAW=4 has been found.

@coveralls
Copy link

Coverage Status

Coverage remained the same at ?% when pulling 7fde55b on fix/acer-consis into d0cab91 on master.

@whaeck
Copy link
Member Author

whaeck commented Mar 25, 2020

@jlconlin If you don't have enough pull requests, here's another :-)

Copy link
Member

@jlconlin jlconlin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I love it when we reduce the amount of code in a project. The few functions you have introduced here help.

@@ -4794,7 +4850,7 @@ subroutine acelod(nin,nedis,suff,matd,tempd,newfor,mcnpx,ismooth)
use util ! repoz,dater,error,skiprz,sigfig
use endf ! provides endf routines and variables
! externals
integer::nin,nedis,matd,newfor,mcnpx,ismooth
integer::nin,matd,newfor,mcnpx,ismooth
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we need to include suff in this list? (I'm not a Fortran guy)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

suff is a real, declared on the line below.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So you are right <feelingSheepish> I don't know why I missed that earlier.


!--law 4 and law 44
! this piece is NOT compatible with LAW=44 (need two additional
! arrays for r and a), potential issue?
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know if we have test that exercises LAW=44. Why is this a problem?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to the ACE specifications, there is supposed to be a difference.

Copy link
Contributor

@nathangibson14 nathangibson14 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe I'm just slow, but I'm struggling to see how the new epmax algorithm is different than the old one. The old one looks like it set a bunch of ep values to epmax and adjusted PDF/CDF values so that integrating to epmax instead of max(ep) worked right. The new algorithm does this (perhaps cleaner?) by only having one value at epmax. But how does this avoid whatever undesired behavior was observed before?

My apologies if I'm missing something that should be clear...

xss(next+nn+n)=pl ! probability
xss(next+2*nn+n)=yn ! cumulative probability
xx=elow ! 1e-5
do while (xx.lt.test1) ! as long as xx is smaller than 0.99999
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm fine with these comments, particularly in legacy code. But it's a little cringy to hardcode numbers (even in comments) that could be changed elsewhere later. :)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

True, these are comments I added trying to understand the code.

enddo

return
end subroutine write_real_list
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like these simple functions wrapping functions with less clear meaning. 👍

do i=1,ntr
call typen(l,nout,1)
l=l+1
call advance_to_locator(nout,l,sig+nint(xss(rlocator))-1) ! sig=jxs(7)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At first glance, without knowing exactly what the code is doing here, it looks like the logic has changed here. Was this a bug previously?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With the change in consis, it is possible that the xss array has a "gap" in them when we expunge them. To solve this, I decided to verify the locators in the xss arrays since the original code here assumed that the entire xss array was continuous (i.e. no gaps). By not doing this check, the start of a gap in the xss array would have been interpreted as the beginning of a new set of data. Since only doing this where I introduced gaps was a bad idea, I decided to perform thes checks for the entire xss array.

do j=3*nnew+1,3*nn
xss(loci+j)=0.0
enddo
endif
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Didn't work through the math here, but it certainly reads like the description of the changed algorithm you provided in the PR comment. And it looks about as clean as could be expected, given you had to fit into what was in this file already.

l=l+1
enddo
call write_integer(nout,l) ! NYP
call write_integer_list(nout,l,nyh) ! MTY (NYP values)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It was hard not to fall asleep looking at these repetitive changes. Props to you for working through all of them (potentially without falling asleep yourself)!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You have no idea how many times I verified myself to be sure not to screw this up.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm glad, because I have no faith in my review of these, haha.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I actually found an error in one of the tests because of this (I had lots of print outs to the screen). One test runs multiple ace files and there appeared to be junk in it because the array was not erased at the end of each run. The locator check was going ballistic over it.

enddo
if (nnew.ne.nn) then ! move values if required
do i=1,4
do j=i*nnew+1,5*nn
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doesn't this upper limit only need to go to i*nn? Aren't you copying junk from above the limit of this dataset here, then zeroing it out later?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. But it feels safer to it this way. The xss array we're changing essentially looks like this:

x1 x2 x3 y1 y2 y3 z1 z2 z3 f1 f2 f3 g1 g2 g3

And we want to make it look like this:

x1 x2y1 y2 z1 z2 f1 f2 g1 g2 0 0 0 0 0

So, the first time I go through this, the xss array will be:

x1 x2 y1 y2 y3 z1 z2 z3 f1 f2 f3 g1 g2 g3 g3

The second time:

x1 x2 y1 y2 z1 z2 z3 f1 f2 f3 g1 g2 g3 g3 g3

and so on.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not worried about the number of flops, although my suggestion would save a few. But what I am concerned about is what happens when j=5*nn. Your first pass would actually result in:
x1 x2 y1 y2 y3 z1 z2 z3 f1 f2 f3 g1 g2 g3 JUNK
Subsequent passes would increase the number of JUNK values.

But fortunately, xss is bigger than this dataset, so chances of segfaulting are minimal.

@whaeck
Copy link
Member Author

whaeck commented Mar 25, 2020

@nathangibson14 The reason why we go from the old treatment in consis to the new treatment was to clean up some of the graphs we were producing. If you had 2 eprime points that were above the calculated epmax (let's say this is 1e-10 MeV), njoy set the last two eprimes to 9.999999999998e-11 and 9.999999999999e-11 (give or take a few 9s). The delta E for these bins is so small that the probability explodes. As a result, the graph only shows spikes. See #130 for an example.

@nathangibson14
Copy link
Contributor

Aha, I see what I was missing. The sigfig routine is changing the value from epmax. And to be clear, this shouldn't meaningfully change the values produced in downstream calculations sampling, it just makes prettier plots. I'm happy with that.

@whaeck
Copy link
Member Author

whaeck commented Mar 26, 2020

One could say it is a cosmetic change with a little more changes behind it ;-)

It did allow me to write up the remainder of the continuous energy ACE specification, maybe I should do that pull request next.

@jlconlin
Copy link
Member

Do we know if @jchsublet has had a chance to try this and whether or not it fixes issue #130

@jchsublet
Copy link

@whaeck @jlconlin yes you should know that jchsublet tested parts of this branch as he unadvertantly use it as is njoy2016.53 code when released see my last comment on #130. Having said that I'll be happy to test any other release when it occurs or even a test one

p.s. your timeline is too long between bug report and correction to be efficient

@jlconlin jlconlin merged commit 21e8d64 into master Mar 31, 2020
@jlconlin
Copy link
Member

Thanks @jchsublet We need to get better at speeding up.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Consistency checks introduces spikes in the edges of the secondary energy distributions
5 participants