Source activity

docs/source/usersguide/settings.rst

@@ -272,6 +272,32 @@ option::
   settings.source = [src1, src2]
   settings.uniform_source_sampling = True
 
+When source strengths represent physical activity rates in Bq, the
+:class:`openmc.stats.PoissonProcess` distribution can be assigned to a source's
+time distribution to generate per-source Poisson timestamps. Each source
+generates an independent ordered sequence of exponentially-distributed
+inter-arrival times at the given rate::
+
+  src1 = openmc.IndependentSource()
+  src1.strength = 1.0e6   # 1 MBq
+  src1.time = openmc.stats.PoissonProcess(rate=1.0e6)
+  ...
+
+  src2 = openmc.IndependentSource()
+  src2.strength = 5.0e5   # 0.5 MBq
+  src2.time = openmc.stats.PoissonProcess(rate=5.0e5)
+  ...
+
+  settings.source = [src1, src2]
+
+As a convenience, :meth:`IndependentSource.set_activity` sets the time
+distribution to a Poisson process using the source strength as the rate::
+
+  src1.set_activity()  # uses src1.strength as rate
+
+This is only used in fixed-source mode. The assigned timestamps can be tallied
+using :class:`openmc.TimeFilter`.
+
 Additionally, sampling from an :class:`openmc.IndependentSource` may be biased
 for local or global variance reduction by modifying the
 :attr:`~openmc.IndependentSource.bias` attribute of each of its four main

include/openmc/distribution.h

@@ -394,6 +394,27 @@ class Mixture : public Distribution {
   double integral_;       //!< Integral of distribution
 };
 
+//==============================================================================
+//! Poisson process distribution (exponential inter-arrival times)
+//==============================================================================
+
+class PoissonProcess : public Distribution {
+public:
+  explicit PoissonProcess(pugi::xml_node node);
+  PoissonProcess(double rate) : rate_ {rate} {};
+
+  double rate() const { return rate_; }
+
+protected:
+  //! Sample a value (unbiased) from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  double sample_unbiased(uint64_t* seed) const override;
+
+private:
+  double rate_; //!< Activity rate in Bq
+};
+
 } // namespace openmc
 
 #endif // OPENMC_DISTRIBUTION_H

include/openmc/simulation.h

@@ -48,6 +48,7 @@ extern const RegularMesh* ufs_mesh;
 
 extern vector<double> k_generation;
 extern vector<int64_t> work_index;
+extern vector<double> decay_times;
 
 } // namespace simulation
 

openmc/source.py

@@ -60,7 +60,7 @@ class SourceBase(ABC):
     def __init__(
         self,
         strength: float | None = 1.0,
-        constraints: dict[str, Any] | None = None
+        constraints: dict[str, Any] | None = None,
     ):
         self.strength = strength
         self.constraints = constraints
@@ -323,7 +323,7 @@ def __init__(
         particle: str | int | ParticleType = 'neutron',
         domains: Sequence[openmc.Cell | openmc.Material |
                           openmc.Universe] | None = None,
-        constraints: dict[str, Any] | None = None
+        constraints: dict[str, Any] | None = None,
     ):
         if domains is not None:
             warnings.warn("The 'domains' arguments has been replaced by the "
@@ -412,6 +412,19 @@ def particle(self) -> ParticleType:
     def particle(self, particle):
         self._particle = ParticleType(particle)
 
+    def set_activity(self, rate=None):
+        """Set time distribution to a Poisson process.
+
+        Parameters
+        ----------
+        rate : float, optional
+            Activity rate in Bq. If not provided, uses the source strength.
+
+        """
+        if rate is None:
+            rate = self.strength
+        self.time = openmc.stats.PoissonProcess(rate)
+
     def populate_xml_element(self, element):
         """Add necessary source information to an XML element
 

openmc/stats/univariate.py

@@ -113,6 +113,8 @@ def from_xml_element(cls, elem):
             return Legendre.from_xml_element(elem)
         elif distribution == 'mixture':
             return Mixture.from_xml_element(elem)
+        elif distribution == 'poisson':
+            return PoissonProcess.from_xml_element(elem)
 
     @abstractmethod
     def _sample_unbiased(self, n_samples: int = 1, seed: int | None = None):
@@ -1948,6 +1950,92 @@ def combine_distributions(
             return Mixture(cont_probs, cont_dists)
 
 
+class PoissonProcess(Univariate):
+    """Poisson process distribution for activity-based timing.
+
+    Samples inter-arrival times from an exponential distribution with rate
+    parameter ``rate`` (in Bq). Each sample returns ``(-1/rate) * ln(1-U)``
+    where ``U`` is a uniform random variate.
+
+    Parameters
+    ----------
+    rate : float
+        Activity rate in Bq (decays per second)
+
+    Attributes
+    ----------
+    rate : float
+        Activity rate in Bq
+
+    """
+
+    def __init__(self, rate: float):
+        self.rate = rate
+        super().__init__(bias=None)
+
+    def __len__(self):
+        return 1
+
+    @property
+    def rate(self):
+        return self._rate
+
+    @rate.setter
+    def rate(self, rate):
+        cv.check_type('Poisson process rate', rate, Real)
+        cv.check_greater_than('Poisson process rate', rate, 0.0)
+        self._rate = rate
+
+    @property
+    def support(self):
+        return (0.0, float('inf'))
+
+    def _sample_unbiased(self, n_samples=1, seed=None):
+        rng = np.random.RandomState(seed)
+        return rng.exponential(1.0 / self.rate, n_samples)
+
+    def evaluate(self, x):
+        x = np.asarray(x, dtype=float)
+        return np.where(x >= 0, self.rate * np.exp(-self.rate * x), 0.0)
+
+    def to_xml_element(self, element_name: str):
+        """Return XML representation of the Poisson process distribution
+
+        Parameters
+        ----------
+        element_name : str
+            XML element name
+
+        Returns
+        -------
+        element : lxml.etree._Element
+            XML element containing Poisson process distribution data
+
+        """
+        element = ET.Element(element_name)
+        element.set("type", "poisson")
+        element.set("parameters", str(self.rate))
+        return element
+
+    @classmethod
+    def from_xml_element(cls, elem: ET.Element):
+        """Generate Poisson process distribution from an XML element
+
+        Parameters
+        ----------
+        elem : lxml.etree._Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.PoissonProcess
+            Poisson process distribution generated from XML element
+
+        """
+        params = get_elem_list(elem, "parameters", float)
+        return cls(params[0])
+
+
 def check_bias_support(parent: Univariate, bias: Univariate | None):
     """Ensure that bias distributions share the support of the univariate
     distribution they are biasing.

src/distribution.cpp

@@ -639,6 +639,20 @@ double Mixture::sample_unbiased(uint64_t* seed) const
   return distribution_[idx]->sample(seed).first;
 }
 
+//==============================================================================
+// PoissonProcess implementation
+//==============================================================================
+
+PoissonProcess::PoissonProcess(pugi::xml_node node)
+{
+  rate_ = std::stod(get_node_value(node, "parameters"));
+}
+
+double PoissonProcess::sample_unbiased(uint64_t* seed) const
+{
+  return -std::log(1.0 - prn(seed)) / rate_;
+}
+
 //==============================================================================
 // Helper function
 //==============================================================================
@@ -669,6 +683,8 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
     dist = UPtrDist {new Tabular(node)};
   } else if (type == "mixture") {
     dist = UPtrDist {new Mixture(node)};
+  } else if (type == "poisson") {
+    dist = UPtrDist {new PoissonProcess(node)};
   } else if (type == "muir") {
     openmc::fatal_error(
       "'muir' distributions are now specified using the openmc.stats.muir() "

src/simulation.cpp

@@ -4,6 +4,7 @@
 #include "openmc/capi.h"
 #include "openmc/collision_track.h"
 #include "openmc/container_util.h"
+#include "openmc/distribution.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
 #include "openmc/event.h"
@@ -40,6 +41,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <limits>
 #include <string>
 
 //==============================================================================
@@ -322,6 +324,7 @@ const RegularMesh* ufs_mesh {nullptr};
 
 vector<double> k_generation;
 vector<int64_t> work_index;
+vector<double> decay_times;
 
 } // namespace simulation
 
@@ -356,6 +359,79 @@ void allocate_banks()
   }
 }
 
+void compute_decay_times()
+{
+  int n_sources = model::external_sources.size();
+
+  // Check each source for a PoissonProcess time distribution
+  vector<double> mean_dts(n_sources, 0.0);
+  bool any_poisson = false;
+  for (int s = 0; s < n_sources; ++s) {
+    auto* indep =
+      dynamic_cast<IndependentSource*>(model::external_sources[s].get());
+    if (indep) {
+      auto* poisson = dynamic_cast<const PoissonProcess*>(indep->time());
+      if (poisson) {
+        double rate = poisson->rate();
+        mean_dts[s] = (rate > 0.0) ? 1.0 / rate : 0.0;
+        any_poisson = true;
+      }
+    }
+  }
+
+  // If no sources use PoissonProcess timing, clear and return
+  if (!any_poisson) {
+    simulation::decay_times.clear();
+    return;
+  }
+
+  int64_t global_start = simulation::work_index[mpi::rank];
+  int64_t global_end = global_start + simulation::work_per_rank;
+
+  // Per-source Poisson process state: RNG seed and cumulative time
+  vector<uint64_t> source_seeds(n_sources);
+  vector<double> cumulative_times(n_sources, 0.0);
+  for (int s = 0; s < n_sources; ++s) {
+    source_seeds[s] = init_seed(
+      simulation::current_batch * (n_sources + 1) + s + 1, STREAM_SOURCE);
+  }
+
+  // Single pass over all global particles up to this rank's end.
+  // For each particle, replicate the source selection, advance that
+  // source's Poisson clock, and store the time if it belongs to this rank.
+  // Non-Poisson sources get NaN so initialize_history() can skip them.
+  simulation::decay_times.resize(simulation::work_per_rank);
+  for (int64_t i = 0; i < global_end; ++i) {
+    // Replicate source selection RNG from sample_external_source
+    int64_t id = (simulation::total_gen + overall_generation() - 1) *
+                   settings::n_particles +
+                 i + 1;
+    uint64_t seed = init_seed(id, STREAM_SOURCE);
+
+    int src_idx;
+    if (settings::uniform_source_sampling) {
+      src_idx = static_cast<int>(prn(&seed) * n_sources);
+      if (src_idx >= n_sources)
+        src_idx = n_sources - 1;
+    } else {
+      src_idx = model::external_sources_probability.sample(&seed);
+    }
+
+    // Advance this source's Poisson clock (or mark NaN for non-Poisson)
+    double decay_time = std::numeric_limits<double>::quiet_NaN();
+    if (mean_dts[src_idx] > 0.0) {
+      cumulative_times[src_idx] +=
+        -mean_dts[src_idx] * std::log(1.0 - prn(&source_seeds[src_idx]));
+      decay_time = cumulative_times[src_idx];
+    }
+
+    // Store if this particle belongs to this rank
+    if (i >= global_start) {
+      simulation::decay_times[i - global_start] = decay_time;
+    }
+  }
+}
+
 void initialize_batch()
 {
   // Increment current batch
@@ -370,6 +446,11 @@ void initialize_batch()
     }
   }
 
+  // Pre-compute decay times for activity-based timing (fixed-source only)
+  if (settings::run_mode == RunMode::FIXED_SOURCE) {
+    compute_decay_times();
+  }
+
   // Reset total starting particle weight used for normalizing tallies
   simulation::total_weight = 0.0;
 
@@ -585,6 +666,14 @@ void initialize_history(Particle& p, int64_t index_source)
     // sample from external source distribution or custom library then set
     auto site = sample_external_source(&seed);
     p.from_source(&site);
+    // Override particle time with decay time if Poisson process timing is used
+    if (!simulation::decay_times.empty()) {
+      double decay_time = simulation::decay_times[index_source - 1];
+      if (!std::isnan(decay_time)) {
+        p.time() = decay_time;
+        p.time_last() = decay_time;
+      }
+    }
   }
   p.current_work() = index_source;
 
@@ -799,6 +888,7 @@ void free_memory_simulation()
 {
   simulation::k_generation.clear();
   simulation::entropy.clear();
+  simulation::decay_times.clear();
 }
 
 void transport_history_based_single_particle(Particle& p)