Survival Normalization PR

docs/source/io_formats/settings.rst

@@ -79,6 +79,12 @@ you care. This element has the following attributes/sub-elements:
 
     *Default*: 1.0
 
+  :survival_normalization:
+    If this element is set to "true", this will enable the use of survival normalization, whereby the weight parameters, weight and weight_avg, are multiplied per 
+    history by the start weight of said history. Currently only implemented for phase-space sources.
+
+    *Default*: true
+
   :energy_neutron:
     The energy under which neutrons will be killed.
 

docs/source/methods/neutron_physics.rst

@@ -1678,6 +1678,23 @@ default, the cutoff weight in OpenMC is :math:`w_c = 0.25` and the survival
 weight is :math:`w_s = 1.0`. These parameters vary from one Monte Carlo code to
 another.
 
+.. _survival_normalization:
+
+Survival Normalization
+----------------
+
+This parameter can be used to make the adjustable parameters of survival biasing,
+:math:`w_c` and :math:`w_s`, relative to the starting weight of a source particle.
+
+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,
+:math:`w_c` and :math:`w_s`, are multiplied by the start weight of the current
+history.
+
+For most problems this normalization is found to improve Figure of Merits and is 
+recommended.
+
 Weight Windows
 --------------
 

include/openmc/particle_data.h

@@ -264,6 +264,8 @@ class ParticleData {
 
   // Other physical data
   double wgt_ {1.0};       //!< particle weight
+  double wgt_cutoff_;      //!< particle weight cutoff
+  double wgt_survive_;     //!< particle weight survive
   double mu_;              //!< angle of scatter
   double time_ {0.0};      //!< time in [s]
   double time_last_ {0.0}; //!< previous time in [s]
@@ -397,6 +399,12 @@ 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& mu() { return mu_; }
   const double& mu() const { return mu_; }
   double& time() { return time_; }

include/openmc/settings.h

@@ -1,4 +1,4 @@
 #ifndef OPENMC_SETTINGS_H
 #define OPENMC_SETTINGS_H
 
 //! \file settings.h
@@ -61,6 +61,8 @@
 extern "C" bool weight_windows_on; //!< are weight windows are enabled?
 extern bool write_all_tracks;      //!< write track files for every particle?
 extern bool write_initial_source;  //!< write out initial source file?
+extern bool survival_normalization;//!< using survival normalization?
+extern bool source_file;           //!< is there a source file?
 
 // Paths to various files
 extern std::string path_cross_sections; //!< path to cross_sections.xml

src/physics.cpp

@@ -1,4 +1,4 @@
 #include "openmc/physics.h"
 
 #include "openmc/bank.h"
 #include "openmc/bremsstrahlung.h"
@@ -157,7 +157,12 @@
 
   // Play russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
-    if (p.wgt() < settings::weight_cutoff) {
+    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);
     }
   }

src/physics_common.cpp

@@ -1,4 +1,4 @@
 #include "openmc/physics_common.h"
 
 #include "openmc/random_lcg.h"
 #include "openmc/settings.h"
@@ -11,11 +11,15 @@
 
 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

@@ -1,4 +1,4 @@
 #include "openmc/physics_mg.h"
 
 #include <stdexcept>
 
@@ -64,7 +64,12 @@
 
   // Play Russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
-    if (p.wgt() < settings::weight_cutoff) {
+    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);
     }
   }

src/settings.cpp

@@ -1,4 +1,4 @@
 #include "openmc/settings.h"
 
 #include <cmath>  // for ceil, pow
 #include <limits> // for numeric_limits
@@ -75,6 +75,8 @@
 bool weight_windows_on {false};
 bool write_all_tracks {false};
 bool write_initial_source {false};
+bool survival_normalization {true};
+bool source_file {false};
 
 std::string path_cross_sections;
 std::string path_input;
@@ -442,6 +444,7 @@
   for (pugi::xml_node node : root.children("source")) {
     if (check_for_node(node, "file")) {
       auto path = get_node_value(node, "file", false, true);
+      source_file = true;
       if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) {
         auto sites = mcpl_source_sites(path);
         model::external_sources.push_back(make_unique<FileSource>(sites));
@@ -513,6 +516,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

@@ -1,4 +1,4 @@
 #include "openmc/simulation.h"
 
 #include "openmc/bank.h"
 #include "openmc/capi.h"
@@ -533,6 +533,14 @@
     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