Statistical weights in IndependentSource

include/openmc/source.h

@@ -92,13 +92,16 @@ class Source {
   bool satisfies_spatial_constraints(Position r) const;
   bool satisfies_energy_constraints(double E) const;
   bool satisfies_time_constraints(double time) const;
+  bool satisfies_weight_constraints(double weight) const;
 
   // Data members
   double strength_ {1.0};                  //!< Source strength
   std::unordered_set<int32_t> domain_ids_; //!< Domains to reject from
   DomainType domain_type_;                 //!< Domain type for rejection
   std::pair<double, double> time_bounds_ {-std::numeric_limits<double>::max(),
     std::numeric_limits<double>::max()}; //!< time limits
+  std::pair<double, double> weight_bounds_ {
+    0, std::numeric_limits<double>::max()}; //!< weight limits
   std::pair<double, double> energy_bounds_ {
     0, std::numeric_limits<double>::max()}; //!< energy limits
   bool only_fissionable_ {
@@ -115,8 +118,8 @@ class Source {
 class IndependentSource : public Source {
 public:
   // Constructors
-  IndependentSource(
-    UPtrSpace space, UPtrAngle angle, UPtrDist energy, UPtrDist time);
+  IndependentSource(UPtrSpace space, UPtrAngle angle, UPtrDist energy,
+    UPtrDist time, double weight);
   explicit IndependentSource(pugi::xml_node node);
 
   //! Sample from the external source distribution
@@ -132,6 +135,7 @@ class IndependentSource : public Source {
   UnitSphereDistribution* angle() const { return angle_.get(); }
   Distribution* energy() const { return energy_.get(); }
   Distribution* time() const { return time_.get(); }
+  const double weight() const { return weight_; }
 
   // Make domain type and ids available
   DomainType domain_type() const { return domain_type_; }
@@ -148,6 +152,7 @@ class IndependentSource : public Source {
   UPtrAngle angle_;                               //!< Angular distribution
   UPtrDist energy_;                               //!< Energy distribution
   UPtrDist time_;                                 //!< Time distribution
+  double weight_;                                 //!< weight value
 };
 
 //==============================================================================

openmc/source.py

@@ -258,6 +258,8 @@ class IndependentSource(SourceBase):
         Energy distribution of source sites
     time : openmc.stats.Univariate
         time distribution of source sites
+    weight : float
+        weight value of source sites
     strength : float
         Strength of the source
     particle : {'neutron', 'photon'}
@@ -292,6 +294,8 @@ class IndependentSource(SourceBase):
         Energy distribution of source sites
     time : openmc.stats.Univariate or None
         time distribution of source sites
+    weight : float or None
+        weight value of source sites
     strength : float
         Strength of the source
     type : str
@@ -314,6 +318,7 @@ def __init__(
         angle: openmc.stats.UnitSphere | None = None,
         energy: openmc.stats.Univariate | None = None,
         time: openmc.stats.Univariate | None = None,
+        weight: float | None = None,
         strength: float = 1.0,
         particle: str = 'neutron',
         domains: Sequence[openmc.Cell | openmc.Material | openmc.Universe] | None = None,
@@ -330,6 +335,7 @@ def __init__(
         self._angle = None
         self._energy = None
         self._time = None
+        self._weight = None
 
         if space is not None:
             self.space = space
@@ -339,6 +345,8 @@ def __init__(
             self.energy = energy
         if time is not None:
             self.time = time
+        if weight is not None:
+            self.weight = weight
         self.particle = particle
 
     @property
@@ -398,6 +406,15 @@ def time(self, time):
         cv.check_type('time distribution', time, Univariate)
         self._time = time
 
+    @property
+    def weight(self):
+        return self._weight
+
+    @weight.setter
+    def weight(self, weight):
+        cv.check_type('weight value', weight, Real)
+        self._weight = weight
+
     @property
     def particle(self):
         return self._particle
@@ -425,6 +442,8 @@ def populate_xml_element(self, element):
             element.append(self.energy.to_xml_element('energy'))
         if self.time is not None:
             element.append(self.time.to_xml_element('time'))
+        if self.weight is not None:
+            element.set("weight", str(self.weight))
 
     @classmethod
     def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase:
@@ -471,6 +490,10 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase:
         if time is not None:
             source.time = Univariate.from_xml_element(time)
 
+        weight = elem.find('weight')
+        if weight is not None:
+            source.weight = float(weight)
+
         return source
 
 

src/settings.cpp

@@ -567,7 +567,7 @@ void read_settings_xml(pugi::xml_node root)
     model::external_sources.push_back(make_unique<IndependentSource>(
       UPtrSpace {new SpatialPoint({0.0, 0.0, 0.0})},
       UPtrAngle {new Isotropic()}, UPtrDist {new Watt(0.988e6, 2.249e-6)},
-      UPtrDist {new Discrete(T, p, 1)}));
+      UPtrDist {new Discrete(T, p, 1)}, double {1.0}));
   }
 
   // Check if we want to write out source

src/source.cpp

@@ -171,7 +171,8 @@ SourceSite Source::sample_with_constraints(uint64_t* seed) const
       // Check whether sampled site satisfies constraints
       accepted = satisfies_spatial_constraints(site.r) &&
                  satisfies_energy_constraints(site.E) &&
-                 satisfies_time_constraints(site.time);
+                 satisfies_time_constraints(site.time) &&
+                 satisfies_weight_constraints(site.wgt);
       if (!accepted) {
         ++n_reject;
         if (n_reject >= EXTSRC_REJECT_THRESHOLD &&
@@ -207,6 +208,11 @@ bool Source::satisfies_time_constraints(double time) const
   return time > time_bounds_.first && time < time_bounds_.second;
 }
 
+bool Source::satisfies_weight_constraints(double weight) const
+{
+  return weight > weight_bounds_.first && weight < weight_bounds_.second;
+}
+
 bool Source::satisfies_spatial_constraints(Position r) const
 {
   GeometryState geom_state;
@@ -255,10 +261,11 @@ bool Source::satisfies_spatial_constraints(Position r) const
 // IndependentSource implementation
 //==============================================================================
 
-IndependentSource::IndependentSource(
-  UPtrSpace space, UPtrAngle angle, UPtrDist energy, UPtrDist time)
+IndependentSource::IndependentSource(UPtrSpace space, UPtrAngle angle,
+  UPtrDist energy, UPtrDist time, double weight)
   : space_ {std::move(space)}, angle_ {std::move(angle)},
-    energy_ {std::move(energy)}, time_ {std::move(time)}
+    energy_ {std::move(energy)}, time_ {std::move(time)}, weight_ {
+                                                            std::move(weight)}
 {}
 
 IndependentSource::IndependentSource(pugi::xml_node node) : Source(node)
@@ -324,6 +331,17 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node)
       double p[] {1.0};
       time_ = UPtrDist {new Discrete {T, p, 1}};
     }
+
+    // Determine external source weight
+    if (check_for_node(node, "weight")) {
+      weight_ = std::stod(get_node_value(node, "weight"));
+      if (weight_ < 0.0) {
+        fatal_error("Source weight is negative.");
+      }
+    } else {
+      // Default to a Constant weight wgt=1.0
+      weight_ = double {1.0};
+    }
   }
 }
 
@@ -394,6 +412,8 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
 
     // Sample particle creation time
     site.time = time_->sample(seed);
+    // Set particle creation weight
+    site.wgt = weight_;
   }
 
   // Increment number of accepted samples
@@ -568,7 +588,8 @@ SourceSite MeshSource::sample(uint64_t* seed) const
 
     // Apply other rejections
     if (satisfies_energy_constraints(site.E) &&
-        satisfies_time_constraints(site.time)) {
+        satisfies_time_constraints(site.time) &&
+        satisfies_weight_constraints(site.wgt)) {
       break;
     }
   }

tests/unit_tests/test_source_weight.py

@@ -0,0 +1,37 @@
+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),
+        weight=100
+    )
+    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
+
+