From ac495ef0f332019f92e906e8f69e0b12d63f5ff7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 8 Jan 2022 15:26:02 -0600 Subject: [PATCH] Move main weight window function to weight_windows.cpp --- include/openmc/weight_windows.h | 4 ++ src/physics.cpp | 69 +-------------------------------- src/weight_windows.cpp | 64 ++++++++++++++++++++++++++++++ 3 files changed, 70 insertions(+), 67 deletions(-) diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 7a7b2177f..12059f1cb 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -19,6 +19,10 @@ namespace openmc { // Non-member functions //============================================================================== +//! Apply weight windows to a particle +//! \param[in] p Particle to apply weight windows to +void apply_weight_windows(Particle& p); + //! Free memory associated with weight windows void free_memory_weight_windows(); diff --git a/src/physics.cpp b/src/physics.cpp index a2190698f..c8f1131f8 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -46,13 +46,13 @@ void collision(Particle& p) case ParticleType::neutron: sample_neutron_reaction(p); if (settings::weight_windows_on) { - split_particle(p); + apply_weight_windows(p); } break; case ParticleType::photon: sample_photon_reaction(p); if (settings::weight_windows_on) { - split_particle(p); + apply_weight_windows(p); } break; case ParticleType::electron: @@ -1204,69 +1204,4 @@ void sample_secondary_photons(Particle& p, int i_nuclide) } } -// split a particle based on the statistical weight of the current location -void split_particle(Particle& p) -{ - // skip dead or no energy - if (p.E() <= 0 || !p.alive()) - return; - - bool in_domain = false; - // TODO: this is a linear search - should do something more clever - WeightWindow weight_window; - for (const auto& ww : variance_reduction::weight_windows) { - weight_window = ww->get_weight_window(p); - if (weight_window) - break; - } - // particle is not in any of the ww domains, do nothing - if (!weight_window) - return; - - // get the paramters - double weight = p.wgt(); - p.wgt_last() = weight; - - // first check to see if particle should be killed for weight cutoff - if (p.wgt() < weight_window.weight_cutoff) { - p.alive() = false; - p.wgt() = 0.0; - return; - } - - // if particle's weight is above the weight window split until they are within - // the window - if (weight > weight_window.upper_weight) { - // do not further split the particle if above the limit - if (p.n_split() >= settings::max_splits) - return; - - double n_split = std::ceil(weight / weight_window.upper_weight); - double max_split = weight_window.max_split; - n_split = std::min(n_split, max_split); - - p.n_split() += n_split; - - // Create secondaries and divide weight among all particles - int i_split = std::round(n_split); - for (int l = 0; l < i_split - 1; l++) { - p.create_secondary(weight / n_split, p.u(), p.E(), p.type()); - } - // TODO: maybe weight should be weight - sum of child weight - p.wgt() = weight / n_split; - - } else if (weight <= weight_window.lower_weight) { - // if the particle weight is below the window, play Russian roulette - double weight_survive = - std::min(weight * weight_window.max_split, weight_window.survival_weight); - if (weight_survive * prn(p.current_seed()) <= weight) { - p.wgt() = weight_survive; - } else { - p.alive() = false; - p.wgt() = 0.0; - } - // else particle is in the window, continue as normal - } -} - } // namespace openmc diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 68b2d25e8..b62725d28 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -28,6 +28,70 @@ openmc::vector> weight_windows; // Non-member functions //============================================================================== +void apply_weight_windows(Particle& p) +{ + // skip dead or no energy + if (p.E() <= 0 || !p.alive()) + return; + + bool in_domain = false; + // TODO: this is a linear search - should do something more clever + WeightWindow weight_window; + for (const auto& ww : variance_reduction::weight_windows) { + weight_window = ww->get_weight_window(p); + if (weight_window) + break; + } + // particle is not in any of the ww domains, do nothing + if (!weight_window) + return; + + // get the paramters + double weight = p.wgt(); + p.wgt_last() = weight; + + // first check to see if particle should be killed for weight cutoff + if (p.wgt() < weight_window.weight_cutoff) { + p.alive() = false; + p.wgt() = 0.0; + return; + } + + // if particle's weight is above the weight window split until they are within + // the window + if (weight > weight_window.upper_weight) { + // do not further split the particle if above the limit + if (p.n_split() >= settings::max_splits) + return; + + double n_split = std::ceil(weight / weight_window.upper_weight); + double max_split = weight_window.max_split; + n_split = std::min(n_split, max_split); + + p.n_split() += n_split; + + // Create secondaries and divide weight among all particles + int i_split = std::round(n_split); + for (int l = 0; l < i_split - 1; l++) { + p.create_secondary(weight / n_split, p.u(), p.E(), p.type()); + } + // TODO: maybe weight should be weight - sum of child weight + p.wgt() = weight / n_split; + + } else if (weight <= weight_window.lower_weight) { + // if the particle weight is below the window, play Russian roulette + double weight_survive = + std::min(weight * weight_window.max_split, weight_window.survival_weight); + if (weight_survive * prn(p.current_seed()) <= weight) { + p.wgt() = weight_survive; + } else { + p.alive() = false; + p.wgt() = 0.0; + } + // else particle is in the window, continue as normal + } +} + void free_memory_weight_windows() { variance_reduction::ww_map.clear();