diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 5a04fedd70a..6a71e259eb2 100644 --- a/docs/source/usersguide/settings.rst +++ b/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 diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 80fe70baed8..9840564e9b8 100644 --- a/include/openmc/distribution.h +++ b/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 diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 9a6cf1b2131..13e627d60ca 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -48,6 +48,7 @@ extern const RegularMesh* ufs_mesh; extern vector k_generation; extern vector work_index; +extern vector decay_times; } // namespace simulation diff --git a/openmc/source.py b/openmc/source.py index a11cfd6c8d8..3678934ced5 100644 --- a/openmc/source.py +++ b/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 diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 94258c833c4..8c6e459a8e2 100644 --- a/openmc/stats/univariate.py +++ b/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. diff --git a/src/distribution.cpp b/src/distribution.cpp index 537a56171d9..35d05bac4cd 100644 --- a/src/distribution.cpp +++ b/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() " diff --git a/src/simulation.cpp b/src/simulation.cpp index 18e40a8bc73..f344fa553b6 100644 --- a/src/simulation.cpp +++ b/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 #include +#include #include //============================================================================== @@ -322,6 +324,7 @@ const RegularMesh* ufs_mesh {nullptr}; vector k_generation; vector work_index; +vector 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 mean_dts(n_sources, 0.0); + bool any_poisson = false; + for (int s = 0; s < n_sources; ++s) { + auto* indep = + dynamic_cast(model::external_sources[s].get()); + if (indep) { + auto* poisson = dynamic_cast(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 source_seeds(n_sources); + vector 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(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::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) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index bb8a1b78528..0298ebaec6d 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -125,6 +125,46 @@ def test_source_dlopen(): assert 'library' in elem.attrib +def test_poisson_process_xml_roundtrip(): + """Test XML round-trip for PoissonProcess distribution.""" + rate = 1.0e6 + dist = openmc.stats.PoissonProcess(rate) + elem = dist.to_xml_element('time') + + assert elem.get('type') == 'poisson' + assert float(elem.get('parameters')) == rate + + new_dist = openmc.stats.Univariate.from_xml_element(elem) + assert isinstance(new_dist, openmc.stats.PoissonProcess) + assert new_dist.rate == rate + + +def test_poisson_process_on_source(): + """Test setting PoissonProcess as time distribution on a source.""" + src = openmc.IndependentSource(strength=1.0e6) + src.time = openmc.stats.PoissonProcess(1.0e6) + elem = src.to_xml_element() + + new_src = openmc.IndependentSource.from_xml_element(elem) + assert isinstance(new_src.time, openmc.stats.PoissonProcess) + assert new_src.time.rate == 1.0e6 + + +def test_set_activity(): + """Test IndependentSource.set_activity() convenience method.""" + # Using explicit rate + src = openmc.IndependentSource(strength=500.0) + src.set_activity(rate=1.0e6) + assert isinstance(src.time, openmc.stats.PoissonProcess) + assert src.time.rate == 1.0e6 + + # Using strength as default rate + src2 = openmc.IndependentSource(strength=2.0e6) + src2.set_activity() + assert isinstance(src2.time, openmc.stats.PoissonProcess) + assert src2.time.rate == 2.0e6 + + def test_source_xml_roundtrip(): # Create a source and write to an XML element space = openmc.stats.Box([-5., -5., -5.], [5., 5., 5.])