Resolve conflict with weight windows and global russian roulette (#3751)

Co-authored-by: Patrick Shriwise <pshriwise@gmail.com>

Author: GuySten

include/openmc/physics_common.h

@@ -13,5 +13,9 @@ namespace openmc {
 //! \param[in] weight_survive Weight assigned to particles that survive
 void russian_roulette(Particle& p, double weight_survive);
 
+//! \brief Performs the global russian roulette operation on a particle
+//! \param[in,out] p  Particle object
+void apply_russian_roulette(Particle& p);
+
 } // namespace openmc
 #endif // OPENMC_PHYSICS_COMMON_H

include/openmc/weight_windows.h

@@ -25,17 +25,6 @@ enum class WeightWindowUpdateMethod { MAGIC, FW_CADIS };
 
 constexpr double DEFAULT_WEIGHT_CUTOFF {1.0e-38}; // default low weight cutoff
 
-//==============================================================================
-// 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();
-
 //==============================================================================
 // Global variables
 //==============================================================================
@@ -137,7 +126,7 @@ class WeightWindows {
 
   //! Retrieve the weight window for a particle
   //! \param[in] p  Particle to get weight window for
-  WeightWindow get_weight_window(const Particle& p) const;
+  std::pair<bool, WeightWindow> get_weight_window(const Particle& p) const;
 
   std::array<int, 2> bounds_size() const;
 
@@ -236,6 +225,26 @@ class WeightWindowsGenerator {
   double ratio_ {5.0};     //<! ratio of lower to upper weight window bounds
 };
 
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+//! Apply weight windows to a particle
+//! \param[in] p  Particle to apply weight windows to
+void apply_weight_windows(Particle& p);
+
+//! Apply weight window to a particle
+//! \param[in] p  Particle to apply weight window to
+//! \param[in] weight_window WeightWindow to apply
+void apply_weight_window(Particle& p, WeightWindow weight_window);
+
+//! Free memory associated with weight windows
+void free_memory_weight_windows();
+
+//! Search weight window that apply to a particle
+//! \param[in]  p  Particle to search weight window for
+std::pair<bool, WeightWindow> search_weight_window(const Particle& p);
+
 //! Finalize variance reduction objects after all inputs have been read
 void finalize_variance_reduction();
 

src/physics.cpp

@@ -64,8 +64,17 @@ void collision(Particle& p)
     fatal_error("Unsupported particle PDG for collision sampling.");
   }
 
-  if (settings::weight_window_checkpoint_collision)
-    apply_weight_windows(p);
+  if (settings::weight_windows_on) {
+    auto [ww_found, ww] = search_weight_window(p);
+    if (!ww_found && p.type() == ParticleType::neutron()) {
+      // if the weight window is not valid, apply russian roulette for neutrons
+      // (regardless of weight window collision checkpoint setting)
+      apply_russian_roulette(p);
+    } else if (settings::weight_window_checkpoint_collision) {
+      // if collision checkpointing is on, apply weight window
+      apply_weight_window(p, ww);
+    }
+  }
 
   // Kill particle if energy falls below cutoff
   int type = p.type().transport_index();
@@ -156,18 +165,9 @@ void sample_neutron_reaction(Particle& p)
     advance_prn_seed(data::nuclides.size(), &p.seeds(STREAM_URR_PTABLE));
   }
 
-  // Play russian roulette if survival biasing is turned on
-  if (settings::survival_biasing) {
-    // if survival normalization is on, use normalized weight cutoff and
-    // normalized weight survive
-    if (settings::survival_normalization) {
-      if (p.wgt() < settings::weight_cutoff * p.wgt_born()) {
-        russian_roulette(p, settings::weight_survive * p.wgt_born());
-      }
-    } else if (p.wgt() < settings::weight_cutoff) {
-      russian_roulette(p, settings::weight_survive);
-    }
-  }
+  // Play russian roulette if there are no weight windows
+  if (!settings::weight_windows_on)
+    apply_russian_roulette(p);
 }
 
 void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)

src/physics_common.cpp

@@ -18,4 +18,20 @@ void russian_roulette(Particle& p, double weight_survive)
   }
 }
 
+void apply_russian_roulette(Particle& p)
+{
+  // Exit if survival biasing is turned off
+  if (!settings::survival_biasing)
+    return;
+
+  // if survival normalization is on, use normalized weight cutoff and
+  // normalized weight survive
+  if (settings::survival_normalization) {
+    if (p.wgt() < settings::weight_cutoff * p.wgt_born()) {
+      russian_roulette(p, settings::weight_survive * p.wgt_born());
+    }
+  } else if (p.wgt() < settings::weight_cutoff) {
+    russian_roulette(p, settings::weight_survive);
+  }
+}
 } // namespace openmc

src/physics_mg.cpp

@@ -32,8 +32,17 @@ void collision_mg(Particle& p)
   // Sample the reaction type
   sample_reaction(p);
 
-  if (settings::weight_window_checkpoint_collision)
-    apply_weight_windows(p);
+  if (settings::weight_windows_on) {
+    auto [ww_found, ww] = search_weight_window(p);
+    if (!ww_found && p.type() == ParticleType::neutron()) {
+      // if the weight window is not valid, apply russian roulette
+      // (regardless of weight window collision checkpoint setting)
+      apply_russian_roulette(p);
+    } else if (settings::weight_window_checkpoint_collision) {
+      // if collision checkpointing is on, apply weight window
+      apply_weight_window(p, ww);
+    }
+  }
 
   // Display information about collision
   if ((settings::verbosity >= 10) || p.trace()) {
@@ -67,18 +76,9 @@ void sample_reaction(Particle& p)
   // Sample a scattering event to determine the energy of the exiting neutron
   scatter(p);
 
-  // Play Russian roulette if survival biasing is turned on
-  if (settings::survival_biasing) {
-    // if survival normalization is applicable, use normalized weight cutoff and
-    // normalized weight survive
-    if (settings::survival_normalization) {
-      if (p.wgt() < settings::weight_cutoff * p.wgt_born()) {
-        russian_roulette(p, settings::weight_survive * p.wgt_born());
-      }
-    } else if (p.wgt() < settings::weight_cutoff) {
-      russian_roulette(p, settings::weight_survive);
-    }
-  }
+  // Play russian roulette if there are no weight windows
+  if (!settings::weight_windows_on)
+    apply_russian_roulette(p);
 }
 
 void scatter(Particle& p)

src/random_ray/flat_source_domain.cpp

@@ -843,7 +843,7 @@ void FlatSourceDomain::output_to_vtk() const
           voxel_positions[z * Ny * Nx + y * Nx + x] = sample;
 
           if (variance_reduction::weight_windows.size() == 1) {
-            WeightWindow ww =
+            auto [ww_found, ww] =
               variance_reduction::weight_windows[0]->get_weight_window(p);
             float weight = ww.lower_weight;
             weight_windows[z * Ny * Nx + y * Nx + x] = weight;

src/settings.cpp

@@ -1266,6 +1266,13 @@ void read_settings_xml(pugi::xml_node root)
     }
   }
 
+  if (weight_windows_on) {
+    if (!weight_window_checkpoint_surface &&
+        !weight_window_checkpoint_collision)
+      fatal_error(
+        "Weight Windows are enabled but there are no valid checkpoints.");
+  }
+
   if (check_for_node(root, "use_decay_photons")) {
     settings::use_decay_photons =
       get_node_value_bool(root, "use_decay_photons");