Adding otf ww adjustment using lower bound ratio.

include/openmc/particle_data.h

@@ -309,7 +309,9 @@ 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 +437,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,10 @@ 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 weight window's lower bound to particle weight.
+        A factor will be applied to raise the weight window to the particle's weight
+        during transport if exceeded. Only applied if > 1.0. (Default: 1.0).
     max_split : int
         Maximum allowable number of particles when splitting
     weight_cutoff : float
@@ -76,6 +80,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 weight window's lower bound to particle
+        weight. (Default: 1.0)
     max_split : int
         Maximum allowable number of particles when splitting
     weight_cutoff : float
@@ -91,7 +98,7 @@ 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_lb_ratio=1.0, max_split=10, weight_cutoff=1.e-38, id=None):
         self.mesh = mesh
         self.id = id
         self.particle_type = particle_type
@@ -121,6 +128,7 @@ class WeightWindows(IDManagerMixin):
                              'do not match')
 
         self.survival_ratio = survival_ratio
+        self.max_lower_bound_ratio = max_lb_ratio
         self.max_split = max_split
         self.weight_cutoff = weight_cutoff
 
@@ -192,6 +200,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, True)
+        self._max_lower_bound_ratio = val
+
     @property
     def max_split(self):
         return self._max_split
@@ -241,6 +259,9 @@ class WeightWindows(IDManagerMixin):
         subelement = ET.SubElement(element, 'survival_ratio')
         subelement.text = str(self.survival_ratio)
 
+        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 +299,7 @@ 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_lb_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 +311,7 @@ class WeightWindows(IDManagerMixin):
             energy_bins=ebins,
             particle_type=particle_type,
             survival_ratio=survival_ratio,
+            max_lb_ratio=max_lb_ratio,
             max_split=max_split,
             weight_cutoff=weight_cutoff,
             id=id
@@ -318,6 +341,7 @@ class WeightWindows(IDManagerMixin):
         lower_ww_bounds = group['lower_ww_bounds'][()]
         upper_ww_bounds = group['upper_ww_bounds'][()]
         survival_ratio = group['survival_ratio'][()]
+        max_lb_ratio = group['max_lower_bound_ratio'][()]
         max_split = group['max_split'][()]
         weight_cutoff = group['weight_cutoff'][()]
 
@@ -328,6 +352,7 @@ class WeightWindows(IDManagerMixin):
             energy_bins=ebins,
             particle_type=ptype,
             survival_ratio=survival_ratio,
+            max_lb_ratio=max_lb_ratio,
             max_split=max_split,
             weight_cutoff=weight_cutoff,
             id=id

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"));
@@ -254,6 +275,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/test.py

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