From 2946473c3c44ee1d82404ccd9f7100fed6232779 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 24 Jun 2026 16:50:48 -0500 Subject: [PATCH 1/8] Only use single cache in ParticleData --- include/openmc/particle_data.h | 23 +++++++++++++++++- src/particle.cpp | 1 + src/particle_data.cpp | 44 ++++++++++++++++++++++++++++++---- src/particle_restart.cpp | 1 + 4 files changed, 64 insertions(+), 5 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index f72948f6eb4..b68a9569090 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -1,6 +1,8 @@ #ifndef OPENMC_PARTICLE_DATA_H #define OPENMC_PARTICLE_DATA_H +#include + #include "openmc/array.h" #include "openmc/constants.h" #include "openmc/particle_type.h" @@ -584,12 +586,24 @@ class ParticleData : public GeometryState { // Cross section caches NuclideMicroXS& neutron_xs(int i) { + ensure_neutron_xs_cache(); return neutron_xs_[i]; } // Microscopic neutron cross sections const NuclideMicroXS& neutron_xs(int i) const { return neutron_xs_[i]; } // Microscopic photon cross sections - ElementMicroXS& photon_xs(int i) { return photon_xs_[i]; } + ElementMicroXS& photon_xs(int i) + { + ensure_photon_xs_cache(); + return photon_xs_[i]; + } + const ElementMicroXS& photon_xs(int i) const { return photon_xs_[i]; } + + void ensure_xs_cache_for_type(); + void ensure_neutron_xs_cache(); + void ensure_photon_xs_cache(); + size_t neutron_xs_cache_size() const { return neutron_xs_.size(); } + size_t photon_xs_cache_size() const { return photon_xs_.size(); } // Macroscopic cross sections MacroXS& macro_xs() { return macro_xs_; } @@ -774,6 +788,9 @@ class ParticleData : public GeometryState { //! Force recalculation of neutron xs by setting last energy to zero void invalidate_neutron_xs() { + if (!type().is_neutron()) + return; + ensure_neutron_xs_cache(); for (auto& micro : neutron_xs_) micro.last_E = 0.0; } @@ -794,6 +811,10 @@ class ParticleData : public GeometryState { d = 0; } } + +private: + void release_neutron_xs_cache(); + void release_photon_xs_cache(); }; } // namespace openmc diff --git a/src/particle.cpp b/src/particle.cpp index 65ff2ad6192..ef788c1ef98 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -160,6 +160,7 @@ void Particle::from_source(const SourceSite* src) // Copy attributes from source bank site type() = src->particle; + ensure_xs_cache_for_type(); wgt() = src->wgt; wgt_last() = src->wgt; r() = src->r; diff --git a/src/particle_data.cpp b/src/particle_data.cpp index 370ca12e469..8f4849356d5 100644 --- a/src/particle_data.cpp +++ b/src/particle_data.cpp @@ -105,16 +105,52 @@ ParticleData::ParticleData() // Allocate space for tally filter matches filter_matches_.resize(model::tally_filters.size()); - // Create microscopic cross section caches - neutron_xs_.resize(data::nuclides.size()); - photon_xs_.resize(data::elements.size()); - // Creates the pulse-height storage for the particle if (!model::pulse_height_cells.empty()) { pht_storage_.resize(model::pulse_height_cells.size(), 0.0); } } +void ParticleData::ensure_xs_cache_for_type() +{ + if (type().is_neutron()) { + ensure_neutron_xs_cache(); + release_photon_xs_cache(); + } else if (type().is_photon()) { + ensure_photon_xs_cache(); + release_neutron_xs_cache(); + } else { + release_neutron_xs_cache(); + release_photon_xs_cache(); + } +} + +void ParticleData::ensure_neutron_xs_cache() +{ + if (neutron_xs_.size() != data::nuclides.size()) { + neutron_xs_.resize(data::nuclides.size()); + } +} + +void ParticleData::ensure_photon_xs_cache() +{ + if (photon_xs_.size() != data::elements.size()) { + photon_xs_.resize(data::elements.size()); + } +} + +void ParticleData::release_neutron_xs_cache() +{ + vector empty; + neutron_xs_.swap(empty); +} + +void ParticleData::release_photon_xs_cache() +{ + vector empty; + photon_xs_.swap(empty); +} + TrackState ParticleData::get_track_state() const { TrackState state; diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index c226d51ec2a..5b12eb54056 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -58,6 +58,7 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode) read_dataset(file_id, "type", type); p.type() = legacy_particle_codes ? legacy_particle_index_to_type(type) : ParticleType {type}; + p.ensure_xs_cache_for_type(); read_dataset(file_id, "weight", p.wgt()); read_dataset(file_id, "energy", p.E()); read_dataset(file_id, "xyz", p.r()); From d2bf6b14d0827214c758872393d61beed5caa2ac Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 24 Jun 2026 16:52:05 -0500 Subject: [PATCH 2/8] Don't use ensure_*_cache in neutron_xs / photon_xs --- include/openmc/particle_data.h | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index b68a9569090..26ffea7c6d6 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -584,19 +584,12 @@ class ParticleData : public GeometryState { // Methods and accessors // Cross section caches - NuclideMicroXS& neutron_xs(int i) - { - ensure_neutron_xs_cache(); - return neutron_xs_[i]; - } // Microscopic neutron cross sections + // Microscopic neutron cross sections + NuclideMicroXS& neutron_xs(int i) { return neutron_xs_[i]; } const NuclideMicroXS& neutron_xs(int i) const { return neutron_xs_[i]; } // Microscopic photon cross sections - ElementMicroXS& photon_xs(int i) - { - ensure_photon_xs_cache(); - return photon_xs_[i]; - } + ElementMicroXS& photon_xs(int i) { return photon_xs_[i]; } const ElementMicroXS& photon_xs(int i) const { return photon_xs_[i]; } void ensure_xs_cache_for_type(); From 6b435bccdca9ff126ca68847e71561076e2e9494 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 24 Jun 2026 18:12:42 -0500 Subject: [PATCH 3/8] Move neutron cache invalidation --- src/particle.cpp | 7 +++++++ src/simulation.cpp | 5 ----- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index ef788c1ef98..953366f23e0 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -192,6 +192,13 @@ void Particle::from_source(const SourceSite* src) wgt_born() = src->wgt_born; wgt_ww_born() = src->wgt_ww_born; n_split() = src->n_split; + + // Force CE neutron cross sections to be recalculated after source state is + // known. For shared secondary transport, this avoids allocating a neutron + // cache before a secondary's actual particle type has been applied. + if (settings::run_CE) { + invalidate_neutron_xs(); + } } void Particle::event_calculate_xs() diff --git a/src/simulation.cpp b/src/simulation.cpp index 900a383ab33..923803d50de 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -681,11 +681,6 @@ void initialize_particle_track( simulation::total_weight += p.wgt(); } - // Force calculation of cross-sections by setting last energy to zero - if (settings::run_CE) { - p.invalidate_neutron_xs(); - } - // Prepare to write out particle track. if (p.write_track()) add_particle_track(p); From a9992c4086183d3c065f0f4560af15e40f62b3aa Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 25 Jun 2026 17:44:47 -0500 Subject: [PATCH 4/8] Factor out common skip-ahead in init_particle_seeds --- src/random_lcg.cpp | 72 ++++++++++++++++++++++++++++------------------ 1 file changed, 44 insertions(+), 28 deletions(-) diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp index 29457569b94..f2a81fc1c7e 100644 --- a/src/random_lcg.cpp +++ b/src/random_lcg.cpp @@ -12,6 +12,45 @@ constexpr uint64_t prn_mult {6364136223846793005ULL}; // multiplication constexpr uint64_t prn_add {1442695040888963407ULL}; // additive factor, c uint64_t prn_stride {DEFAULT_STRIDE}; // stride between particles +namespace { + +struct SkipAheadCoefficients { + uint64_t multiplier; + uint64_t increment; +}; + +SkipAheadCoefficients future_seed_coefficients(uint64_t n) +{ + // The algorithm here to determine the parameters used to skip ahead is + // described in F. Brown, "Random Number Generation with Arbitrary Stride," + // Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in + // O(log2(N)) operations instead of O(N). Basically, it computes parameters G + // and C which can then be used to find x_N = G*x_0 + C mod 2^M. + + // Initialize constants + uint64_t g {prn_mult}; + uint64_t c {prn_add}; + uint64_t g_new {1}; + uint64_t c_new {0}; + + while (n > 0) { + // Check if the least significant bit is 1. + if (n & 1) { + g_new *= g; + c_new = c_new * g + c; + } + c *= (g + 1); + g *= g; + + // Move bits right, dropping least significant bit. + n >>= 1; + } + + return {g_new, c_new}; +} + +} // namespace + //============================================================================== // PRN //============================================================================== @@ -69,9 +108,10 @@ uint64_t init_seed(int64_t id, int offset) void init_particle_seeds(int64_t id, uint64_t* seeds) { + auto [multiplier, increment] = + future_seed_coefficients(static_cast(id) * prn_stride); for (int i = 0; i < N_STREAMS; i++) { - seeds[i] = - future_seed(static_cast(id) * prn_stride, master_seed + i); + seeds[i] = multiplier * (master_seed + i) + increment; } } @@ -90,33 +130,9 @@ void advance_prn_seed(int64_t n, uint64_t* seed) uint64_t future_seed(uint64_t n, uint64_t seed) { - // The algorithm here to determine the parameters used to skip ahead is - // described in F. Brown, "Random Number Generation with Arbitrary Stride," - // Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in - // O(log2(N)) operations instead of O(N). Basically, it computes parameters G - // and C which can then be used to find x_N = G*x_0 + C mod 2^M. - - // Initialize constants - uint64_t g {prn_mult}; - uint64_t c {prn_add}; - uint64_t g_new {1}; - uint64_t c_new {0}; - - while (n > 0) { - // Check if the least significant bit is 1. - if (n & 1) { - g_new *= g; - c_new = c_new * g + c; - } - c *= (g + 1); - g *= g; - - // Move bits right, dropping least significant bit. - n >>= 1; - } - // With G and C, we can now find the new seed. - return g_new * seed + c_new; + auto [multiplier, increment] = future_seed_coefficients(n); + return multiplier * seed + increment; } //============================================================================== From 966433be50a8e9df5cd1edf22f2521f72f9fe64e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 25 Jun 2026 17:57:19 -0500 Subject: [PATCH 5/8] Avoid atomics in event_death as much as possible --- src/particle.cpp | 32 ++++++++++++++++++++++++-------- 1 file changed, 24 insertions(+), 8 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 953366f23e0..8dd36c9d9b3 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -562,15 +562,30 @@ void Particle::event_death() finalize_particle_track(*this); } -// Contribute tally reduction variables to global accumulator + // Contribute tally reduction variables to global accumulator + const auto k_absorption = keff_tally_absorption(); + const auto k_collision = keff_tally_collision(); + const auto k_tracklength = keff_tally_tracklength(); + const auto leakage = keff_tally_leakage(); + + if (settings::run_mode == RunMode::EIGENVALUE) { + if (k_absorption != 0.0) { #pragma omp atomic - global_tally_absorption += keff_tally_absorption(); + global_tally_absorption += k_absorption; + } + if (k_collision != 0.0) { #pragma omp atomic - global_tally_collision += keff_tally_collision(); + global_tally_collision += k_collision; + } + if (k_tracklength != 0.0) { #pragma omp atomic - global_tally_tracklength += keff_tally_tracklength(); + global_tally_tracklength += k_tracklength; + } + } + if (leakage != 0.0) { #pragma omp atomic - global_tally_leakage += keff_tally_leakage(); + global_tally_leakage += leakage; + } // Reset particle tallies once accumulated keff_tally_absorption() = 0.0; @@ -590,9 +605,10 @@ void Particle::event_death() // Record the number of progeny created by this particle. // This data will be used to efficiently sort the fission bank. - if (settings::run_mode == RunMode::EIGENVALUE || - settings::use_shared_secondary_bank) { - simulation::progeny_per_particle[current_work()] = n_progeny(); + const auto progeny = n_progeny(); + if (progeny != 0 && (settings::run_mode == RunMode::EIGENVALUE || + settings::use_shared_secondary_bank)) { + simulation::progeny_per_particle[current_work()] = progeny; } } From 13ddbd36795b5bf30b0b4766ba0cfd2d21484886 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 27 Jun 2026 17:21:00 -0500 Subject: [PATCH 6/8] Reuse Particle for each OpenMP thread --- src/simulation.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/simulation.cpp b/src/simulation.cpp index 923803d50de..e64c226d00c 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -918,15 +918,16 @@ void transport_history_based_shared_secondary() #pragma omp parallel { vector thread_bank; + Particle p; #pragma omp for schedule(runtime) for (int64_t i = 1; i <= simulation::work_per_rank; i++) { - Particle p; initialize_particle_track(p, i, false); transport_history_based_single_particle(p); for (auto& site : p.local_secondary_bank()) { thread_bank.push_back(site); } + p.local_secondary_bank().clear(); } // Drain thread-local bank into the shared secondary bank (once per thread) @@ -983,11 +984,11 @@ void transport_history_based_shared_secondary() #pragma omp parallel { vector thread_bank; + Particle p; #pragma omp for schedule(runtime) for (int64_t i = 1; i <= simulation::shared_secondary_bank_read.size(); i++) { - Particle p; initialize_particle_track(p, i, true); SourceSite& site = simulation::shared_secondary_bank_read[i - 1]; p.event_revive_from_secondary(site); @@ -995,6 +996,7 @@ void transport_history_based_shared_secondary() for (auto& secondary_site : p.local_secondary_bank()) { thread_bank.push_back(secondary_site); } + p.local_secondary_bank().clear(); } // Drain thread-local bank into the shared secondary bank (once per From 2c18425ca05b30ed2a4435be5aff929df8b8ba49 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 28 Jun 2026 12:43:07 -0500 Subject: [PATCH 7/8] Revert xs cache changes (unnecessary when Particle is reused) --- include/openmc/particle_data.h | 16 ------------- src/particle.cpp | 8 ------- src/particle_data.cpp | 44 ++++------------------------------ src/particle_restart.cpp | 1 - src/simulation.cpp | 5 ++++ 5 files changed, 9 insertions(+), 65 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 26ffea7c6d6..44e82fd2350 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -1,8 +1,6 @@ #ifndef OPENMC_PARTICLE_DATA_H #define OPENMC_PARTICLE_DATA_H -#include - #include "openmc/array.h" #include "openmc/constants.h" #include "openmc/particle_type.h" @@ -590,13 +588,6 @@ class ParticleData : public GeometryState { // Microscopic photon cross sections ElementMicroXS& photon_xs(int i) { return photon_xs_[i]; } - const ElementMicroXS& photon_xs(int i) const { return photon_xs_[i]; } - - void ensure_xs_cache_for_type(); - void ensure_neutron_xs_cache(); - void ensure_photon_xs_cache(); - size_t neutron_xs_cache_size() const { return neutron_xs_.size(); } - size_t photon_xs_cache_size() const { return photon_xs_.size(); } // Macroscopic cross sections MacroXS& macro_xs() { return macro_xs_; } @@ -781,9 +772,6 @@ class ParticleData : public GeometryState { //! Force recalculation of neutron xs by setting last energy to zero void invalidate_neutron_xs() { - if (!type().is_neutron()) - return; - ensure_neutron_xs_cache(); for (auto& micro : neutron_xs_) micro.last_E = 0.0; } @@ -804,10 +792,6 @@ class ParticleData : public GeometryState { d = 0; } } - -private: - void release_neutron_xs_cache(); - void release_photon_xs_cache(); }; } // namespace openmc diff --git a/src/particle.cpp b/src/particle.cpp index 8dd36c9d9b3..412710e0b21 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -160,7 +160,6 @@ void Particle::from_source(const SourceSite* src) // Copy attributes from source bank site type() = src->particle; - ensure_xs_cache_for_type(); wgt() = src->wgt; wgt_last() = src->wgt; r() = src->r; @@ -192,13 +191,6 @@ void Particle::from_source(const SourceSite* src) wgt_born() = src->wgt_born; wgt_ww_born() = src->wgt_ww_born; n_split() = src->n_split; - - // Force CE neutron cross sections to be recalculated after source state is - // known. For shared secondary transport, this avoids allocating a neutron - // cache before a secondary's actual particle type has been applied. - if (settings::run_CE) { - invalidate_neutron_xs(); - } } void Particle::event_calculate_xs() diff --git a/src/particle_data.cpp b/src/particle_data.cpp index 8f4849356d5..370ca12e469 100644 --- a/src/particle_data.cpp +++ b/src/particle_data.cpp @@ -105,52 +105,16 @@ ParticleData::ParticleData() // Allocate space for tally filter matches filter_matches_.resize(model::tally_filters.size()); + // Create microscopic cross section caches + neutron_xs_.resize(data::nuclides.size()); + photon_xs_.resize(data::elements.size()); + // Creates the pulse-height storage for the particle if (!model::pulse_height_cells.empty()) { pht_storage_.resize(model::pulse_height_cells.size(), 0.0); } } -void ParticleData::ensure_xs_cache_for_type() -{ - if (type().is_neutron()) { - ensure_neutron_xs_cache(); - release_photon_xs_cache(); - } else if (type().is_photon()) { - ensure_photon_xs_cache(); - release_neutron_xs_cache(); - } else { - release_neutron_xs_cache(); - release_photon_xs_cache(); - } -} - -void ParticleData::ensure_neutron_xs_cache() -{ - if (neutron_xs_.size() != data::nuclides.size()) { - neutron_xs_.resize(data::nuclides.size()); - } -} - -void ParticleData::ensure_photon_xs_cache() -{ - if (photon_xs_.size() != data::elements.size()) { - photon_xs_.resize(data::elements.size()); - } -} - -void ParticleData::release_neutron_xs_cache() -{ - vector empty; - neutron_xs_.swap(empty); -} - -void ParticleData::release_photon_xs_cache() -{ - vector empty; - photon_xs_.swap(empty); -} - TrackState ParticleData::get_track_state() const { TrackState state; diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index 5b12eb54056..c226d51ec2a 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -58,7 +58,6 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode) read_dataset(file_id, "type", type); p.type() = legacy_particle_codes ? legacy_particle_index_to_type(type) : ParticleType {type}; - p.ensure_xs_cache_for_type(); read_dataset(file_id, "weight", p.wgt()); read_dataset(file_id, "energy", p.E()); read_dataset(file_id, "xyz", p.r()); diff --git a/src/simulation.cpp b/src/simulation.cpp index e64c226d00c..7baa387929a 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -681,6 +681,11 @@ void initialize_particle_track( simulation::total_weight += p.wgt(); } + // Force calculation of cross-sections by setting last energy to zero + if (settings::run_CE) { + p.invalidate_neutron_xs(); + } + // Prepare to write out particle track. if (p.write_track()) add_particle_track(p); From 9457ab4e73627dd822435f8cb36b2992839eaeae Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 3 Jul 2026 16:58:40 -0500 Subject: [PATCH 8/8] Remove n_progeny change --- src/particle.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 412710e0b21..10e8f213dd4 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -597,10 +597,9 @@ void Particle::event_death() // Record the number of progeny created by this particle. // This data will be used to efficiently sort the fission bank. - const auto progeny = n_progeny(); - if (progeny != 0 && (settings::run_mode == RunMode::EIGENVALUE || - settings::use_shared_secondary_bank)) { - simulation::progeny_per_particle[current_work()] = progeny; + if (settings::run_mode == RunMode::EIGENVALUE || + settings::use_shared_secondary_bank) { + simulation::progeny_per_particle[current_work()] = n_progeny(); } }