From cacf6e17bd3b56433b5f05a7625f51bf09d9c44c Mon Sep 17 00:00:00 2001 From: Christopher Date: Sun, 8 Feb 2026 10:08:56 -0500 Subject: [PATCH 1/9] sources get an activity st their time distribution is sampled from a Poisson distribution --- include/openmc/simulation.h | 1 + include/openmc/source.h | 2 ++ openmc/source.py | 32 +++++++++++++++++-- src/simulation.cpp | 55 +++++++++++++++++++++++++++++++++ src/source.cpp | 8 +++++ tests/unit_tests/test_source.py | 24 ++++++++++++++ 6 files changed, 119 insertions(+), 3 deletions(-) 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/include/openmc/source.h b/include/openmc/source.h index 2f32aa2a058..e20608dc02a 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -63,6 +63,7 @@ class Source { // Methods that can be overridden virtual double strength() const { return strength_; } + double activity() const { return activity_; } //! Sample a source site and apply constraints // @@ -97,6 +98,7 @@ class Source { // Data members double strength_ {1.0}; //!< Source strength + double activity_ {0.0}; //!< Source activity rate [Bq] std::unordered_set domain_ids_; //!< Domains to reject from DomainType domain_type_; //!< Domain type for rejection std::pair time_bounds_ {-std::numeric_limits::max(), diff --git a/openmc/source.py b/openmc/source.py index a11cfd6c8d8..427077769a6 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -30,6 +30,9 @@ class SourceBase(ABC): ---------- strength : float Strength of the source + activity : float or None + Activity rate of the source in Bq. When set, each history is assigned + a unique timestamp sampled from a Poisson process with this rate. constraints : dict Constraints on sampled source particles. Valid keys include 'domains', 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. @@ -50,6 +53,8 @@ class SourceBase(ABC): Indicator of source type. strength : float Strength of the source + activity : float or None + Activity rate of the source in Bq constraints : dict Constraints on sampled source particles. Valid keys include 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', @@ -60,10 +65,12 @@ class SourceBase(ABC): def __init__( self, strength: float | None = 1.0, - constraints: dict[str, Any] | None = None + constraints: dict[str, Any] | None = None, + activity: float | None = None ): self.strength = strength self.constraints = constraints + self.activity = activity @property def strength(self): @@ -76,6 +83,17 @@ def strength(self, strength): cv.check_greater_than('source strength', strength, 0.0, True) self._strength = strength + @property + def activity(self): + return self._activity + + @activity.setter + def activity(self, activity): + cv.check_type('source activity', activity, Real, none_ok=True) + if activity is not None: + cv.check_greater_than('source activity', activity, 0.0) + self._activity = activity + @property def constraints(self) -> dict[str, Any]: return self._constraints @@ -138,6 +156,8 @@ def to_xml_element(self) -> ET.Element: element.set("type", self.type) if self.strength is not None: element.set("strength", str(self.strength)) + if self.activity is not None: + element.set("activity", str(self.activity)) self.populate_xml_element(element) constraints = self.constraints if constraints: @@ -323,14 +343,16 @@ 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, + activity: float | None = None ): if domains is not None: warnings.warn("The 'domains' arguments has been replaced by the " "'constraints' argument.", FutureWarning) constraints = {'domains': domains} - super().__init__(strength=strength, constraints=constraints) + super().__init__(strength=strength, constraints=constraints, + activity=activity) self._space = None self._angle = None @@ -456,6 +478,10 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: if strength is not None: source.strength = float(strength) + activity = elem.get('activity') + if activity is not None: + source.activity = float(activity) + particle = get_text(elem, 'particle') if particle is not None: source.particle = particle diff --git a/src/simulation.cpp b/src/simulation.cpp index 18e40a8bc73..c74f5e2f98d 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -322,6 +322,7 @@ const RegularMesh* ufs_mesh {nullptr}; vector k_generation; vector work_index; +vector decay_times; } // namespace simulation @@ -356,6 +357,48 @@ void allocate_banks() } } +void compute_decay_times() +{ + // Sum activities across all external sources + double total_activity = 0.0; + for (const auto& src : model::external_sources) { + total_activity += src->activity(); + } + + // If no activity is set, clear decay times and return + if (total_activity <= 0.0) { + simulation::decay_times.clear(); + return; + } + + // Mean inter-arrival time in seconds + double mean_dt = 1.0 / total_activity; + + // Compute global start index for this MPI rank + int64_t global_start = simulation::work_index[mpi::rank]; + + // Use a deterministic seed based on batch number + uint64_t seed = init_seed(simulation::current_batch, STREAM_SOURCE); + + // Generate exponential inter-arrival times sequentially from the start + // to ensure reproducibility across MPI ranks + double cumulative_time = 0.0; + + // Skip to this rank's starting index + for (int64_t i = 0; i < global_start; ++i) { + double dt = -mean_dt * std::log(1.0 - prn(&seed)); + cumulative_time += dt; + } + + // Compute decay times for this rank's particles + simulation::decay_times.resize(simulation::work_per_rank); + for (int64_t i = 0; i < simulation::work_per_rank; ++i) { + double dt = -mean_dt * std::log(1.0 - prn(&seed)); + cumulative_time += dt; + simulation::decay_times[i] = cumulative_time; + } +} + void initialize_batch() { // Increment current batch @@ -370,6 +413,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 +633,12 @@ 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 activity is set + if (!simulation::decay_times.empty()) { + double decay_time = simulation::decay_times[index_source - 1]; + p.time() = decay_time; + p.time_last() = decay_time; + } } p.current_work() = index_source; @@ -799,6 +853,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/src/source.cpp b/src/source.cpp index 8bd9c789352..57a19a652a6 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -77,6 +77,14 @@ Source::Source(pugi::xml_node node) } } + // Check for source activity rate + if (check_for_node(node, "activity")) { + activity_ = std::stod(get_node_value(node, "activity")); + if (activity_ <= 0.0) { + fatal_error("Source activity must be positive."); + } + } + // Check for additional defined constraints read_constraints(node); } diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index bb8a1b78528..e0accb9b9df 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -125,6 +125,30 @@ def test_source_dlopen(): assert 'library' in elem.attrib +def test_activity_independent_source_roundtrip(): + """Test XML round-trip for IndependentSource with activity.""" + src = openmc.IndependentSource( + energy=openmc.stats.Discrete([1.0e6], [1.0]), + strength=2.0, + activity=1.0e6 + ) + elem = src.to_xml_element() + assert elem.get('activity') == str(1.0e6) + + new_src = openmc.IndependentSource.from_xml_element(elem) + assert new_src.activity == approx(1.0e6) + assert new_src.strength == approx(2.0) + + +def test_activity_none_by_default(): + """Test that activity is None by default.""" + src = openmc.IndependentSource() + assert src.activity is None + elem = src.to_xml_element() + assert elem.get('activity') is None + + + def test_source_xml_roundtrip(): # Create a source and write to an XML element space = openmc.stats.Box([-5., -5., -5.], [5., 5., 5.]) From b63f99f0561af83dca9c1aec479ebd605abaa2c2 Mon Sep 17 00:00:00 2001 From: Christopher Date: Wed, 11 Feb 2026 10:02:21 +0100 Subject: [PATCH 2/9] use source strength as activity --- include/openmc/settings.h | 2 + include/openmc/source.h | 2 - openmc/settings.py | 25 +++++++ openmc/source.py | 28 +------- src/settings.cpp | 8 +++ src/simulation.cpp | 115 ++++++++++++++++++++++++-------- src/source.cpp | 8 --- tests/unit_tests/test_source.py | 36 ++++------ 8 files changed, 136 insertions(+), 88 deletions(-) diff --git a/include/openmc/settings.h b/include/openmc/settings.h index b369c99fef8..57ba1239926 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -104,6 +104,8 @@ extern bool weight_window_checkpoint_collision; //!< enable weight window check //!< upon collision? extern bool write_all_tracks; //!< write track files for every particle? extern bool write_initial_source; //!< write out initial source file? +extern bool activity_based_timing; //!< use source strengths as activity rates + //!< for Poisson timestamp generation // Paths to various files extern std::string path_cross_sections; //!< path to cross_sections.xml diff --git a/include/openmc/source.h b/include/openmc/source.h index e20608dc02a..2f32aa2a058 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -63,7 +63,6 @@ class Source { // Methods that can be overridden virtual double strength() const { return strength_; } - double activity() const { return activity_; } //! Sample a source site and apply constraints // @@ -98,7 +97,6 @@ class Source { // Data members double strength_ {1.0}; //!< Source strength - double activity_ {0.0}; //!< Source activity rate [Bq] std::unordered_set domain_ids_; //!< Domains to reject from DomainType domain_type_; //!< Domain type for rejection std::pair time_bounds_ {-std::numeric_limits::max(), diff --git a/openmc/settings.py b/openmc/settings.py index 3cf662e311e..efaff1b1a1c 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -327,6 +327,9 @@ class Settings: uniform_source_sampling : bool Whether to sampling among multiple sources uniformly, applying their strengths as weights to sampled particles. + activity_based_timing : bool + Whether to use source strengths as activity rates (Bq) for generating + per-source Poisson timestamps assigned to each particle. ufs_mesh : openmc.RegularMesh Mesh to be used for redistributing source sites via the uniform fission site (UFS) method. @@ -393,6 +396,7 @@ def __init__(self, **kwargs): self._plot_seed = None self._ptables = None self._uniform_source_sampling = None + self._activity_based_timing = None self._seed = None self._stride = None self._survival_biasing = None @@ -676,6 +680,15 @@ def uniform_source_sampling(self, uniform_source_sampling: bool): cv.check_type('strength as weights', uniform_source_sampling, bool) self._uniform_source_sampling = uniform_source_sampling + @property + def activity_based_timing(self) -> bool: + return self._activity_based_timing + + @activity_based_timing.setter + def activity_based_timing(self, activity_based_timing: bool): + cv.check_type('activity-based timing', activity_based_timing, bool) + self._activity_based_timing = activity_based_timing + @property def plot_seed(self): return self._plot_seed @@ -1498,6 +1511,11 @@ def _create_uniform_source_sampling_subelement(self, root): element = ET.SubElement(root, "uniform_source_sampling") element.text = str(self._uniform_source_sampling).lower() + def _create_activity_based_timing_subelement(self, root): + if self._activity_based_timing is not None: + element = ET.SubElement(root, "activity_based_timing") + element.text = str(self._activity_based_timing).lower() + def _create_sourcepoint_subelement(self, root): if self._sourcepoint: element = ET.SubElement(root, "source_point") @@ -2104,6 +2122,11 @@ def _uniform_source_sampling_from_xml_element(self, root): if text is not None: self.uniform_source_sampling = text in ('true', '1') + def _activity_based_timing_from_xml_element(self, root): + text = get_text(root, 'activity_based_timing') + if text is not None: + self.activity_based_timing = text in ('true', '1') + def _plot_seed_from_xml_element(self, root): text = get_text(root, 'plot_seed') if text is not None: @@ -2427,6 +2450,7 @@ def to_xml_element(self, mesh_memo=None): self._create_max_order_subelement(element) self._create_photon_transport_subelement(element) self._create_uniform_source_sampling_subelement(element) + self._create_activity_based_timing_subelement(element) self._create_plot_seed_subelement(element) self._create_ptables_subelement(element) self._create_seed_subelement(element) @@ -2541,6 +2565,7 @@ def from_xml_element(cls, elem, meshes=None): settings._max_order_from_xml_element(elem) settings._photon_transport_from_xml_element(elem) settings._uniform_source_sampling_from_xml_element(elem) + settings._activity_based_timing_from_xml_element(elem) settings._plot_seed_from_xml_element(elem) settings._ptables_from_xml_element(elem) settings._seed_from_xml_element(elem) diff --git a/openmc/source.py b/openmc/source.py index 427077769a6..08ada7298f2 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -30,9 +30,6 @@ class SourceBase(ABC): ---------- strength : float Strength of the source - activity : float or None - Activity rate of the source in Bq. When set, each history is assigned - a unique timestamp sampled from a Poisson process with this rate. constraints : dict Constraints on sampled source particles. Valid keys include 'domains', 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. @@ -53,8 +50,6 @@ class SourceBase(ABC): Indicator of source type. strength : float Strength of the source - activity : float or None - Activity rate of the source in Bq constraints : dict Constraints on sampled source particles. Valid keys include 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', @@ -66,11 +61,9 @@ def __init__( self, strength: float | None = 1.0, constraints: dict[str, Any] | None = None, - activity: float | None = None ): self.strength = strength self.constraints = constraints - self.activity = activity @property def strength(self): @@ -83,17 +76,6 @@ def strength(self, strength): cv.check_greater_than('source strength', strength, 0.0, True) self._strength = strength - @property - def activity(self): - return self._activity - - @activity.setter - def activity(self, activity): - cv.check_type('source activity', activity, Real, none_ok=True) - if activity is not None: - cv.check_greater_than('source activity', activity, 0.0) - self._activity = activity - @property def constraints(self) -> dict[str, Any]: return self._constraints @@ -156,8 +138,6 @@ def to_xml_element(self) -> ET.Element: element.set("type", self.type) if self.strength is not None: element.set("strength", str(self.strength)) - if self.activity is not None: - element.set("activity", str(self.activity)) self.populate_xml_element(element) constraints = self.constraints if constraints: @@ -344,15 +324,13 @@ def __init__( domains: Sequence[openmc.Cell | openmc.Material | openmc.Universe] | None = None, constraints: dict[str, Any] | None = None, - activity: float | None = None ): if domains is not None: warnings.warn("The 'domains' arguments has been replaced by the " "'constraints' argument.", FutureWarning) constraints = {'domains': domains} - super().__init__(strength=strength, constraints=constraints, - activity=activity) + super().__init__(strength=strength, constraints=constraints) self._space = None self._angle = None @@ -478,10 +456,6 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: if strength is not None: source.strength = float(strength) - activity = elem.get('activity') - if activity is not None: - source.activity = float(activity) - particle = get_text(elem, 'particle') if particle is not None: source.particle = particle diff --git a/src/settings.cpp b/src/settings.cpp index 5b472468fc4..54ff4688e2c 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -87,6 +87,7 @@ bool weight_window_checkpoint_surface {false}; bool weight_window_checkpoint_collision {true}; bool write_all_tracks {false}; bool write_initial_source {false}; +bool activity_based_timing {false}; std::string path_cross_sections; std::string path_input; @@ -887,6 +888,13 @@ void read_settings_xml(pugi::xml_node root) get_node_value_bool(root, "uniform_source_sampling"); } + // Check if the user wants activity-based timing (per-source Poisson + // timestamps using strength as activity rate in Bq) + if (check_for_node(root, "activity_based_timing")) { + activity_based_timing = + get_node_value_bool(root, "activity_based_timing"); + } + // Check if the user has specified to write surface source if (check_for_node(root, "surf_source_write")) { surf_source_write = true; diff --git a/src/simulation.cpp b/src/simulation.cpp index c74f5e2f98d..c33362ddf78 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -359,43 +359,100 @@ void allocate_banks() void compute_decay_times() { - // Sum activities across all external sources - double total_activity = 0.0; - for (const auto& src : model::external_sources) { - total_activity += src->activity(); - } - - // If no activity is set, clear decay times and return - if (total_activity <= 0.0) { + // If activity-based timing is not enabled, clear and return + if (!settings::activity_based_timing) { simulation::decay_times.clear(); return; } - // Mean inter-arrival time in seconds - double mean_dt = 1.0 / total_activity; + int n_sources = model::external_sources.size(); + int64_t n_particles = simulation::work_per_rank; - // Compute global start index for this MPI rank - int64_t global_start = simulation::work_index[mpi::rank]; + // Pass 1: Determine which source each particle comes from by replicating + // the source selection RNG logic from sample_external_source + vector source_index(n_particles); + vector n_per_source(n_sources, 0); - // Use a deterministic seed based on batch number - uint64_t seed = init_seed(simulation::current_batch, STREAM_SOURCE); + for (int64_t i = 0; i < n_particles; ++i) { + // Replicate the seed derivation from initialize_history (line 629-632) + int64_t id = (simulation::total_gen + overall_generation() - 1) * + settings::n_particles + + simulation::work_index[mpi::rank] + (i + 1); + uint64_t seed = init_seed(id, STREAM_SOURCE); + + // Source selection is the first RNG call in sample_external_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); + } + source_index[i] = src_idx; + n_per_source[src_idx]++; + } + + // Pass 2: For each source, generate cumulative exponential inter-arrival + // times with rate = strength (activity in Bq), using a deterministic + // per-source seed + vector> per_source_times(n_sources); + for (int s = 0; s < n_sources; ++s) { + if (n_per_source[s] == 0) continue; + + double activity = model::external_sources[s]->strength(); + if (activity <= 0.0) { + // If strength is zero or negative, assign time 0 + per_source_times[s].assign(n_per_source[s], 0.0); + continue; + } + double mean_dt = 1.0 / activity; + + // Deterministic per-source seed based on batch and source index + uint64_t seed = init_seed( + simulation::current_batch * (n_sources + 1) + s + 1, STREAM_SOURCE); + + // Compute global start index for this rank's contribution to this source + // We need to skip past particles assigned to this source on earlier ranks + int64_t global_skip = 0; + for (int64_t rank = 0; rank < mpi::rank; ++rank) { + int64_t rank_start = simulation::work_index[rank]; + int64_t rank_end = simulation::work_index[rank + 1]; + for (int64_t j = rank_start; j < rank_end; ++j) { + // Replicate source selection for particles on earlier ranks + int64_t pid = (simulation::total_gen + overall_generation() - 1) * + settings::n_particles + j + 1; + uint64_t pseed = init_seed(pid, STREAM_SOURCE); + int pidx; + if (settings::uniform_source_sampling) { + pidx = static_cast(prn(&pseed) * n_sources); + if (pidx >= n_sources) pidx = n_sources - 1; + } else { + pidx = model::external_sources_probability.sample(&pseed); + } + if (pidx == s) global_skip++; + } + } - // Generate exponential inter-arrival times sequentially from the start - // to ensure reproducibility across MPI ranks - double cumulative_time = 0.0; + // Generate and skip cumulative times for earlier ranks + double cumulative_time = 0.0; + for (int64_t k = 0; k < global_skip; ++k) { + cumulative_time += -mean_dt * std::log(1.0 - prn(&seed)); + } - // Skip to this rank's starting index - for (int64_t i = 0; i < global_start; ++i) { - double dt = -mean_dt * std::log(1.0 - prn(&seed)); - cumulative_time += dt; + // Generate times for this rank's particles assigned to this source + per_source_times[s].resize(n_per_source[s]); + for (int64_t k = 0; k < n_per_source[s]; ++k) { + cumulative_time += -mean_dt * std::log(1.0 - prn(&seed)); + per_source_times[s][k] = cumulative_time; + } } - // Compute decay times for this rank's particles - simulation::decay_times.resize(simulation::work_per_rank); - for (int64_t i = 0; i < simulation::work_per_rank; ++i) { - double dt = -mean_dt * std::log(1.0 - prn(&seed)); - cumulative_time += dt; - simulation::decay_times[i] = cumulative_time; + // Pass 3: Assign per-particle timestamps using per-source counters + simulation::decay_times.resize(n_particles); + vector source_counter(n_sources, 0); + for (int64_t i = 0; i < n_particles; ++i) { + int s = source_index[i]; + simulation::decay_times[i] = per_source_times[s][source_counter[s]++]; } } @@ -633,8 +690,8 @@ 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 activity is set - if (!simulation::decay_times.empty()) { + // Override particle time with decay time if activity-based timing is on + if (settings::activity_based_timing) { double decay_time = simulation::decay_times[index_source - 1]; p.time() = decay_time; p.time_last() = decay_time; diff --git a/src/source.cpp b/src/source.cpp index 57a19a652a6..8bd9c789352 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -77,14 +77,6 @@ Source::Source(pugi::xml_node node) } } - // Check for source activity rate - if (check_for_node(node, "activity")) { - activity_ = std::stod(get_node_value(node, "activity")); - if (activity_ <= 0.0) { - fatal_error("Source activity must be positive."); - } - } - // Check for additional defined constraints read_constraints(node); } diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index e0accb9b9df..b94a8d79d3b 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -125,28 +125,20 @@ def test_source_dlopen(): assert 'library' in elem.attrib -def test_activity_independent_source_roundtrip(): - """Test XML round-trip for IndependentSource with activity.""" - src = openmc.IndependentSource( - energy=openmc.stats.Discrete([1.0e6], [1.0]), - strength=2.0, - activity=1.0e6 - ) - elem = src.to_xml_element() - assert elem.get('activity') == str(1.0e6) - - new_src = openmc.IndependentSource.from_xml_element(elem) - assert new_src.activity == approx(1.0e6) - assert new_src.strength == approx(2.0) - - -def test_activity_none_by_default(): - """Test that activity is None by default.""" - src = openmc.IndependentSource() - assert src.activity is None - elem = src.to_xml_element() - assert elem.get('activity') is None - +def test_activity_based_timing_setting_roundtrip(): + """Test XML round-trip for activity_based_timing setting.""" + settings = openmc.Settings() + settings.activity_based_timing = True + settings.particles = 100 + settings.batches = 1 + settings.run_mode = 'fixed source' + elem = settings.to_xml_element() + text = elem.find('activity_based_timing') + assert text is not None + assert text.text == 'true' + + new_settings = openmc.Settings.from_xml_element(elem) + assert new_settings.activity_based_timing is True def test_source_xml_roundtrip(): From 2e60ad34906126d95103f781e8deabaa710535ee Mon Sep 17 00:00:00 2001 From: Christopher Date: Wed, 11 Feb 2026 10:03:39 +0100 Subject: [PATCH 3/9] clang and typo --- include/openmc/settings.h | 2 +- src/settings.cpp | 3 +-- src/simulation.cpp | 15 ++++++++++----- 3 files changed, 12 insertions(+), 8 deletions(-) diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 57ba1239926..627742aca27 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -105,7 +105,7 @@ extern bool weight_window_checkpoint_collision; //!< enable weight window check extern bool write_all_tracks; //!< write track files for every particle? extern bool write_initial_source; //!< write out initial source file? extern bool activity_based_timing; //!< use source strengths as activity rates - //!< for Poisson timestamp generation + //!< for Poisson timestamp generation? // Paths to various files extern std::string path_cross_sections; //!< path to cross_sections.xml diff --git a/src/settings.cpp b/src/settings.cpp index 54ff4688e2c..9231ffd978f 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -891,8 +891,7 @@ void read_settings_xml(pugi::xml_node root) // Check if the user wants activity-based timing (per-source Poisson // timestamps using strength as activity rate in Bq) if (check_for_node(root, "activity_based_timing")) { - activity_based_timing = - get_node_value_bool(root, "activity_based_timing"); + activity_based_timing = get_node_value_bool(root, "activity_based_timing"); } // Check if the user has specified to write surface source diff --git a/src/simulation.cpp b/src/simulation.cpp index c33362ddf78..105d5f2cc04 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -384,7 +384,8 @@ void compute_decay_times() 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; + if (src_idx >= n_sources) + src_idx = n_sources - 1; } else { src_idx = model::external_sources_probability.sample(&seed); } @@ -397,7 +398,8 @@ void compute_decay_times() // per-source seed vector> per_source_times(n_sources); for (int s = 0; s < n_sources; ++s) { - if (n_per_source[s] == 0) continue; + if (n_per_source[s] == 0) + continue; double activity = model::external_sources[s]->strength(); if (activity <= 0.0) { @@ -420,16 +422,19 @@ void compute_decay_times() for (int64_t j = rank_start; j < rank_end; ++j) { // Replicate source selection for particles on earlier ranks int64_t pid = (simulation::total_gen + overall_generation() - 1) * - settings::n_particles + j + 1; + settings::n_particles + + j + 1; uint64_t pseed = init_seed(pid, STREAM_SOURCE); int pidx; if (settings::uniform_source_sampling) { pidx = static_cast(prn(&pseed) * n_sources); - if (pidx >= n_sources) pidx = n_sources - 1; + if (pidx >= n_sources) + pidx = n_sources - 1; } else { pidx = model::external_sources_probability.sample(&pseed); } - if (pidx == s) global_skip++; + if (pidx == s) + global_skip++; } } From 83b9188dbe8a4b8cebe88322836334664666640d Mon Sep 17 00:00:00 2001 From: Christopher Date: Wed, 11 Feb 2026 10:10:05 +0100 Subject: [PATCH 4/9] shorten --- src/simulation.cpp | 100 ++++++++++++--------------------------------- 1 file changed, 26 insertions(+), 74 deletions(-) diff --git a/src/simulation.cpp b/src/simulation.cpp index 105d5f2cc04..74fc470964c 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -366,21 +366,31 @@ void compute_decay_times() } int n_sources = model::external_sources.size(); - int64_t n_particles = simulation::work_per_rank; + int64_t global_start = simulation::work_index[mpi::rank]; + int64_t global_end = global_start + simulation::work_per_rank; - // Pass 1: Determine which source each particle comes from by replicating - // the source selection RNG logic from sample_external_source - vector source_index(n_particles); - vector n_per_source(n_sources, 0); + // Per-source Poisson process state: RNG seed and cumulative time + vector source_seeds(n_sources); + vector cumulative_times(n_sources, 0.0); + vector mean_dts(n_sources); + for (int s = 0; s < n_sources; ++s) { + source_seeds[s] = init_seed( + simulation::current_batch * (n_sources + 1) + s + 1, STREAM_SOURCE); + double activity = model::external_sources[s]->strength(); + mean_dts[s] = (activity > 0.0) ? 1.0 / activity : 0.0; + } - for (int64_t i = 0; i < n_particles; ++i) { - // Replicate the seed derivation from initialize_history (line 629-632) + // 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. + 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 + - simulation::work_index[mpi::rank] + (i + 1); + i + 1; uint64_t seed = init_seed(id, STREAM_SOURCE); - // Source selection is the first RNG call in sample_external_source int src_idx; if (settings::uniform_source_sampling) { src_idx = static_cast(prn(&seed) * n_sources); @@ -389,75 +399,17 @@ void compute_decay_times() } else { src_idx = model::external_sources_probability.sample(&seed); } - source_index[i] = src_idx; - n_per_source[src_idx]++; - } - - // Pass 2: For each source, generate cumulative exponential inter-arrival - // times with rate = strength (activity in Bq), using a deterministic - // per-source seed - vector> per_source_times(n_sources); - for (int s = 0; s < n_sources; ++s) { - if (n_per_source[s] == 0) - continue; - - double activity = model::external_sources[s]->strength(); - if (activity <= 0.0) { - // If strength is zero or negative, assign time 0 - per_source_times[s].assign(n_per_source[s], 0.0); - continue; - } - double mean_dt = 1.0 / activity; - - // Deterministic per-source seed based on batch and source index - uint64_t seed = init_seed( - simulation::current_batch * (n_sources + 1) + s + 1, STREAM_SOURCE); - // Compute global start index for this rank's contribution to this source - // We need to skip past particles assigned to this source on earlier ranks - int64_t global_skip = 0; - for (int64_t rank = 0; rank < mpi::rank; ++rank) { - int64_t rank_start = simulation::work_index[rank]; - int64_t rank_end = simulation::work_index[rank + 1]; - for (int64_t j = rank_start; j < rank_end; ++j) { - // Replicate source selection for particles on earlier ranks - int64_t pid = (simulation::total_gen + overall_generation() - 1) * - settings::n_particles + - j + 1; - uint64_t pseed = init_seed(pid, STREAM_SOURCE); - int pidx; - if (settings::uniform_source_sampling) { - pidx = static_cast(prn(&pseed) * n_sources); - if (pidx >= n_sources) - pidx = n_sources - 1; - } else { - pidx = model::external_sources_probability.sample(&pseed); - } - if (pidx == s) - global_skip++; - } + // Advance this source's Poisson clock + if (mean_dts[src_idx] > 0.0) { + cumulative_times[src_idx] += + -mean_dts[src_idx] * std::log(1.0 - prn(&source_seeds[src_idx])); } - // Generate and skip cumulative times for earlier ranks - double cumulative_time = 0.0; - for (int64_t k = 0; k < global_skip; ++k) { - cumulative_time += -mean_dt * std::log(1.0 - prn(&seed)); + // Store if this particle belongs to this rank + if (i >= global_start) { + simulation::decay_times[i - global_start] = cumulative_times[src_idx]; } - - // Generate times for this rank's particles assigned to this source - per_source_times[s].resize(n_per_source[s]); - for (int64_t k = 0; k < n_per_source[s]; ++k) { - cumulative_time += -mean_dt * std::log(1.0 - prn(&seed)); - per_source_times[s][k] = cumulative_time; - } - } - - // Pass 3: Assign per-particle timestamps using per-source counters - simulation::decay_times.resize(n_particles); - vector source_counter(n_sources, 0); - for (int64_t i = 0; i < n_particles; ++i) { - int s = source_index[i]; - simulation::decay_times[i] = per_source_times[s][source_counter[s]++]; } } From ee5b4c070ffab2f9d0693624facbe540ed805273 Mon Sep 17 00:00:00 2001 From: Christopher Date: Wed, 11 Feb 2026 10:23:15 +0100 Subject: [PATCH 5/9] documentation --- docs/source/io_formats/settings.rst | 12 ++++++++++++ docs/source/usersguide/settings.rst | 21 +++++++++++++++++++++ 2 files changed, 33 insertions(+) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index ed7c6273aaa..47bb681aa21 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -1418,6 +1418,18 @@ particles. *Default*: False +------------------------------------- +```` Element +------------------------------------- + +The ```` element indicates whether to treat each source's +strength as an activity rate in Bq and assign each particle a timestamp from a +per-source Poisson process. When enabled, each source generates an independent +ordered sequence of exponentially-distributed inter-arrival times at rate equal +to its strength. This is only used in fixed-source mode. + + *Default*: False + ------------------------ ```` Element ------------------------ diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 5a04fedd70a..32aaa37dda6 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -272,6 +272,27 @@ option:: settings.source = [src1, src2] settings.uniform_source_sampling = True +When source strengths represent physical activity rates in Bq, the +:attr:`Settings.activity_based_timing` attribute can be used to assign each +particle a timestamp from an independent per-source Poisson process. Each +source's strength is treated as its decay rate, and particles receive ordered +timestamps with exponentially-distributed inter-arrival times:: + + src1 = openmc.IndependentSource() + src1.strength = 1.0e6 # 1 MBq + ... + + src2 = openmc.IndependentSource() + src2.strength = 5.0e5 # 0.5 MBq + ... + + settings.source = [src1, src2] + settings.activity_based_timing = True + +This is only used in fixed-source mode. The assigned timestamps override any +time distribution set on the source and 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 From 3de2714afbf861c519c83412b6a8dfd5f508b953 Mon Sep 17 00:00:00 2001 From: Christopher Date: Wed, 11 Feb 2026 10:59:10 +0100 Subject: [PATCH 6/9] clang --- include/openmc/settings.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 627742aca27..6ce46fa64e5 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -102,8 +102,8 @@ extern bool weight_window_checkpoint_surface; //!< enable weight window check //!< upon surface crossing? extern bool weight_window_checkpoint_collision; //!< enable weight window check //!< upon collision? -extern bool write_all_tracks; //!< write track files for every particle? -extern bool write_initial_source; //!< write out initial source file? +extern bool write_all_tracks; //!< write track files for every particle? +extern bool write_initial_source; //!< write out initial source file? extern bool activity_based_timing; //!< use source strengths as activity rates //!< for Poisson timestamp generation? From 8aa476ef7cf6da7bd4ecceda0d2cb1aa3a34fdc5 Mon Sep 17 00:00:00 2001 From: Christopher Date: Thu, 12 Feb 2026 15:10:30 +0100 Subject: [PATCH 7/9] new activity approach based on Poisson as a new stochastic process --- docs/source/io_formats/settings.rst | 12 ---- docs/source/usersguide/settings.rst | 21 ++++--- include/openmc/distribution.h | 21 +++++++ include/openmc/settings.h | 4 +- openmc/settings.py | 25 -------- openmc/source.py | 13 +++++ openmc/stats/univariate.py | 88 +++++++++++++++++++++++++++++ src/distribution.cpp | 16 ++++++ src/settings.cpp | 8 +-- src/simulation.cpp | 45 +++++++++++---- tests/unit_tests/test_source.py | 52 ++++++++++++----- 11 files changed, 224 insertions(+), 81 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 47bb681aa21..ed7c6273aaa 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -1418,18 +1418,6 @@ particles. *Default*: False -------------------------------------- -```` Element -------------------------------------- - -The ```` element indicates whether to treat each source's -strength as an activity rate in Bq and assign each particle a timestamp from a -per-source Poisson process. When enabled, each source generates an independent -ordered sequence of exponentially-distributed inter-arrival times at rate equal -to its strength. This is only used in fixed-source mode. - - *Default*: False - ------------------------ ```` Element ------------------------ diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 32aaa37dda6..6a71e259eb2 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -273,25 +273,30 @@ option:: settings.uniform_source_sampling = True When source strengths represent physical activity rates in Bq, the -:attr:`Settings.activity_based_timing` attribute can be used to assign each -particle a timestamp from an independent per-source Poisson process. Each -source's strength is treated as its decay rate, and particles receive ordered -timestamps with exponentially-distributed inter-arrival times:: +: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] - settings.activity_based_timing = True -This is only used in fixed-source mode. The assigned timestamps override any -time distribution set on the source and can be tallied using -:class:`openmc.TimeFilter`. +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 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/settings.h b/include/openmc/settings.h index 6ce46fa64e5..8d0503190c8 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -103,9 +103,7 @@ extern bool weight_window_checkpoint_surface; //!< enable weight window check extern bool weight_window_checkpoint_collision; //!< enable weight window check //!< upon collision? extern bool write_all_tracks; //!< write track files for every particle? -extern bool write_initial_source; //!< write out initial source file? -extern bool activity_based_timing; //!< use source strengths as activity rates - //!< for Poisson timestamp generation? +extern bool write_initial_source; //!< write out initial source file? // Paths to various files extern std::string path_cross_sections; //!< path to cross_sections.xml diff --git a/openmc/settings.py b/openmc/settings.py index efaff1b1a1c..3cf662e311e 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -327,9 +327,6 @@ class Settings: uniform_source_sampling : bool Whether to sampling among multiple sources uniformly, applying their strengths as weights to sampled particles. - activity_based_timing : bool - Whether to use source strengths as activity rates (Bq) for generating - per-source Poisson timestamps assigned to each particle. ufs_mesh : openmc.RegularMesh Mesh to be used for redistributing source sites via the uniform fission site (UFS) method. @@ -396,7 +393,6 @@ def __init__(self, **kwargs): self._plot_seed = None self._ptables = None self._uniform_source_sampling = None - self._activity_based_timing = None self._seed = None self._stride = None self._survival_biasing = None @@ -680,15 +676,6 @@ def uniform_source_sampling(self, uniform_source_sampling: bool): cv.check_type('strength as weights', uniform_source_sampling, bool) self._uniform_source_sampling = uniform_source_sampling - @property - def activity_based_timing(self) -> bool: - return self._activity_based_timing - - @activity_based_timing.setter - def activity_based_timing(self, activity_based_timing: bool): - cv.check_type('activity-based timing', activity_based_timing, bool) - self._activity_based_timing = activity_based_timing - @property def plot_seed(self): return self._plot_seed @@ -1511,11 +1498,6 @@ def _create_uniform_source_sampling_subelement(self, root): element = ET.SubElement(root, "uniform_source_sampling") element.text = str(self._uniform_source_sampling).lower() - def _create_activity_based_timing_subelement(self, root): - if self._activity_based_timing is not None: - element = ET.SubElement(root, "activity_based_timing") - element.text = str(self._activity_based_timing).lower() - def _create_sourcepoint_subelement(self, root): if self._sourcepoint: element = ET.SubElement(root, "source_point") @@ -2122,11 +2104,6 @@ def _uniform_source_sampling_from_xml_element(self, root): if text is not None: self.uniform_source_sampling = text in ('true', '1') - def _activity_based_timing_from_xml_element(self, root): - text = get_text(root, 'activity_based_timing') - if text is not None: - self.activity_based_timing = text in ('true', '1') - def _plot_seed_from_xml_element(self, root): text = get_text(root, 'plot_seed') if text is not None: @@ -2450,7 +2427,6 @@ def to_xml_element(self, mesh_memo=None): self._create_max_order_subelement(element) self._create_photon_transport_subelement(element) self._create_uniform_source_sampling_subelement(element) - self._create_activity_based_timing_subelement(element) self._create_plot_seed_subelement(element) self._create_ptables_subelement(element) self._create_seed_subelement(element) @@ -2565,7 +2541,6 @@ def from_xml_element(cls, elem, meshes=None): settings._max_order_from_xml_element(elem) settings._photon_transport_from_xml_element(elem) settings._uniform_source_sampling_from_xml_element(elem) - settings._activity_based_timing_from_xml_element(elem) settings._plot_seed_from_xml_element(elem) settings._ptables_from_xml_element(elem) settings._seed_from_xml_element(elem) diff --git a/openmc/source.py b/openmc/source.py index 08ada7298f2..3678934ced5 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -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/settings.cpp b/src/settings.cpp index 9231ffd978f..fb17493241d 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -87,7 +87,7 @@ bool weight_window_checkpoint_surface {false}; bool weight_window_checkpoint_collision {true}; bool write_all_tracks {false}; bool write_initial_source {false}; -bool activity_based_timing {false}; + std::string path_cross_sections; std::string path_input; @@ -888,12 +888,6 @@ void read_settings_xml(pugi::xml_node root) get_node_value_bool(root, "uniform_source_sampling"); } - // Check if the user wants activity-based timing (per-source Poisson - // timestamps using strength as activity rate in Bq) - if (check_for_node(root, "activity_based_timing")) { - activity_based_timing = get_node_value_bool(root, "activity_based_timing"); - } - // Check if the user has specified to write surface source if (check_for_node(root, "surf_source_write")) { surf_source_write = true; diff --git a/src/simulation.cpp b/src/simulation.cpp index 74fc470964c..ebc94edaaea 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -2,6 +2,7 @@ #include "openmc/bank.h" #include "openmc/capi.h" +#include "openmc/distribution.h" #include "openmc/collision_track.h" #include "openmc/container_util.h" #include "openmc/eigenvalue.h" @@ -40,6 +41,7 @@ #include #include +#include #include //============================================================================== @@ -359,30 +361,45 @@ void allocate_banks() void compute_decay_times() { - // If activity-based timing is not enabled, clear and return - if (!settings::activity_based_timing) { + 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; } - int n_sources = model::external_sources.size(); 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); - vector mean_dts(n_sources); for (int s = 0; s < n_sources; ++s) { source_seeds[s] = init_seed( simulation::current_batch * (n_sources + 1) + s + 1, STREAM_SOURCE); - double activity = model::external_sources[s]->strength(); - mean_dts[s] = (activity > 0.0) ? 1.0 / activity : 0.0; } // 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 @@ -400,15 +417,17 @@ void compute_decay_times() src_idx = model::external_sources_probability.sample(&seed); } - // Advance this source's Poisson clock + // 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] = cumulative_times[src_idx]; + simulation::decay_times[i - global_start] = decay_time; } } } @@ -647,11 +666,13 @@ 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 activity-based timing is on - if (settings::activity_based_timing) { + // 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]; - p.time() = decay_time; - p.time_last() = decay_time; + if (!std::isnan(decay_time)) { + p.time() = decay_time; + p.time_last() = decay_time; + } } } p.current_work() = index_source; diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index b94a8d79d3b..0298ebaec6d 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -125,20 +125,44 @@ def test_source_dlopen(): assert 'library' in elem.attrib -def test_activity_based_timing_setting_roundtrip(): - """Test XML round-trip for activity_based_timing setting.""" - settings = openmc.Settings() - settings.activity_based_timing = True - settings.particles = 100 - settings.batches = 1 - settings.run_mode = 'fixed source' - elem = settings.to_xml_element() - text = elem.find('activity_based_timing') - assert text is not None - assert text.text == 'true' - - new_settings = openmc.Settings.from_xml_element(elem) - assert new_settings.activity_based_timing is True +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(): From dd362849ecde7f3fa6113e454f0426b674c8a53b Mon Sep 17 00:00:00 2001 From: Christopher Date: Thu, 12 Feb 2026 15:10:41 +0100 Subject: [PATCH 8/9] clang --- include/openmc/settings.h | 2 +- src/settings.cpp | 1 - src/simulation.cpp | 2 +- 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 8d0503190c8..d58b0a09c03 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -103,7 +103,7 @@ extern bool weight_window_checkpoint_surface; //!< enable weight window check extern bool weight_window_checkpoint_collision; //!< enable weight window check //!< upon collision? extern bool write_all_tracks; //!< write track files for every particle? -extern bool write_initial_source; //!< write out initial source file? +extern bool write_initial_source; //!< write out initial source file? // Paths to various files extern std::string path_cross_sections; //!< path to cross_sections.xml diff --git a/src/settings.cpp b/src/settings.cpp index fb17493241d..5b472468fc4 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -88,7 +88,6 @@ bool weight_window_checkpoint_collision {true}; bool write_all_tracks {false}; bool write_initial_source {false}; - std::string path_cross_sections; std::string path_input; std::string path_output; diff --git a/src/simulation.cpp b/src/simulation.cpp index ebc94edaaea..f344fa553b6 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -2,9 +2,9 @@ #include "openmc/bank.h" #include "openmc/capi.h" -#include "openmc/distribution.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" From 1fad4a96b9d29997c289ae08c6602c3f6d5afebd Mon Sep 17 00:00:00 2001 From: Christopher Date: Thu, 12 Feb 2026 16:53:19 +0100 Subject: [PATCH 9/9] clang again --- include/openmc/settings.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/openmc/settings.h b/include/openmc/settings.h index d58b0a09c03..b369c99fef8 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -102,8 +102,8 @@ extern bool weight_window_checkpoint_surface; //!< enable weight window check //!< upon surface crossing? extern bool weight_window_checkpoint_collision; //!< enable weight window check //!< upon collision? -extern bool write_all_tracks; //!< write track files for every particle? -extern bool write_initial_source; //!< write out initial source file? +extern bool write_all_tracks; //!< write track files for every particle? +extern bool write_initial_source; //!< write out initial source file? // Paths to various files extern std::string path_cross_sections; //!< path to cross_sections.xml