Weight Window Birth Scaling (#3459)

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

include/openmc/particle_data.h

@@ -435,6 +435,7 @@ private:
 
   double wgt_ {1.0};
   double wgt_born_ {1.0};
+  double wgt_ww_born_ {-1.0};
   double mu_;
   double time_ {0.0};
   double time_last_ {0.0};
@@ -545,6 +546,10 @@ public:
   double& wgt_born() { return wgt_born_; }
   double wgt_born() const { return wgt_born_; }
 
+  // Weight window value at birth
+  double& wgt_ww_born() { return wgt_ww_born_; }
+  const double& wgt_ww_born() const { return wgt_ww_born_; }
+
   // Statistic weight of particle at last collision
   double& wgt_last() { return wgt_last_; }
   const double& wgt_last() const { return wgt_last_; }

src/simulation.cpp

@@ -614,6 +614,10 @@ void initialize_history(Particle& p, int64_t index_source)
   // Set particle track.
   p.write_track() = check_track_criteria(p);
 
+  // Set the particle's initial weight window value.
+  p.wgt_ww_born() = -1.0;
+  apply_weight_windows(p);
+
   // Display message if high verbosity or trace is on
   if (settings::verbosity >= 9 || p.trace()) {
     write_message("Simulating Particle {}", p.id());

src/weight_windows.cpp

@@ -73,10 +73,26 @@ void apply_weight_windows(Particle& p)
     if (weight_window.is_valid())
       break;
   }
+
+  // If particle has not yet had its birth weight window value set, set it to
+  // the current weight window (or 1.0 if not born in a weight window).
+  if (p.wgt_ww_born() == -1.0) {
+    if (weight_window.is_valid()) {
+      p.wgt_ww_born() =
+        (weight_window.lower_weight + weight_window.upper_weight) / 2;
+    } else {
+      p.wgt_ww_born() = 1.0;
+    }
+  }
+
   // particle is not in any of the ww domains, do nothing
   if (!weight_window.is_valid())
     return;
 
+  // Normalize weight windows based on particle's starting weight
+  // and the value of the weight window the particle was born in.
+  weight_window.scale(p.wgt_born() / p.wgt_ww_born());
+
   // get the paramters
   double weight = p.wgt();