Statistical weights in IndependentSource

docs/source/usersguide/settings.rst

@@ -183,6 +183,7 @@ source distributions and has four main attributes that one can set:
 :attr:`IndependentSource.energy`, which defines the energy distribution, and
 :attr:`IndependentSource.time`, which defines the time distribution.
 
+
 The spatial distribution can be set equal to a sub-class of
 :class:`openmc.stats.Spatial`; common choices are :class:`openmc.stats.Point` or
 :class:`openmc.stats.Box`. To independently specify distributions in the
@@ -223,6 +224,7 @@ distribution. This could be a probability mass function
 (:class:`openmc.stats.Tabular`). By default, if no time distribution is
 specified, particles are started at :math:`t=0`.
 
+
 As an example, to create an isotropic, 10 MeV monoenergetic source uniformly
 distributed over a cube centered at the origin with an edge length of 10 cm, and
 emitting a pulse of particles from 0 to 10 µs, one
@@ -250,6 +252,29 @@ sampled 70% of the time and another that should be sampled 30% of the time::
 
   settings.source = [src1, src2]
 
+When the relative strengths are several orders of magnitude different, it may 
+happen that not enough statistics are obtained from the lower strength source. 
+This can be improved by sampling particles with different weights with 
+probability 1, instead of sampling particles with weight 1 with different 
+probability. The :attr:`Settings.strengths_as_weights` attribute can be used 
+to enable this option::
+
+  settings.strengths_as_weights = True
+
+Then, the strength values will be used as statistical weights:: 
+
+
+  src1 = openmc.IndependentSource()
+  src1.strength = 100.0
+  ...
+
+  src2 = openmc.IndependentSource()
+  src2.strength = 1.0
+  ...
+
+  settings.strengths_as_weights = True
+  settings.source = [src1, src2]
+
 Finally, the :attr:`IndependentSource.particle` attribute can be used to
 indicate the source should be composed of particles other than neutrons. For
 example, the following would generate a photon source::

include/openmc/settings.h

@@ -57,6 +57,7 @@ extern bool source_latest;         //!< write latest source at each batch?
 extern bool source_separate;       //!< write source to separate file?
 extern bool source_write;          //!< write source in HDF5 files?
 extern bool source_mcpl_write;     //!< write source in mcpl files?
+extern bool strengths_as_weights;  //!< convert strength to statistical weight?
 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?

openmc/settings.py

@@ -202,6 +202,8 @@ class Settings:
         Options for writing state points. Acceptable keys are:
 
         :batches: list of batches at which to write statepoint files
+    strengths_as_weights : bool
+        Whether to convert strength to statistical weight.
     surf_source_read : dict
         Options for reading surface source points. Acceptable keys are:
 
@@ -325,6 +327,7 @@ def __init__(self, **kwargs):
         self._photon_transport = None
         self._plot_seed = None
         self._ptables = None
+        self._strengths_as_weights = None
         self._seed = None
         self._survival_biasing = None
 
@@ -572,6 +575,15 @@ def photon_transport(self, photon_transport: bool):
         cv.check_type('photon transport', photon_transport, bool)
         self._photon_transport = photon_transport
 
+    @property
+    def strengths_as_weights(self) -> bool:
+        return self._strengths_as_weights
+
+    @strengths_as_weights.setter
+    def strengths_as_weights(self, strengths_as_weights: bool):
+        cv.check_type('strength as weights', strengths_as_weights, bool)
+        self._strengths_as_weights = strengths_as_weights
+
     @property
     def plot_seed(self):
         return self._plot_seed
@@ -1216,6 +1228,11 @@ def _create_statepoint_subelement(self, root):
                 subelement.text = ' '.join(
                     str(x) for x in self._statepoint['batches'])
 
+    def _create_strengths_as_weights_subelement(self, root):
+        if self._strengths_as_weights is not None:
+            element = ET.SubElement(root, "strengths_as_weights")
+            element.text = str(self._strengths_as_weights).lower()
+
     def _create_sourcepoint_subelement(self, root):
         if self._sourcepoint:
             element = ET.SubElement(root, "source_point")
@@ -1697,6 +1714,11 @@ def _photon_transport_from_xml_element(self, root):
         if text is not None:
             self.photon_transport = text in ('true', '1')
 
+    def _strengths_as_weights_from_xml_element(self, root):
+        text = get_text(root, 'strengths_as_weights')
+        if text is not None:
+            self.strengths_as_weights = text in ('true', '1')
+
     def _plot_seed_from_xml_element(self, root):
         text = get_text(root, 'plot_seed')
         if text is not None:
@@ -1948,6 +1970,7 @@ def to_xml_element(self, mesh_memo=None):
         self._create_energy_mode_subelement(element)
         self._create_max_order_subelement(element)
         self._create_photon_transport_subelement(element)
+        self._create_strengths_as_weights_subelement(element)
         self._create_plot_seed_subelement(element)
         self._create_ptables_subelement(element)
         self._create_seed_subelement(element)
@@ -2054,6 +2077,7 @@ def from_xml_element(cls, elem, meshes=None):
         settings._energy_mode_from_xml_element(elem)
         settings._max_order_from_xml_element(elem)
         settings._photon_transport_from_xml_element(elem)
+        settings._strengths_as_weights_from_xml_element(elem)
         settings._plot_seed_from_xml_element(elem)
         settings._ptables_from_xml_element(elem)
         settings._seed_from_xml_element(elem)

src/settings.cpp

@@ -65,6 +65,7 @@ bool source_latest {false};
 bool source_separate {false};
 bool source_write {true};
 bool source_mcpl_write {false};
+bool strengths_as_weights {false};
 bool surf_source_write {false};
 bool surf_mcpl_write {false};
 bool surf_source_read {false};
@@ -782,6 +783,11 @@ void read_settings_xml(pugi::xml_node root)
     sourcepoint_batch = statepoint_batch;
   }
 
+  // Check is the user specified to convert strength to statistical weight
+  if (check_for_node(root, "strengths_as_weights")) {
+    strengths_as_weights = get_node_value_bool(root, "strengths_as_weights");
+  }
+
   // Check if the user has specified to write surface source
   if (check_for_node(root, "surf_source_write")) {
     surf_source_write = true;

src/source.cpp

@@ -394,6 +394,10 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
 
     // Sample particle creation time
     site.time = time_->sample(seed);
+    // Set particle creation weight
+    if (openmc::settings::strengths_as_weights) {
+      site.wgt = strength();
+    }
   }
 
   // Increment number of accepted samples
@@ -611,15 +615,23 @@ SourceSite sample_external_source(uint64_t* seed)
   // Determine total source strength
   double total_strength = 0.0;
   for (auto& s : model::external_sources)
-    total_strength += s->strength();
+    if (openmc::settings::strengths_as_weights) {
+      total_strength += 1.0;
+    } else {
+      total_strength += s->strength();
+    }
 
   // Sample from among multiple source distributions
   int i = 0;
   if (model::external_sources.size() > 1) {
     double xi = prn(seed) * total_strength;
     double c = 0.0;
     for (; i < model::external_sources.size(); ++i) {
-      c += model::external_sources[i]->strength();
+      if (openmc::settings::strengths_as_weights) {
+        c += 1.0;
+      } else {
+        c += model::external_sources[i]->strength();
+      }
       if (xi < c)
         break;
     }

tests/unit_tests/test_source_weight.py

@@ -0,0 +1,38 @@
+import openmc
+import openmc.stats
+import h5py
+
+
+def test_musurface(run_in_tmpdir):
+    sphere = openmc.Sphere(r=1.0, boundary_type='vacuum')
+    cell = openmc.Cell(region=-sphere, fill=None)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell])
+    model.settings.particles = 100
+    model.settings.batches = 1
+    E = 1.0
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point(),
+        angle=openmc.stats.Isotropic(),
+        energy=openmc.stats.delta_function(E),
+        strength=100.0
+    )
+    model.settings.strengths_as_weights = True
+    model.settings.run_mode = "fixed source"
+    model.settings.surf_source_write = {
+        "max_particles": 100,
+    }
+
+    # Run OpenMC
+    sp_filename = model.run()
+
+    # All contributions should show up in last bin
+    with h5py.File("surface_source.h5", "r") as f:
+        source = f["source_bank"]
+
+        assert len(source) == 100
+
+        for point in source:
+            assert point["wgt"] == 100.0
+
+