Survival Normalization 0.15.1

docs/source/io_formats/settings.rst

@@ -80,6 +80,12 @@ time.
     roulette.
 
     *Default*: 1.0
+  :survival_normalization:
+    If this element is set to "true", this will enable the use of survival biasing source normalization, whereby the weight parameters, weight and weight_avg, are multiplied per 
+    history by the start weight of said history. 
+
+    *Default*: false
+
 
   :energy_neutron:
     The energy under which neutrons will be killed.

include/openmc/particle_data.h

@@ -406,6 +406,7 @@ class ParticleData : public GeometryState {
   int g_last_;
 
   double wgt_ {1.0};
+  double wgt0_ {1.0};
   double mu_;
   double time_ {0.0};
   double time_last_ {0.0};
@@ -507,6 +508,8 @@ class ParticleData : public GeometryState {
   // indicates that the particle is dead.
   double& wgt() { return wgt_; }
   double wgt() const { return wgt_; }
+  double& wgt0() { return wgt0_; }
+  double wgt0() const { return wgt0_; }
   double& wgt_last() { return wgt_last_; }
   const double& wgt_last() const { return wgt_last_; }
   bool alive() const { return wgt_ != 0.0; }

include/openmc/settings.h

@@ -1,4 +1,4 @@
 #ifndef OPENMC_SETTINGS_H
 #define OPENMC_SETTINGS_H
 
 //! \file settings.h
@@ -61,6 +61,7 @@
 extern bool surf_mcpl_write;       //!< write surface mcpl file?
 extern bool surf_source_read;      //!< read surface source file?
 extern bool survival_biasing;      //!< use survival biasing?
+extern bool survival_normalization;//!< use survival normalization?
 extern bool temperature_multipole; //!< use multipole data?
 extern "C" bool trigger_on;        //!< tally triggers enabled?
 extern bool trigger_predict;       //!< predict batches for triggers?

src/physics.cpp

@@ -1,4 +1,4 @@
 #include "openmc/physics.h"
 
 #include "openmc/bank.h"
 #include "openmc/bremsstrahlung.h"
@@ -152,8 +152,13 @@
 
   // Play russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
-    if (p.wgt() < settings::weight_cutoff) {
-      russian_roulette(p, settings::weight_survive);
+    // if survival normalization is on, use normalized weight cutoff and normalized weight survive
+    if (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_mg.cpp

@@ -1,4 +1,4 @@
 #include "openmc/physics_mg.h"
 
 #include <stdexcept>
 
@@ -64,8 +64,13 @@
 
   // Play Russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
-    if (p.wgt() < settings::weight_cutoff) {
-      russian_roulette(p, settings::weight_survive);
+    // if survival normalization is applicable, use normalized weight cutoff and normalized weight survive
+    if (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/settings.cpp

@@ -1,4 +1,4 @@
 #include "openmc/settings.h"
 
 #include <cmath>  // for ceil, pow
 #include <limits> // for numeric_limits
@@ -69,6 +69,7 @@
 bool surf_mcpl_write {false};
 bool surf_source_read {false};
 bool survival_biasing {false};
+bool survival_normalization {false};
 bool temperature_multipole {false};
 bool trigger_on {false};
 bool trigger_predict {false};
@@ -566,6 +567,9 @@
     if (check_for_node(node_cutoff, "weight_avg")) {
       weight_survive = std::stod(get_node_value(node_cutoff, "weight_avg"));
     }
+    if (check_for_node(node_cutoff, "survival_normalization")){ 
+      survival_normalization = get_node_value_bool(node_cutoff, "survival_normalization");
+    }
     if (check_for_node(node_cutoff, "energy_neutron")) {
       energy_cutoff[0] =
         std::stod(get_node_value(node_cutoff, "energy_neutron"));

src/simulation.cpp

@@ -567,6 +567,9 @@ void initialize_history(Particle& p, int64_t index_source)
   // set identifier for particle
   p.id() = simulation::work_index[mpi::rank] + index_source;
 
+  // set particle history start weight
+  p.wgt0() = p.wgt(); // p.wgt0(p.wgt());
+
   // set progeny count to zero
   p.n_progeny() = 0;
 
@@ -603,7 +606,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();