diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index f72948f6eb4..44e82fd2350 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -582,10 +582,8 @@ class ParticleData : public GeometryState { // Methods and accessors // Cross section caches - NuclideMicroXS& neutron_xs(int i) - { - 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 diff --git a/src/particle.cpp b/src/particle.cpp index 65ff2ad6192..10e8f213dd4 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -554,15 +554,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; 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; } //============================================================================== diff --git a/src/simulation.cpp b/src/simulation.cpp index 900a383ab33..7baa387929a 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -923,15 +923,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) @@ -988,11 +989,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); @@ -1000,6 +1001,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