Shared secondary bank performance optimizations (#3995)

include/openmc/particle_data.h

@@ -582,10 +582,8 @@ public:
   // 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

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;

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<uint64_t>(id) * prn_stride);
   for (int i = 0; i < N_STREAMS; i++) {
-    seeds[i] =
-      future_seed(static_cast<uint64_t>(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;
 }
 
 //==============================================================================

src/simulation.cpp

@@ -886,11 +886,14 @@ void transport_history_based_single_particle(Particle& p)
 
 void transport_history_based()
 {
-#pragma omp parallel for schedule(runtime)
-  for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
+#pragma omp parallel
+  {
     Particle p;
-    initialize_particle_track(p, i_work, false);
-    transport_history_based_single_particle(p);
+#pragma omp for schedule(runtime)
+    for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
+      initialize_particle_track(p, i_work, false);
+      transport_history_based_single_particle(p);
+    }
   }
 }
 
@@ -923,15 +926,16 @@ void transport_history_based_shared_secondary()
 #pragma omp parallel
   {
     vector<SourceSite> 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 +992,11 @@ void transport_history_based_shared_secondary()
 #pragma omp parallel
     {
       vector<SourceSite> 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 +1004,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