Add option for survival biasing source normalization (#3070)

Co-authored-by: Patrick Shriwise <[email protected]>
Co-authored-by: Paul Romano <[email protected]>

Author: 3 people

docs/source/io_formats/settings.rst

@@ -81,6 +81,13 @@ time.
 
     *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

@@ -434,6 +434,7 @@ class ParticleData : public GeometryState {
   int g_last_;
 
   double wgt_ {1.0};
+  double wgt_born_ {1.0};
   double mu_;
   double time_ {0.0};
   double time_last_ {0.0};
@@ -532,12 +533,20 @@ class ParticleData : public GeometryState {
   int& g_last() { return g_last_; }
   const int& g_last() const { return g_last_; }
 
-  // Statistic weight of particle. Setting to zero
-  // indicates that the particle is dead.
+  // Statistic weight of particle. Setting to zero indicates that the particle
+  // is dead.
   double& wgt() { return wgt_; }
   double wgt() const { return wgt_; }
+
+  // Statistic weight of particle at birth
+  double& wgt_born() { return wgt_born_; }
+  double wgt_born() const { return wgt_born_; }
+
+  // Statistic weight of particle at last collision
   double& wgt_last() { return wgt_last_; }
   const double& wgt_last() const { return wgt_last_; }
+
+  // Whether particle is alive
   bool alive() const { return wgt_ != 0.0; }
 
   // Polar scattering angle after a collision

include/openmc/settings.h

@@ -61,6 +61,7 @@ extern bool surf_source_write;       //!< write surface source file?
 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?

openmc/settings.py

@@ -50,15 +50,19 @@ class Settings:
         Indicate whether fission neutrons should be created or not.
     cutoff : dict
         Dictionary defining weight cutoff, energy cutoff and time cutoff. The
-        dictionary may have ten keys, 'weight', 'weight_avg', 'energy_neutron',
-        'energy_photon', 'energy_electron', 'energy_positron', 'time_neutron',
-        'time_photon', 'time_electron', and 'time_positron'. Value for 'weight'
-        should be a float indicating weight cutoff below which particle undergo
-        Russian roulette. Value for 'weight_avg' should be a float indicating
-        weight assigned to particles that are not killed after Russian roulette.
-        Value of energy should be a float indicating energy in eV below which
-        particle type will be killed. Value of time should be a float in
-        seconds. Particles will be killed exactly at the specified time.
+        dictionary may have the following keys, 'weight', 'weight_avg',
+        'survival_normalization', 'energy_neutron', 'energy_photon',
+        'energy_electron', 'energy_positron', 'time_neutron', 'time_photon',
+        'time_electron', and 'time_positron'. Value for 'weight' should be a
+        float indicating weight cutoff below which particle undergo Russian
+        roulette. Value for 'weight_avg' should be a float indicating weight
+        assigned to particles that are not killed after Russian roulette. Value
+        of energy should be a float indicating energy in eV below which particle
+        type will be killed. Value of time should be a float in seconds.
+        Particles will be killed exactly at the specified time. Value for
+        'survival_normalization' is a bool indicating whether or not the weight
+        cutoff parameters will be applied relative to the particle's starting
+        weight or to its current weight.
     delayed_photon_scaling : bool
         Indicate whether to scale the fission photon yield by (EGP + EGD)/EGP
         where EGP is the energy release of prompt photons and EGD is the energy
@@ -162,20 +166,20 @@ class Settings:
             'naive', 'simulation_averaged', or 'hybrid'.
             The default is 'hybrid'.
         :source_shape:
-            Assumed shape of the source distribution within each source
-            region. Options are 'flat' (default), 'linear', or 'linear_xy'.
+            Assumed shape of the source distribution within each source region.
+            Options are 'flat' (default), 'linear', or 'linear_xy'.
         :volume_normalized_flux_tallies:
-            Whether to normalize flux tallies by volume (bool). The default
-            is 'False'. When enabled, flux tallies will be reported in units of
-            cm/cm^3. When disabled, flux tallies will be reported in units
-            of cm (i.e., total distance traveled by neutrons in the spatial
-            tally region).
+            Whether to normalize flux tallies by volume (bool). The default is
+            'False'. When enabled, flux tallies will be reported in units of
+            cm/cm^3. When disabled, flux tallies will be reported in units of cm
+            (i.e., total distance traveled by neutrons in the spatial tally
+            region).
         :adjoint:
             Whether to run the random ray solver in adjoint mode (bool). The
             default is 'False'.
         :sample_method:
-            Sampling method for the ray starting location and direction of travel.
-            Options are `prng` (default) or 'halton`.
+            Sampling method for the ray starting location and direction of
+            travel. Options are `prng` (default) or 'halton`.
 
         .. versionadded:: 0.15.0
     resonance_scattering : dict
@@ -899,6 +903,8 @@ def cutoff(self, cutoff: dict):
                 cv.check_type('average survival weight', cutoff[key], Real)
                 cv.check_greater_than('average survival weight',
                                       cutoff[key], 0.0)
+            elif key == 'survival_normalization':
+                cv.check_type('survival normalization', cutoff[key], bool)
             elif key in ['energy_neutron', 'energy_photon', 'energy_electron',
                          'energy_positron']:
                 cv.check_type('energy cutoff', cutoff[key], Real)
@@ -1364,7 +1370,8 @@ def _create_cutoff_subelement(self, root):
             element = ET.SubElement(root, "cutoff")
             for key, value in self._cutoff.items():
                 subelement = ET.SubElement(element, key)
-                subelement.text = str(value)
+                subelement.text = str(value) if key != 'survival_normalization' \
+                    else str(value).lower()
 
     def _create_entropy_mesh_subelement(self, root, mesh_memo=None):
         if self.entropy_mesh is None:
@@ -1797,10 +1804,14 @@ def _cutoff_from_xml_element(self, root):
             self.cutoff = {}
             for key in ('energy_neutron', 'energy_photon', 'energy_electron',
                         'energy_positron', 'weight', 'weight_avg', 'time_neutron',
-                        'time_photon', 'time_electron', 'time_positron'):
+                        'time_photon', 'time_electron', 'time_positron',
+                        'survival_normalization'):
                 value = get_text(elem, key)
                 if value is not None:
-                    self.cutoff[key] = float(value)
+                    if key == 'survival_normalization':
+                        self.cutoff[key] = value in ('true', '1')
+                    else:
+                        self.cutoff[key] = float(value)
 
     def _entropy_mesh_from_xml_element(self, root, meshes):
         text = get_text(root, 'entropy_mesh')

src/physics.cpp

@@ -153,7 +153,13 @@ void sample_neutron_reaction(Particle& p)
 
   // Play russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
-    if (p.wgt() < settings::weight_cutoff) {
+    // if survival normalization is on, use normalized weight cutoff and
+    // normalized weight survive
+    if (settings::survival_normalization) {
+      if (p.wgt() < settings::weight_cutoff * p.wgt_born()) {
+        russian_roulette(p, settings::weight_survive * p.wgt_born());
+      }
+    } else if (p.wgt() < settings::weight_cutoff) {
       russian_roulette(p, settings::weight_survive);
     }
   }

src/physics_mg.cpp

@@ -68,7 +68,13 @@ void sample_reaction(Particle& p)
 
   // Play Russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
-    if (p.wgt() < settings::weight_cutoff) {
+    // if survival normalization is applicable, use normalized weight cutoff and
+    // normalized weight survive
+    if (settings::survival_normalization) {
+      if (p.wgt() < settings::weight_cutoff * p.wgt_born()) {
+        russian_roulette(p, settings::weight_survive * p.wgt_born());
+      }
+    } else if (p.wgt() < settings::weight_cutoff) {
       russian_roulette(p, settings::weight_survive);
     }
   }

src/settings.cpp

@@ -69,6 +69,7 @@ bool surf_source_write {false};
 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};
@@ -624,6 +625,10 @@ void read_settings_xml(pugi::xml_node root)
     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

@@ -586,6 +586,9 @@ void initialize_history(Particle& p, int64_t index_source)
   // Reset weight window ratio
   p.ww_factor() = 0.0;
 
+  // set particle history start weight
+  p.wgt_born() = p.wgt();
+
   // Reset pulse_height_storage
   std::fill(p.pht_storage().begin(), p.pht_storage().end(), 0);
 
@@ -610,7 +613,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();
 

tests/unit_tests/test_settings.py

@@ -26,6 +26,7 @@ def test_export_to_xml(run_in_tmpdir):
     s.plot_seed = 100
     s.survival_biasing = True
     s.cutoff = {'weight': 0.25, 'weight_avg': 0.5, 'energy_neutron': 1.0e-5,
+                'survival_normalization': True,
                 'energy_photon': 1000.0, 'energy_electron': 1.0e-5,
                 'energy_positron': 1.0e-5, 'time_neutron': 1.0e-5,
                 'time_photon': 1.0e-5, 'time_electron': 1.0e-5,
@@ -99,6 +100,7 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.seed == 17
     assert s.survival_biasing
     assert s.cutoff == {'weight': 0.25, 'weight_avg': 0.5,
+                        'survival_normalization': True,
                         'energy_neutron': 1.0e-5, 'energy_photon': 1000.0,
                         'energy_electron': 1.0e-5, 'energy_positron': 1.0e-5,
                         'time_neutron': 1.0e-5, 'time_photon': 1.0e-5,