Merge pull request #1956 from pshriwise/otf_ww_adj

On-the-fly adjustments for weight windows

docs/source/io_formats/settings.rst

@@ -1050,6 +1050,12 @@ sub-elements/attributes:
 
     *Default*: 3.0
 
+  :max_lower_bound_ratio:
+    Maximum allowed ratio of a particle's weight to the weight window's lower
+    bound. A factor will be applied to raise the weight window to be lower than
+    the particle's weight by a factor of max_lower_bound_ratio during transport
+    if exceeded.
+
   :max_split:
     Maximum allowable number of particles when splitting
 

include/openmc/particle_data.h

@@ -309,7 +309,8 @@ private:
   int n_event_ {0}; // number of events executed in this particle's history
 
   // Weight window information
-  int n_split_ {0}; // Number of splits this particle has undergone
+  int n_split_ {0}; // Number of times this particle has been split
+  double ww_factor_ {0.0}; // Particle-specific factor for on-the-fly weight window adjustment
 
 // DagMC state variables
 #ifdef DAGMC
@@ -435,6 +436,9 @@ public:
   int n_split() const { return n_split_; }
   int& n_split() { return n_split_; }
 
+  double ww_factor() const { return ww_factor_; }
+  double& ww_factor() { return ww_factor_; }
+
 #ifdef DAGMC
   moab::DagMC::RayHistory& history() { return history_; }
   Direction& last_dir() { return last_dir_; }

include/openmc/weight_windows.h

@@ -52,12 +52,20 @@ extern vector<unique_ptr<WeightWindows>> weight_windows;
 struct WeightWindow {
   double lower_weight {-1}; // -1 indicates invalid state
   double upper_weight {1};
+  double max_lb_ratio {1};
   double survival_weight {0.5};
   double weight_cutoff {DEFAULT_WEIGHT_CUTOFF};
   int max_split {1};
 
   //! Whether the weight window is in a valid state
   bool is_valid() const { return lower_weight >= 0.0; }
+
+  //! Adjust the weight window by a constant factor
+  void scale(double factor)
+  {
+    lower_weight *= factor;
+    upper_weight *= factor;
+  }
 };
 
 //==============================================================================
@@ -96,6 +104,7 @@ private:
   vector<double> lower_ww_;        //!< Lower weight window bounds
   vector<double> upper_ww_;        //!< Upper weight window bounds
   double survival_ratio_ {3.0};    //!< Survival weight ratio
+  double max_lb_ratio_ {1.0}; //!< Maximum lower bound to particle weight ratio
   double weight_cutoff_ {DEFAULT_WEIGHT_CUTOFF}; //!< Weight cutoff
   int max_split_ {10};             //!< Maximum value for particle splitting
   int32_t mesh_idx_;               //!< index in meshes vector

openmc/weight_windows.py

@@ -48,6 +48,11 @@ class WeightWindows(IDManagerMixin):
     survival_ratio : float
         Ratio of the survival weight to the lower weight window bound for
         rouletting
+    max_lower_bound_ratio : float
+        Maximum allowed ratio of a particle's weight to the weight window's
+        lower bound. A factor will be applied to raise the weight window to be lower
+        than the particle's weight by a factor of max_lower_bound_ratio
+        during transport if exceeded.
     max_split : int
         Maximum allowable number of particles when splitting
     weight_cutoff : float
@@ -76,6 +81,9 @@ class WeightWindows(IDManagerMixin):
     survival_ratio : float
         Ratio of the survival weight to the lower weight window bound for
         rouletting
+    max_lower_bound_ratio: float
+        Maximum allowed ratio of a particle's weight to the weight window's
+        lower bound. (Default: 1.0)
     max_split : int
         Maximum allowable number of particles when splitting
     weight_cutoff : float
@@ -91,7 +99,8 @@ class WeightWindows(IDManagerMixin):
 
     def __init__(self, mesh, lower_ww_bounds, upper_ww_bounds=None,
                  upper_bound_ratio=None, energy_bins=None, particle_type='neutron',
-                 survival_ratio=3, max_split=10, weight_cutoff=1.e-38, id=None):
+                 survival_ratio=3, max_lower_bound_ratio=None, max_split=10,
+                 weight_cutoff=1.e-38, id=None):
         self.mesh = mesh
         self.id = id
         self.particle_type = particle_type
@@ -121,6 +130,11 @@ class WeightWindows(IDManagerMixin):
                              'do not match')
 
         self.survival_ratio = survival_ratio
+
+        self._max_lower_bound_ratio = None
+        if max_lower_bound_ratio is not None:
+            self.max_lower_bound_ratio = max_lower_bound_ratio
+
         self.max_split = max_split
         self.weight_cutoff = weight_cutoff
 
@@ -192,6 +206,16 @@ class WeightWindows(IDManagerMixin):
         cv.check_greater_than('Survival ratio', val, 1.0, True)
         self._survival_ratio = val
 
+    @property
+    def max_lower_bound_ratio(self):
+        return self._max_lower_bound_ratio
+
+    @max_lower_bound_ratio.setter
+    def max_lower_bound_ratio(self, val):
+        cv.check_type('Maximum lower bound ratio', val, Real)
+        cv.check_greater_than('Maximum lower bound ratio', val, 1.0)
+        self._max_lower_bound_ratio = val
+
     @property
     def max_split(self):
         return self._max_split
@@ -241,6 +265,10 @@ class WeightWindows(IDManagerMixin):
         subelement = ET.SubElement(element, 'survival_ratio')
         subelement.text = str(self.survival_ratio)
 
+        if self.max_lower_bound_ratio is not None:
+            subelement = ET.SubElement(element, 'max_lower_bound_ratio')
+            subelement.text = str(self.max_lower_bound_ratio)
+
         subelement = ET.SubElement(element, 'max_split')
         subelement.text = str(self.max_split)
 
@@ -278,6 +306,11 @@ class WeightWindows(IDManagerMixin):
         ebins = [float(b) for b in get_text(elem, 'energy_bins').split()]
         particle_type = get_text(elem, 'particle_type')
         survival_ratio = float(get_text(elem, 'survival_ratio'))
+
+        max_lower_bound_ratio = None
+        if get_text(elem, 'max_lower_bound_ratio'):
+            max_lower_bound_ratio = float(get_text(elem, 'max_lower_bound_ratio'))
+
         max_split = int(get_text(elem, 'max_split'))
         weight_cutoff = float(get_text(elem, 'weight_cutoff'))
         id = int(get_text(elem, 'id'))
@@ -289,6 +322,7 @@ class WeightWindows(IDManagerMixin):
             energy_bins=ebins,
             particle_type=particle_type,
             survival_ratio=survival_ratio,
+            max_lower_bound_ratio=max_lower_bound_ratio,
             max_split=max_split,
             weight_cutoff=weight_cutoff,
             id=id
@@ -318,6 +352,11 @@ class WeightWindows(IDManagerMixin):
         lower_ww_bounds = group['lower_ww_bounds'][()]
         upper_ww_bounds = group['upper_ww_bounds'][()]
         survival_ratio = group['survival_ratio'][()]
+
+        max_lower_bound_ratio = None
+        if group.get('max_lower_bound_ratio') is not None:
+            max_lower_bound_ratio = group['max_lower_bound_ratio'][()]
+
         max_split = group['max_split'][()]
         weight_cutoff = group['weight_cutoff'][()]
 
@@ -328,6 +367,7 @@ class WeightWindows(IDManagerMixin):
             energy_bins=ebins,
             particle_type=ptype,
             survival_ratio=survival_ratio,
+            max_lower_bound_ratio=max_lower_bound_ratio,
             max_split=max_split,
             weight_cutoff=weight_cutoff,
             id=id

src/simulation.cpp

@@ -490,6 +490,9 @@ void initialize_history(Particle& p, int64_t index_source)
   // Reset split counter
   p.n_split() = 0;
 
+  // Reset weight window ratio
+  p.ww_factor() = 0.0;
+
   // set random number seed
   int64_t particle_seed =
     (simulation::total_gen + overall_generation() - 1) * settings::n_particles +

src/weight_windows.cpp

@@ -56,6 +56,19 @@ void apply_weight_windows(Particle& p)
     return;
   }
 
+  // check if particle is far above current weight window
+  // only do this if the factor is not already set on the particle and a
+  // maximum lower bound ratio is specified
+  if (p.ww_factor() == 0.0 && weight_window.max_lb_ratio > 1.0 &&
+      p.wgt() > weight_window.lower_weight * weight_window.max_lb_ratio) {
+    p.ww_factor() =
+      p.wgt() / (weight_window.lower_weight * weight_window.max_lb_ratio);
+  }
+
+  // move weight window closer to the particle weight if needed
+  if (p.ww_factor() > 1.0)
+    weight_window.scale(p.ww_factor());
+
   // if particle's weight is above the weight window split until they are within
   // the window
   if (weight > weight_window.upper_weight) {
@@ -139,6 +152,14 @@ WeightWindows::WeightWindows(pugi::xml_node node)
                   "and less than the upper to lower weight window ratio.");
   }
 
+  // get the max lower bound ratio - optional
+  if (check_for_node(node, "max_lower_bound_ratio")) {
+    max_lb_ratio_ = std::stod(get_node_value(node, "max_lower_bound_ratio"));
+    if (max_lb_ratio_ < 1.0) {
+      fatal_error("Maximum lower bound ratio must be larger than 1");
+    }
+  }
+
   // get the max split - optional
   if (check_for_node(node, "max_split")) {
     max_split_ = std::stod(get_node_value(node, "max_split"));
@@ -239,6 +260,7 @@ WeightWindow WeightWindows::get_weight_window(const Particle& p) const
   ww.lower_weight = lower_ww_[ww_index];
   ww.upper_weight = upper_ww_[ww_index];
   ww.survival_weight = ww.lower_weight * survival_ratio_;
+  ww.max_lb_ratio = max_lb_ratio_;
   ww.max_split = max_split_;
   ww.weight_cutoff = weight_cutoff_;
   return ww;
@@ -254,6 +276,7 @@ void WeightWindows::to_hdf5(hid_t group) const
   write_dataset(ww_group, "lower_ww_bounds", lower_ww_);
   write_dataset(ww_group, "upper_ww_bounds", upper_ww_);
   write_dataset(ww_group, "survival_ratio", survival_ratio_);
+  write_dataset(ww_group, "max_lower_bound_ratio", max_lb_ratio_);
   write_dataset(ww_group, "max_split", max_split_);
   write_dataset(ww_group, "weight_cutoff", weight_cutoff_);
   write_dataset(ww_group, "mesh", this->mesh().id_);

tests/regression_tests/weightwindows/results_true.dat

@@ -1 +1 @@
-c77af605fb23bc3f03bb304b82cc2d0d9efe7826167f43202fa0ecba17e6a77fc1a91bf9f4b23f175ae98c04c69067bffb96d7d68e30a126f5fa7c0080ea8367
+f10c722e27d1f0a69f700bc72c4b7751f375dc0a74e049c9abb4b261d2265155c0e516e7e38f9751b83a24eac07e944e51740d398b720524cf8cd6bf0b8c51fc

tests/regression_tests/weightwindows/test.py

@@ -92,13 +92,15 @@ def model():
                                 ww_n_lower_bnds,
                                 None,
                                 10.0,
-                                e_bnds)
+                                e_bnds,
+                                max_lower_bound_ratio=1.5)
 
     ww_p = openmc.WeightWindows(ww_mesh,
                                 ww_p_lower_bnds,
                                 None,
                                 10.0,
-                                e_bnds)
+                                e_bnds,
+                                max_lower_bound_ratio=1.5)
 
     model.settings.weight_windows = [ww_n, ww_p]