Resolve conflict with weight windows and global russian roulette

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
 //==============================================================================
@@ -135,10 +124,19 @@ class WeightWindows {
   //! \param[in] group  HDF5 group to write to
   void to_hdf5(hid_t group) const;
 
+  //! Get particle bin inside weight window mesh
+  //! \param[in] p  Particle to get weight window bin for
+  const int get_mesh_bin(const Particle& p) const;
+
   //! Retrieve the weight window for a particle
   //! \param[in] p  Particle to get weight window for
   WeightWindow get_weight_window(const Particle& p) const;
 
+  //! Retrieve the weight window for a mesh bin
+  //! \param[in] double E particle energy
+  //! \param[in] int mesh_bin mesh bin
+  WeightWindow get_weight_window(double E, const int mesh_bin) const;
+
   std::array<int, 2> bounds_size() const;
 
   const vector<double>& energy_bounds() const { return energy_bounds_; }
@@ -237,6 +235,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
+WeightWindow search_weight_window(const Particle& p);
+
 //! Finalize variance reduction objects after all inputs have been read
 void finalize_variance_reduction();
 

src/physics.cpp

@@ -61,8 +61,25 @@ void collision(Particle& p)
     break;
   }
 
-  if (settings::weight_window_checkpoint_collision)
-    apply_weight_windows(p);
+  if (settings::weight_windows_on) {
+    if (settings::weight_window_checkpoint_collision) {
+      // If collision checkpoint is enabled, apply weight window
+      // if valid and apply global russian roulette if not.
+      auto ww = search_weight_window(p);
+      if (ww.is_valid()) {
+        apply_weight_window(p, ww);
+      } else if (p.type() == ParticleType::neutron) {
+        apply_russian_roulette(p);
+      }
+    } else if (p.type() == ParticleType::neutron) {
+      // If collision checkpoint is disabled and weight windows are enabled,
+      // apply russian roulette if the particle is a neutron and it is outside
+      // the weight window domain
+      auto ww = search_weight_window(p);
+      if (!ww.is_valid())
+        apply_russian_roulette(p);
+    }
+  }
 
   // Kill particle if energy falls below cutoff
   int type = static_cast<int>(p.type());
@@ -153,18 +170,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

@@ -31,8 +31,15 @@ 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 = search_weight_window(p);
+    if (ww.is_valid()) {
+      if (settings::weight_window_checkpoint_collision)
+        apply_weight_window(p, ww);
+    } else {
+      apply_russian_roulette(p);
+    }
+  }
 
   // Display information about collision
   if ((settings::verbosity >= 10) || p.trace()) {
@@ -66,18 +73,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/settings.cpp

@@ -1264,6 +1264,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");