Merge pull request #18 from yrrepy/Surv-Norm_PR-review

Commits to address review

docs/source/methods/neutron_physics.rst

@@ -1688,12 +1688,12 @@ This parameter can be used to make the adjustable parameters of survival biasing
 
 It is currently only implemented for phase-space sources (MCPL file sources
 and HDF5 surface sources). For these, the normalization of the parameters
-is done per source particle. That is for each history the parameters,
+is done per source particle. That is, for each history, the parameters
 :math:`w_c` and :math:`w_s`, are multiplied by the start weight of the current
 history.
 
 This normalization is recommended for problems where numerous source particles are
-initialized below the survival biasing weight cutoff, :math:`w_c`; to prevent them
+initialized below the survival biasing weight cutoff :math:`w_c`; to prevent them
 from being immediately rouletted. This can notably be the case for biased sources.
 
 Weight Windows

include/openmc/particle_data.h

@@ -264,8 +264,7 @@ class ParticleData {
 
   // Other physical data
   double wgt_ {1.0};       //!< particle weight
-  double wgt_cutoff_;      //!< particle weight cutoff
-  double wgt_survive_;     //!< particle weight survive
+  double wgt0_ {1.0};      //!< particle start weight
   double mu_;              //!< angle of scatter
   double time_ {0.0};      //!< time in [s]
   double time_last_ {0.0}; //!< previous time in [s]
@@ -399,12 +398,9 @@ class ParticleData {
 
   double& wgt() { return wgt_; }
   double wgt() const { return wgt_; }
-  double& wgt_cutoff() { return wgt_cutoff_; }
-  double wgt_cutoff() const { return wgt_cutoff_; }
-  void wgt_cutoff(double cutoff) { wgt_cutoff_ = cutoff; }
-  double& wgt_survive() { return wgt_survive_; }
-  double wgt_survive() const { return wgt_survive_; }
-  void wgt_survive(double survive) { wgt_survive_ = survive; }
+  double& wgt0() { return wgt0_; }
+  double wgt0() const { return wgt0_; }
+
   double& mu() { return mu_; }
   const double& mu() const { return mu_; }
   double& time() { return time_; }

src/physics.cpp

@@ -157,13 +157,13 @@ void sample_neutron_reaction(Particle& p)
 
   // Play russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
-    double weight_cutoff = settings::weight_cutoff;
-    // if survival normalization is applicable, use normalized weight cutoff
-    if((settings::source_file || settings::surf_source_read)&&(settings::survival_normalization)){
-      weight_cutoff  = p.wgt_cutoff();
-    }
-    if (p.wgt() < weight_cutoff) {
-      russian_roulette(p, settings::weight_survive);
+    // if survival normalization is applicable, use normalized weight cutoff and normalized weight survive
+    if ((settings::source_file || settings::surf_source_read)&&(settings::survival_normalization)) {
+      if (p.wgt() < settings::weight_cutoff*p.wgt0()) {
+        russian_roulette(p, settings::weight_survive*p.wgt0());
+      } 
+    } else if (p.wgt() < settings::weight_cutoff) {
+            russian_roulette(p, settings::weight_survive);
     }
   }
 }

src/physics_common.cpp

@@ -11,15 +11,11 @@ namespace openmc {
 
 void russian_roulette(Particle& p, double weight_survive)
 {
-  // if survival normalization is applicable, use normalized weight survive
-  if((settings::source_file || settings::surf_source_read)&&(settings::survival_normalization && settings::survival_biasing)){
-    weight_survive = p.wgt_survive();
-  }
   if (weight_survive * prn(p.current_seed()) < p.wgt()) {
     p.wgt() = weight_survive;
   } else {
     p.wgt() = 0.;
-  } 
+  }
 }
 
 } // namespace openmc

src/physics_mg.cpp

@@ -64,12 +64,7 @@ void sample_reaction(Particle& p)
 
   // Play Russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
-    double weight_cutoff = settings::weight_cutoff;
-    // if survival normalization is applicable, use normalized weight cutoff
-    if((settings::source_file || settings::surf_source_read)&&(settings::survival_normalization)){
-      weight_cutoff  = p.wgt_cutoff();
-    }
-    if (p.wgt() < weight_cutoff) {
+    if (p.wgt() < settings::weight_cutoff) {
       russian_roulette(p, settings::weight_survive);
     }
   }

src/simulation.cpp

@@ -533,19 +533,15 @@ void initialize_history(Particle& p, int64_t index_source)
     auto site = sample_external_source(&seed);
     p.from_source(&site);
   }
-  // normalize biasing weight parameters; multiply them by start weight of initialized history
-  // applicable only for MCPL and HDF5 phase-space sources.
-  if(settings::source_file || settings::surf_source_read){
-    if(settings::survival_normalization && settings::survival_biasing){
-      p.wgt_cutoff(settings::weight_cutoff * p.wgt());
-      p.wgt_survive(settings::weight_survive * p.wgt());   
-    }
-  }
+
   p.current_work() = index_source;
 
   // set identifier for particle
   p.id() = simulation::work_index[mpi::rank] + index_source;
 
+  // set particle history start weight
+  p.wgt0(p.wgt());
+
   // set progeny count to zero
   p.n_progeny() = 0;
 
@@ -582,7 +578,7 @@ void initialize_history(Particle& p, int64_t index_source)
     write_message("Simulating Particle {}", p.id());
   }
 
-// Add paricle's starting weight to count for normalizing tallies later
+// Add particle's starting weight to count for normalizing tallies later
 #pragma omp atomic
   simulation::total_weight += p.wgt();