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

Survival Normalization #2673

Closed
wants to merge 43 commits into from

Conversation

yrrepy
Copy link
Contributor

@yrrepy yrrepy commented Aug 29, 2023

Please close #2650 as I will be attempting to shepherd this through.
Apologies for the long read ahead.

Description

This PR creates the functionality of Survival Normalization, perhaps better described as Survival Biasing Source Normalization. It is meant to complement and enhance OpenMC Survival Biasing in certain situations. This functionality is used to make the weight parameters (weight and weight_avg (known in the C++ code as weight_cutoff and weight_survive) relative to the weight of the source.

Situations where this should be of benefit is when numerous important particles are present on the source at weights below that of the weight_cutoff, normalization prevents them from being rapidly rouletted. This should notably be the case for biased general sources (which OpenMC currently lacks natively).
It is currently only applicable to phase-space sources (.mcpl and .h5). Should/could later be expanded and enabled for general (biased) sources.

Background

In comparisons between MCNP and OpenMC with survival biasing and with phase-spaces I kept seeing significant differences in uncertainties. This was partly complicated by the tracklength heating being disabled in np. Eventually I tracked it down to how the two are respectively handling survival biasing.
Unfortunately my analysis erred and this feature will not improve my work, I do believe however that it is still a useful feature and will aid some. (In my case FOMs decrease by ~2%).

For phase-space sources, the weight parameters are made relative to the starting weight of each initialized history. That is, unique weight parameters per particle history.
This follows the practice in the grand-daddy of Monte Carlos
https://apps.dtic.mil/sti/tr/pdf/ADA231425.pdf
doc page 64/65, pdf page 70/71
code block 375

Here is the best reference and investigation into weight cutoffs (survival biasing parameters) that I have come across.
https://www.osti.gov/servlets/purl/5823239
doc page 87, pdf page 96
paragraphs + Table 2.2

Example

An example application of Survival Biasing with Survival Normalization can be found here:
https://github.com/yrrepy/OMC-SurvNorm_Example
The .xlsx in that repository describes the problem and contains the results.

Case 1 is an SSW/SSR where the user changed the weight parameters between the two calculations. Survival normalization improves the FOM by a minimum of +12% in the outer shell, the most difficult region to score in this case. FOMs elsewhere decrease.

Case 2&3 is a constructed MCPL phase-space source masquerading as a biased isotropic point source.
Survival normalization almost across the board significantly improves FOMs
Case 2&3 Xb use the Survival Normalization feature. Case 2&3 Xc uses a manual multiplication of the weight parameters by the minimum weight of the biased source. The in-code per history Survival Normalization produces better FOMs than with the singular adjustment.
It could well be that the current phase-space condition in the code could be removed and allow Survival Normalization to work for all sources (under the hypothesis that someone produced a biased library source). But whether the per history FOM improvement continued in general I am not sure and could well be not the case.

In term one these examples could be the basis for additional test(s).

Code Changes

We added a few bits of functionality into the settings.cpp, particle_data.h, simulation.cpp and connecting physics files. The code takes the xml Boolean and enables the functionality of survival normalization. The values inside the survival normalization are then stored inside new object variables inside particle class. These values are then used inside all the physics files for calculations.

Settings.cpp

Added survival_normalization and source_file Booleans
Added source_file boolean inside the source xml toggle
Added new xml toggle boolean to enable and disable survival_normalization

Particle_data.h

Added private variables inside the particle class
Added new object variables wgt_cutoff and wgt_survive
Added setters and getters for both new object variables

Simulation.cpp

Toggle to adjust weight cutoff and weight survive by multiplying the current weight
Added the survival normalization functionality inside the given if statements scenario from settings.cpp

Physics.cpp/physics_common.cpp/physics_mg.cpp

Add the object variable into the appropriate weight_cutoff and weight_survive given the if statement scenario

Fixes

This PR seeks to address in part:
https://openmc.discourse.group/t/tally-uncertainty-openmc-vs-mcnp/1872/11

Checklist

  • I have performed a self-review of my own code
  • I have run clang-format on any C++ source files (if applicable)
  • I have followed the style guidelines for Python source files (if applicable)
  • I have made corresponding changes to the documentation (if applicable)
  • I have added tests that prove my fix is effective or that my feature works (if applicable)

To-Do

  • Implement Python API functionality

quantaroo and others added 24 commits August 10, 2023 13:07
Description of Boolean enabling Survival Normalization.

Description briefly of functionality and current limitations.
Described Survival Normalization for phase-space sources. Unfortunately there is little description or justification for this out in the general topic literature. Though I have found it to work better in multiple test cases, notably the FOM in areas with middling statistics (~10%).

In term should add capability of normalization for biased general sources.
Code commenting for actual normalization mechanism
Code commenting
Remove duplicate description of Surv Norm
Survival Normalization by Default Off
Update text on Survival Normalization, default=false
Better explaining when to use survival normalization
@yrrepy yrrepy requested a review from nelsonag as a code owner August 29, 2023 20:23
@paulromano paulromano mentioned this pull request Sep 1, 2023
5 tasks
Copy link
Contributor

@eepeterson eepeterson left a comment

Choose a reason for hiding this comment

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

Thanks @yrrepy for the PR! We do need to implement general source biasing relatively soon, so I would say we might want to rework this in that direction a little bit. Namely, I don't think it would be too difficult to make this work with general sources. My recommendations from a first round of comments are below:

  • Add just the initial source weight to the particle class rather than two other variables
  • multiply relevant weight_survive and weight_cutoff values by initial weight
  • This will have to also be updated in our weight windows implementation

I would think this would make the implementation cleaner, make it ready for when we introduce biased sources, and have minimal performance drawbacks, but curious to hear other peoples opinions. @paulromano @pshriwise

include/openmc/particle_data.h Outdated Show resolved Hide resolved
src/physics.cpp Outdated Show resolved Hide resolved
src/physics_common.cpp Outdated Show resolved Hide resolved
src/physics_mg.cpp Outdated Show resolved Hide resolved
Single variable for Survival Normalization
Use single variable wgt0, and multiply cutoffs here (on the fly), 

rather than having to multply in simulation.cpp.
and then telling the physics modules to get this or that cutoff.
Revert to original, now that cutoffs were multiplied on the fly
Revert with simplified survival normalization
Assign Particle Start Weight
Punctuation
Remove lingering cuttoff variables in particle_data.h
Bracket typo in physics.cpp
Not entering final if statement (russian roulette for not survival normalization)
multiplication by wgt0 if Survival Normalization on
Pedantic capitalization to match physics_mg.cpp
@yrrepy
Copy link
Contributor Author

yrrepy commented Sep 30, 2023

Hi Ethan,
thank you for the review. This is my first time formally PR'ing, so apologies if I don't get all the workflow correct.

Indeed, I find your suggestion to simplify the mechanics of the normalization much cleaner. Having the variable wgt0 in the particle class is arguably much more useful than normalized cutoffs in the class.

I left the if statement in physicsX.cpp for the case where a source is biased (or of varying weights) but the user does not want to normalize. But perhaps there is a way to have the code check instantly there when it multiplies by wgt0?, I'm not terribly familiar with the fancier functions in C++.
The if survival_normalization=false could go in simulation.cpp and tell it to not set wgt0. But seems better to just always set wgt0.


Re: General Sources
I did some testing and in my two simple cases found that the per particle survival normalization was better than normalization by the minimum weight of all the particles.
The current condition that the normalization mechanics only work for phase spaces,
if ((settings::source_file || settings::surf_source_read),
could be removed so that this mechanism works for general sources. But I would recommend testing it on some purpose-built biased sources first.

Reorganize and reorder Boolean flags to account for changes that have occurred in the interim with the recent addition of weight_window_checkpoint_surface and weight_window_checkpoint_collision.

Thus Resolving a conflict in the present PR.
@yrrepy
Copy link
Contributor Author

yrrepy commented Oct 13, 2023

What should I focus on the keep this moving forward?

Some tests? Python bindings?

@paulromano paulromano added this to the v0.14.0 milestone Oct 19, 2023
@paulromano
Copy link
Contributor

@eepeterson @pshriwise and I talked about this PR the other day. I'll let @eepeterson summarize our thoughts as we were all on the same page.

Copy link
Contributor

@eepeterson eepeterson left a comment

Choose a reason for hiding this comment

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

@yrrepy sorry for the delays in the review. I think we can accomplish what you want with minimal changes to the existing code and no new settings or runtime parameters aside from setting the initial particle weight as previously discussed. That way we can rely on the existing test suite to make sure we haven't broken anything and since we haven't introduced any new "features" shouldn't need to write additional tests.

include/openmc/particle_data.h Outdated Show resolved Hide resolved
include/openmc/settings.h Outdated Show resolved Hide resolved
@@ -53,14 +54,19 @@ extern bool surf_source_write; //!< write surface source file?
extern bool surf_mcpl_write; //!< write surface mcpl file?
extern bool surf_source_read; //!< read surface source file?
extern bool survival_biasing; //!< use survival biasing?
extern bool survival_normalization;//!< use survival normalization?
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
extern bool survival_normalization;//!< use survival normalization?

extern bool temperature_multipole; //!< use multipole data?
extern "C" bool trigger_on; //!< tally triggers enabled?
extern bool trigger_predict; //!< predict batches for triggers?
extern bool ufs_on; //!< uniform fission site method on?
extern bool urr_ptables_on; //!< use unresolved resonance prob. tables?
extern "C" bool weight_windows_on; //!< are weight windows are enabled?
extern bool write_all_tracks; //!< write track files for every particle?
extern bool write_initial_source; //!< write out initial source file?
extern bool weight_window_checkpoint_surface; //!< enable weight window check
Copy link
Contributor

Choose a reason for hiding this comment

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

I think these are from another PR - if you rebase on develop these should go away.

src/physics.cpp Outdated
Comment on lines 160 to 166
// if survival normalization is applicable, use normalized weight cutoff and normalized weight survive
if ((settings::source_file || settings::surf_source_read)&&(settings::survival_normalization)) {
if (p.wgt() < settings::weight_cutoff*p.wgt0()) {
russian_roulette(p, settings::weight_survive*p.wgt0());
}
} else if (p.wgt() < settings::weight_cutoff) {
russian_roulette(p, settings::weight_survive);
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
// if survival normalization is applicable, use normalized weight cutoff and normalized weight survive
if ((settings::source_file || settings::surf_source_read)&&(settings::survival_normalization)) {
if (p.wgt() < settings::weight_cutoff*p.wgt0()) {
russian_roulette(p, settings::weight_survive*p.wgt0());
}
} else if (p.wgt() < settings::weight_cutoff) {
russian_roulette(p, settings::weight_survive);
if (p.wgt() < settings::weight_cutoff*p.wgt0()) {
russian_roulette(p, settings::weight_survive*p.wgt0());

I think we can apply this globally for all cases rather than needing additional settings and checks.

Comment on lines 67 to 73
// if survival normalization is applicable, use normalized weight cutoff and normalized weight survive
if ((settings::source_file || settings::surf_source_read)&&(settings::survival_normalization)) {
if (p.wgt() < settings::weight_cutoff*p.wgt0()) {
russian_roulette(p, settings::weight_survive*p.wgt0());
}
} else if (p.wgt() < settings::weight_cutoff) {
russian_roulette(p, settings::weight_survive);
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
// if survival normalization is applicable, use normalized weight cutoff and normalized weight survive
if ((settings::source_file || settings::surf_source_read)&&(settings::survival_normalization)) {
if (p.wgt() < settings::weight_cutoff*p.wgt0()) {
russian_roulette(p, settings::weight_survive*p.wgt0());
}
} else if (p.wgt() < settings::weight_cutoff) {
russian_roulette(p, settings::weight_survive);
if (p.wgt() < settings::weight_cutoff*p.wgt0()) {
russian_roulette(p, settings::weight_survive*p.wgt0());

same comment here as above.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

fine with removing SSR gating, but I do think it is best to be able to turn this feature off

src/settings.cpp Outdated
Copy link
Contributor

Choose a reason for hiding this comment

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

I think we can accomplish the point of this PR with no additional settings, so my recommendation would be to revert those everywhere (settings.cpp, settings.h, etc)

p.current_work() = index_source;

// set identifier for particle
p.id() = simulation::work_index[mpi::rank] + index_source;

// set particle history start weight
p.wgt0(p.wgt());
Copy link
Contributor

Choose a reason for hiding this comment

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

We should set the initial particle weight in the Particle::from_source method instead of here in the same fashion that wgt and wgt_last are set.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I tried what was mentioned, but it didn't work:
p.wgt0(p.wgt());

The legacy one does:
p.wgt0() = p.wgt();

Comment on lines 1681 to 1698
.. _survival_normalization:

Survival Normalization
----------------

This parameter can be used to make the adjustable parameters of survival biasing,
:math:`w_c` and :math:`w_s`, relative to the starting weight of a source particle.

It is currently only implemented for phase-space sources (MCPL file sources
and HDF5 surface sources). For these, the normalization of the parameters
is done per source particle. That is, for each history, the parameters
:math:`w_c` and :math:`w_s`, are multiplied by the start weight of the current
history.

This normalization is recommended for problems where numerous source particles are
initialized below the survival biasing weight cutoff :math:`w_c`; to prevent them
from being immediately rouletted. This can notably be the case for biased sources.

Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
.. _survival_normalization:
Survival Normalization
----------------
This parameter can be used to make the adjustable parameters of survival biasing,
:math:`w_c` and :math:`w_s`, relative to the starting weight of a source particle.
It is currently only implemented for phase-space sources (MCPL file sources
and HDF5 surface sources). For these, the normalization of the parameters
is done per source particle. That is, for each history, the parameters
:math:`w_c` and :math:`w_s`, are multiplied by the start weight of the current
history.
This normalization is recommended for problems where numerous source particles are
initialized below the survival biasing weight cutoff :math:`w_c`; to prevent them
from being immediately rouletted. This can notably be the case for biased sources.

We don't need a new section here, but in the existing survival biasing section it could be worth mentioning the the weight_cutoff and weight_survive used for russian roulette are the respective fractions of the initial particle weight.

@yrrepy
Copy link
Contributor Author

yrrepy commented Oct 27, 2023

Hi,
well glad you guys think this kind of normalization will work in general. From what I looked at, it did seem effective.

I'm happy to remove the Boolean conditions that were limiting it to only phase-spaces, but I do think it would be best to keep the survival_normalization Boolean allowing it to be switched on and off (a-la MCNP where wc1,wc2 can be <0 or >0 determining whether to norm, but here with the Boolean)
I did find some situations where it was better to not normalize the weight cutoffs.

@paulromano paulromano removed this from the v0.14.0 milestone Nov 2, 2023
@yrrepy
Copy link
Contributor Author

yrrepy commented Dec 8, 2023

Bump on the above, and comments on the approach to take for a Boolean switch
(before coding it and the other suggested changes)

yrrepy and others added 5 commits June 25, 2024 14:38
no boolean gating off survival normalization for SSR only

Co-authored-by: Ethan Peterson <[email protected]>
no SSR gating
remover "relevant for only SSR"
@yrrepy
Copy link
Contributor Author

yrrepy commented Jul 4, 2024

Closed in favor of #3070 which is synced with latest OpenMC version

@yrrepy yrrepy closed this Jul 4, 2024
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.

4 participants