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 energy bin sampling in source #941

Merged
merged 7 commits into from Dec 29, 2017
Merged

Conversation

brbass
Copy link
Contributor

@brbass brbass commented Dec 9, 2017

The algorithm to find the appropriate energy group for a sampled particle in source.F90 assumes the energy bounds are ordered from lowest to highest. The current implementation reverses the energy bins to be from highest to lowest in input_xml.F90.

This change reverts the energy groups to their original order, from lowest to highest. This has a side effect of aligning the energy bins with the tally bins, which previously were reversed.

After the change, test_mgxs_library_ce_to_mg.py fails. It looks like the the difference is probably due to a change in order of some random operation on the flipped array of energy bins. As I increase the number of particles, the old and new results converge to the same eigenvalue and same cross sections. I went ahead and updated the results of the test for now.

@paulromano
Copy link
Contributor

Unreversing energy_bins would also have the effect of unreversing energy_bin_avg which is used as a proxy for the particle's energy. @nelsonag Wouldn't this result in tallies with energy filters being wrong?

@nelsonag
Copy link
Contributor

@paulromano Yes. And that only matters if the tally's energy/energout filters bins dont match up with the group structure; if they do match then the particle's group is used to set the filter bin.

@brbass The energy bin indexing is set up as follows:

  • group structure is input similar to tally energy bins: from lowest to highest
  • code reverses this group structure so that a group index is ordered as is usual (group 1 == fast, group G == thermal)
  • The code then transports in normal group ordering
  • The tallying system uses energy/energyout filters ordered from lowest energy to highest energy, so at tally time consistency must be ensured.

After digging deeper and reminding myself of the above, I think instead this PR might want to focus on fixing source.F90 while minimallly perturbing the rest of the source that implements the above energies. That code looks to be simply incorrect (even if there wasnt confusion about the ordering). In addition to the errors you've pointed out, the code should probably be checking for sampled energies outside the bounds of the problem and erroring if found.

@paulromano
Copy link
Contributor

I agree -- easiest thing to do right now is to just fix source.F90 and leave everything else as is.

@paulromano paulromano added this to the v0.10 milestone Dec 11, 2017
@brbass brbass changed the title Remove energy bin reversal in multigroup mode Fix energy bin sampling in source Dec 11, 2017
@brbass
Copy link
Contributor Author

brbass commented Dec 11, 2017

I reversed the energy groups back to their original format. Since the binary search used to find the energy group in source.F90 requires that the data be ordered from smallest to largest, I added the rev_energy_bins to mgxs_header.F90 and correctly initialized it in input_xml.F90. This requires that the energy group be reversed after the search. Finally, I changed the tally_filter_energy.F90 to use the rev_energy_bins array instead of reversing energy_bins.

@brbass
Copy link
Contributor Author

brbass commented Dec 11, 2017

As a side note, it might be a good idea to add a multigroup, fixed-source problem with an energy-dependent source to the tests. Since the incorrect code didn't raise any flags, I'd guess there aren't any tests that would fail if the source particle was always sampled in the group with index one.

@nelsonag
Copy link
Contributor

nelsonag commented Dec 12, 2017 via email

@brbass
Copy link
Contributor Author

brbass commented Dec 14, 2017

@nelsonag I could add the fixed-source MG test sometime in the next few weeks. For the purposes of testing the code we've changed here, a three-group infinite medium problem with an energy-dependent source and an analytic solution would be a good sanity check. I'm not familiar enough with the existing tests to say if others are needed.

@nelsonag
Copy link
Contributor

@brbass The changes so far look good. We might as well wrap in the checking the MG-sampled source to make sure the energy sampled is within the range of the multi-group structure. To get the high-end energy we simply need to set the energy_max_neutron and energy_min_neutron in input_xml.F90 (near where you added back in rev_energy_bins = energy_bins).

The maximum energy will then be checked by the SourceDistribution % sample method.
The minimum energy is a bit more complicated the continuous-energy mode doesn't check the minimum energy of a source. So we either can add it in to SourceDistribution % sample so both CE and MG have a minimum energy, or we can have the code you modified in souce.F90 for this PR check against it so it only affects MG mode.

I personally vote for the latter; i suspect its more likely for this to be an issue in MG mode than CE and therefore dont think the additional if-then is necessary for every source sample in both modes. Thoughts?

@brbass
Copy link
Contributor Author

brbass commented Dec 14, 2017

In terms of code simplicity, putting the minimum energy check inside sample() in SourceDistribution would be preferable, since in the sample_external() function of source.F90 we would then just do the binary search directly without adding an additional check. This would also mean we wouldn't need to add the nuclide header to source.F90. My issue with checking in the source.F90 code for the minimum energy would be that we would have to either (1) decide to throw a fatal exception if the energy was below the minimum energy boundary or (2) put in another "do" loop for resampling the external source if the energy was below the minimum. It seems like the first option would be overly restrictive for continuous distributions, since the user may wish to specify a distribution and then only accept energies inside the MG energy bounds. The second option would be pretty clunky, since there's already a rejection mechanism built into the SourceDistribution code. Would adding a check for minimum energy in the SourceDistribution break anything in continuous energy land?

As an aside, is there some reason why the discrete source is checked every time for incompatible energies in source_header.F90? Unless the discrete sources are changing after reading the input, it would be a lot more efficient to check the sources for compatible energies just after the sources and cross sections are read in.

@paulromano
Copy link
Contributor

All very good points @brbass. I agree with you that it is probably cleaner to have that check in SourceDistribution % sample. In fact, probably what we should do is set the energy_min_neutron and energy_max_neutron global variables (that are currently only used for CE data) when the MG cross section data is read in. That way, the check in sample() only needs to check if energy_min_neutron < site % E .and. site % E < energy_max_neutron and it will work for both CE and MG data. I don't foresee any problems rejecting energies below the the minimum.

Regarding the discrete distribution check, the reason it's probably there is because at the time the source distributions are read in (as part of settings.xml), the cross sections haven't been loaded yet, so the check has to be deferred until afterwards, but having a check on the source at the time the cross sections are loaded might appear a little awkward, so we probably found it a little more logical to have it in source_header.F90

@paulromano
Copy link
Contributor

I'm going to go ahead and merge this as it looks fine. @brbass If you end up writing some tests for this, you can submit a separate PR.

@paulromano paulromano merged commit 9374695 into openmc-dev:develop Dec 29, 2017
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants