Coincident source

docs/source/pythonapi/base.rst

@@ -26,6 +26,7 @@ Simulation Settings
    openmc.FileSource
    openmc.CompiledSource
    openmc.MeshSource
+   openmc.CoincidentSource   
    openmc.SourceParticle
    openmc.VolumeCalculation
    openmc.Settings

include/openmc/source.h

@@ -77,6 +77,23 @@ class Source {
   //! \return Sampled site
   virtual SourceSite sample(uint64_t* seed) const = 0;
 
+  //! Sample one or more source sites (without applying constraints)
+  //
+  //! For most sources this returns a single site. CoincidentSource overrides
+  //! this to return multiple coincident sites.
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Vector of sampled sites
+  virtual vector<SourceSite> sample_sites(uint64_t* seed) const
+  {
+    return {sample(seed)};
+  }
+
+  //! Sample one or more source sites and apply constraints
+  //
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Vector of sampled sites
+  vector<SourceSite> sample_sites_with_constraints(uint64_t* seed) const;
+
   static unique_ptr<Source> create(pugi::xml_node node);
 
 protected:
@@ -250,18 +267,50 @@ class MeshSource : public Source {
   vector<unique_ptr<IndependentSource>> sources_; //!< Source distributions
 };
 
+//==============================================================================
+//! Source that emits multiple coincident particles per source event
+//!
+//! All particles share the same position and time, but each has its own
+//! particle type, energy, and angular distributions. Used for simulating
+//! coincident emissions such as Co-60 summation peaks.
+//==============================================================================
+
+class CoincidentSource : public Source {
+public:
+  // Constructors
+  explicit CoincidentSource(pugi::xml_node node);
+
+  //! Sample a single source site (returns first coincident particle)
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Sampled site
+  SourceSite sample(uint64_t* seed) const override;
+
+  //! Sample all coincident source sites
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Vector of sampled sites sharing the same position and time
+  vector<SourceSite> sample_sites(uint64_t* seed) const override;
+
+private:
+  UPtrSpace space_; //!< Shared spatial distribution
+  UPtrDist time_;   //!< Shared time distribution
+  vector<unique_ptr<IndependentSource>> sources_; //!< Sub-source distributions
+  vector<double> probabilities_;     //!< Emission probability per sub-source
+  double prob_at_least_one_;         //!< P(>=1 emission) = 1 - prod(1-p_i)
+  vector<double> first_success_cdf_; //!< CDF for "first success" distribution
+};
+
 //==============================================================================
 // Functions
 //==============================================================================
 
 //! Initialize source bank from file/distribution
 extern "C" void initialize_source();
 
-//! Sample a site from all external source distributions in proportion to their
-//! source strength
+//! Sample one or more sites from all external source distributions in
+//! proportion to their source strength
 //! \param[inout] seed Pseudorandom seed pointer
-//! \return Sampled source site
-SourceSite sample_external_source(uint64_t* seed);
+//! \return Vector of sampled source sites (multiple for coincident sources)
+vector<SourceSite> sample_external_source(uint64_t* seed);
 
 void free_memory_source();
 

openmc/source.py

@@ -203,6 +203,8 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase:
                 return FileSource.from_xml_element(elem)
             elif source_type == 'mesh':
                 return MeshSource.from_xml_element(elem, meshes)
+            elif source_type == 'coincident':
+                return CoincidentSource.from_xml_element(elem, meshes)
             else:
                 raise ValueError(
                     f'Source type {source_type} is not recognized')
@@ -658,6 +660,194 @@ def from_xml_element(cls, elem: ET.Element, meshes) -> openmc.MeshSource:
         return cls(mesh, sources, constraints=constraints)
 
 
+class CoincidentSource(SourceBase):
+    """A source that emits multiple coincident particles per source event.
+
+    All particles share the same spatial position and time, but each has
+    independent particle type, energy, and angular distributions. This is
+    useful for simulating coincident emissions such as Co-60 summation peaks.
+
+    .. versionadded:: 0.15.4
+
+    Parameters
+    ----------
+    space : openmc.stats.Spatial, optional
+        Shared spatial distribution for all emitted particles
+    time : openmc.stats.Univariate, optional
+        Shared time distribution for all emitted particles
+    sources : sequence of openmc.IndependentSource
+        Sub-sources defining particle type, energy, and angle for each
+        coincident particle. Must contain at least 2 sources.
+    strength : float
+        Strength of the source
+    probabilities : sequence of float, optional
+        Emission probability for each sub-source, values in (0, 1].
+        Default is 1.0 (always emit) for every sub-source.
+    constraints : dict, optional
+        Constraints on sampled source particles.
+
+    Attributes
+    ----------
+    space : openmc.stats.Spatial or None
+        Shared spatial distribution
+    time : openmc.stats.Univariate or None
+        Shared time distribution
+    sources : list of openmc.IndependentSource
+        Sub-source distributions
+    probabilities : list of float
+        Emission probability per sub-source
+    strength : float
+        Strength of the source
+    type : str
+        Indicator of source type: 'coincident'
+
+    """
+
+    def __init__(
+        self,
+        space: openmc.stats.Spatial | None = None,
+        time: openmc.stats.Univariate | None = None,
+        sources: list[IndependentSource] | None = None,
+        strength: float = 1.0,
+        probabilities: list[float] | None = None,
+        constraints: dict[str, Any] | None = None
+    ):
+        super().__init__(strength=strength, constraints=constraints)
+
+        self._space = None
+        self._time = None
+        self._sources = []
+        self._probabilities = None
+
+        if space is not None:
+            self.space = space
+        if time is not None:
+            self.time = time
+        if sources is not None:
+            self.sources = sources
+        if probabilities is not None:
+            self.probabilities = probabilities
+
+    @property
+    def type(self) -> str:
+        return 'coincident'
+
+    @property
+    def space(self):
+        return self._space
+
+    @space.setter
+    def space(self, space):
+        cv.check_type('spatial distribution', space, Spatial)
+        self._space = space
+
+    @property
+    def time(self):
+        return self._time
+
+    @time.setter
+    def time(self, time):
+        cv.check_type('time distribution', time, Univariate)
+        self._time = time
+
+    @property
+    def sources(self):
+        return self._sources
+
+    @sources.setter
+    def sources(self, sources):
+        cv.check_type('sub-sources', sources, Iterable, IndependentSource)
+        if self._probabilities is not None:
+            n = len(self._probabilities)
+            cv.check_length('sub-sources', sources, n, n)
+        else:
+            cv.check_length('sub-sources', sources, 2)
+        self._sources = list(sources)
+
+    @property
+    def probabilities(self):
+        return self._probabilities
+
+    @probabilities.setter
+    def probabilities(self, probabilities):
+        cv.check_type('probabilities', probabilities, Iterable, Real)
+        if self._sources is not None:
+            n = len(self._sources)
+            cv.check_length('probabilities', probabilities, n, n)
+        else:
+            cv.check_length('probabilities', probabilities, 2)
+
+        for p in probabilities:
+            cv.check_less_than('probability', p, 1.0, equality=True)
+            cv.check_greater_than('probability', p, 0.0)
+        self._probabilities = list(probabilities)
+
+    def populate_xml_element(self, element):
+        """Add necessary source information to an XML element
+
+        Returns
+        -------
+        element : lxml.etree._Element
+            XML element containing source data
+
+        """
+        if self.space is not None:
+            element.append(self.space.to_xml_element())
+        if self.time is not None:
+            element.append(self.time.to_xml_element('time'))
+        probs = self.probabilities
+        for i, src in enumerate(self.sources):
+            sub_elem = src.to_xml_element()
+            if probs is not None and probs[i] != 1.0:
+                sub_elem.set("probability", str(probs[i]))
+            element.append(sub_elem)
+
+    @classmethod
+    def from_xml_element(cls, elem: ET.Element, meshes=None) -> CoincidentSource:
+        """Generate coincident source from an XML element
+
+        Parameters
+        ----------
+        elem : lxml.etree._Element
+            XML element
+        meshes : dict, optional
+            Dictionary with mesh IDs as keys and openmc.MeshBase instances as
+            values
+
+        Returns
+        -------
+        openmc.CoincidentSource
+            Source generated from XML element
+
+        """
+        constraints = cls._get_constraints(elem)
+        source = cls(constraints=constraints)
+
+        strength = get_text(elem, 'strength')
+        if strength is not None:
+            source.strength = float(strength)
+
+        space = elem.find('space')
+        if space is not None:
+            source.space = Spatial.from_xml_element(space, meshes)
+
+        time = elem.find('time')
+        if time is not None:
+            source.time = Univariate.from_xml_element(time)
+
+        sub_sources = []
+        probabilities = []
+        for e in elem.iterchildren('source'):
+            sub_sources.append(IndependentSource.from_xml_element(e))
+            probabilities.append(float(e.get('probability', 1.0)))
+        if sub_sources:
+            source.sources = sub_sources
+            if any(p != 1.0 for p in probabilities):
+                source.probabilities = probabilities
+
+        return source
+
+
 def Source(*args, **kwargs):
     """
     A function for backward compatibility of sources. Will be removed in the

src/particle.cpp

@@ -853,8 +853,8 @@ void Particle::write_restart() const
                      settings::n_particles +
                    simulation::work_index[mpi::rank] + i;
       uint64_t seed = init_seed(id, STREAM_SOURCE);
-      // re-sample source site
-      auto site = sample_external_source(&seed);
+      // re-sample source site (use primary site only)
+      auto site = sample_external_source(&seed)[0];
       write_dataset(file_id, "weight", site.wgt);
       write_dataset(file_id, "energy", site.E);
       write_dataset(file_id, "xyz", site.r);

src/simulation.cpp

@@ -583,8 +583,12 @@ void initialize_history(Particle& p, int64_t index_source)
                  simulation::work_index[mpi::rank] + index_source;
     uint64_t seed = init_seed(id, STREAM_SOURCE);
     // sample from external source distribution or custom library then set
-    auto site = sample_external_source(&seed);
-    p.from_source(&site);
+    auto sites = sample_external_source(&seed);
+    p.from_source(&sites[0]);
+    // For coincident sources, add extra particles to secondary bank
+    for (size_t i = 1; i < sites.size(); ++i) {
+      p.secondary_bank().push_back(sites[i]);
+    }
   }
   p.current_work() = index_source;