Shared Secondary Particle Bank (#3863)

Co-authored-by: John Tramm <jtramm@gmail.com>
Co-authored-by: Claude Opus 4.6 <noreply@anthropic.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
Co-authored-by: Copilot <copilot@github.com>

include/openmc/bank.h

@@ -34,18 +34,24 @@ extern vector<vector<double>> ifp_fission_lifetime_bank;
 
 extern vector<int64_t> progeny_per_particle;
 
+extern SharedArray<SourceSite> shared_secondary_bank_read;
+extern SharedArray<SourceSite> shared_secondary_bank_write;
+
 } // namespace simulation
 
 //==============================================================================
 // Non-member functions
 //==============================================================================
 
-void sort_fission_bank();
+void sort_bank(SharedArray<SourceSite>& bank, bool is_fission_bank);
 
 void free_memory_bank();
 
 void init_fission_bank(int64_t max);
 
+int64_t synchronize_global_secondary_bank(
+  SharedArray<SourceSite>& shared_secondary_bank);
+
 } // namespace openmc
 
 #endif // OPENMC_BANK_H

include/openmc/event.h

@@ -112,6 +112,19 @@ void process_collision_events();
 //! \param n_particles The number of particles in the particle buffer
 void process_death_events(int64_t n_particles);
 
+//! Process event queues until all are empty. Each iteration processes the
+//! longest queue first to maximize vectorization efficiency.
+void process_transport_events();
+
+//! Initialize secondary particles from a shared secondary bank for
+//! event-based transport
+//
+//! \param n_particles The number of particles to initialize
+//! \param offset The offset index in the shared secondary bank
+//! \param shared_secondary_bank The shared secondary bank to read from
+void process_init_secondary_events(int64_t n_particles, int64_t offset,
+  const SharedArray<SourceSite>& shared_secondary_bank);
+
 } // namespace openmc
 
 #endif // OPENMC_EVENT_H

include/openmc/particle.h

@@ -71,7 +71,8 @@ public:
   void event_advance();
   void event_cross_surface();
   void event_collide();
-  void event_revive_from_secondary();
+  void event_revive_from_secondary(const SourceSite& site);
+  void event_check_limit_and_revive();
   void event_death();
 
   //! pulse-height recording

include/openmc/particle_data.h

@@ -50,8 +50,11 @@ struct SourceSite {
 
   // Extra attributes that don't show up in source written to file
   int parent_nuclide {-1};
-  int64_t parent_id;
-  int64_t progeny_id;
+  int64_t parent_id {0};
+  int64_t progeny_id {0};
+  double wgt_born {1.0};
+  double wgt_ww_born {-1.0};
+  int64_t n_split {0};
 };
 
 struct CollisionTrackSite {
@@ -533,14 +536,14 @@ private:
   uint64_t seeds_[N_STREAMS];
   int stream_;
 
-  vector<SourceSite> secondary_bank_;
+  vector<SourceSite> local_secondary_bank_;
 
   // Keep track of how many secondary particles were created in the collision
   // and what the starting index is in the secondary bank for this particle
   int n_secondaries_ {0};
   int secondary_bank_index_ {0};
 
-  int64_t current_work_;
+  int64_t current_work_ {0};
 
   vector<double> flux_derivs_;
 
@@ -563,7 +566,9 @@ private:
 
   int n_event_ {0};
 
-  int n_split_ {0};
+  int64_t n_tracks_ {0}; //!< number of tracks in this particle history
+
+  int64_t n_split_ {0};
   double ww_factor_ {0.0};
 
   int64_t n_progeny_ {0};
@@ -693,12 +698,14 @@ public:
   int& stream() { return stream_; }
 
   // secondary particle bank
-  SourceSite& secondary_bank(int i) { return secondary_bank_[i]; }
-  const SourceSite& secondary_bank(int i) const { return secondary_bank_[i]; }
-  decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; }
-  decltype(secondary_bank_) const& secondary_bank() const
+  SourceSite& local_secondary_bank(int i) { return local_secondary_bank_[i]; }
+  const SourceSite& local_secondary_bank(int i) const
   {
-    return secondary_bank_;
+    return local_secondary_bank_[i];
+  }
+  decltype(local_secondary_bank_)& local_secondary_bank()
+  {
+    return local_secondary_bank_;
   }
 
   // Number of secondaries created in a collision
@@ -747,13 +754,16 @@ public:
   int& n_event() { return n_event_; }
 
   // Number of times variance reduction has caused a particle split
-  int n_split() const { return n_split_; }
-  int& n_split() { return n_split_; }
+  int64_t n_split() const { return n_split_; }
+  int64_t& n_split() { return n_split_; }
 
   // Particle-specific factor for on-the-fly weight window adjustment
   double ww_factor() const { return ww_factor_; }
   double& ww_factor() { return ww_factor_; }
 
+  // Number of tracks in this particle history
+  int64_t& n_tracks() { return n_tracks_; }
+
   // Number of progeny produced by this particle
   int64_t& n_progeny() { return n_progeny_; }
 

include/openmc/settings.h

@@ -98,7 +98,8 @@ extern bool uniform_source_sampling; //!< sample sources uniformly?
 extern bool ufs_on;                  //!< uniform fission site method on?
 extern bool urr_ptables_on;          //!< use unresolved resonance prob. tables?
 extern bool use_decay_photons;       //!< use decay photons for D1S
-extern "C" bool weight_windows_on;   //!< are weight windows are enabled?
+extern bool use_shared_secondary_bank; //!< Use shared bank for secondaries
+extern "C" bool weight_windows_on;     //!< are weight windows are enabled?
 extern bool weight_window_checkpoint_surface;   //!< enable weight window check
                                                 //!< upon surface crossing?
 extern bool weight_window_checkpoint_collision; //!< enable weight window check

include/openmc/shared_array.h

@@ -4,6 +4,8 @@
 //! \file shared_array.h
 //! \brief Shared array data structure
 
+#include <algorithm> // for copy_n
+
 #include "openmc/memory.h"
 
 namespace openmc {
@@ -30,14 +32,12 @@ public:
   //! Default constructor.
   SharedArray() = default;
 
-  //! Construct a zero size container with space to hold capacity number of
-  //! elements.
+  //! Construct a container with `size` elements and capacity equal to `size`.
   //
-  //! \param capacity The number of elements for the container to allocate
-  //! space for
-  SharedArray(int64_t capacity) : capacity_(capacity)
+  //! \param size The number of elements to allocate and initialize
+  SharedArray(int64_t size) : size_(size), capacity_(size)
   {
-    data_ = make_unique<T[]>(capacity);
+    data_ = make_unique<T[]>(size);
   }
 
   //==========================================================================
@@ -97,8 +97,26 @@ public:
     capacity_ = 0;
   }
 
+  //! Push back an element to the array, with capacity and reallocation behavior
+  //! as if this were a vector. This does not perform any thread safety checks.
+  //! If the size exceeds the capacity, then the capacity will double just as
+  //! with a vector. Data will be reallocated and moved to a new pointer and
+  //! copied in before the new item is appended. Old data will be freed.
+  void thread_unsafe_append(const T& value)
+  {
+    if (size_ == capacity_) {
+      int64_t new_capacity = capacity_ == 0 ? 8 : 2 * capacity_;
+      unique_ptr<T[]> new_data = make_unique<T[]>(new_capacity);
+      std::copy_n(data_.get(), size_, new_data.get());
+      data_ = std::move(new_data);
+      capacity_ = new_capacity;
+    }
+    data_[size_++] = value;
+  }
+
   //! Return the number of elements in the container
   int64_t size() { return size_; }
+  int64_t size() const { return size_; }
 
   //! Resize the container to contain a specified number of elements. This is
   //! useful in cases where the container is written to in a non-thread safe

include/openmc/simulation.h

@@ -49,6 +49,9 @@ extern const RegularMesh* ufs_mesh;
 extern vector<double> k_generation;
 extern vector<int64_t> work_index;
 
+extern int64_t
+  simulation_tracks_completed; //!< Number of tracks completed on this rank
+
 } // namespace simulation
 
 //==============================================================================
@@ -59,7 +62,7 @@ extern vector<int64_t> work_index;
 void allocate_banks();
 
 //! Determine number of particles to transport per process
-void calculate_work();
+void calculate_work(int64_t n_particles);
 
 //! Initialize nuclear data before a simulation
 void initialize_data();
@@ -70,8 +73,9 @@ void initialize_batch();
 //! Initialize a fission generation
 void initialize_generation();
 
-//! Full initialization of a particle history
-void initialize_history(Particle& p, int64_t index_source);
+//! Full initialization of a particle track
+void initialize_particle_track(
+  Particle& p, int64_t index_source, bool is_secondary);
 
 //! Finalize a batch
 //!
@@ -92,16 +96,35 @@ void broadcast_results();
 
 void free_memory_simulation();
 
-//! Simulate a single particle history (and all generated secondary particles,
-//!  if enabled), from birth to death
+//! Compute unique particle ID from a 1-based source index
+//! \param index_source 1-based source index within this rank's work
+//! \return globally unique particle ID
+int64_t compute_particle_id(int64_t index_source);
+
+//! Compute the transport RNG seed from a particle ID
+//! \param particle_id the particle's globally unique ID
+//! \return seed value passed to init_particle_seeds()
+int64_t compute_transport_seed(int64_t particle_id);
+
+//! Simulate a single particle history from birth to death, inclusive of any
+//! secondary particles. In shared secondary mode, only a single track is
+//! transported and secondaries are deposited into a shared bank instead.
 void transport_history_based_single_particle(Particle& p);
 
 //! Simulate all particle histories using history-based parallelism
 void transport_history_based();
 
+//! Simulate all particles using history-based parallelism, with a shared
+//! secondary bank
+void transport_history_based_shared_secondary();
+
 //! Simulate all particle histories using event-based parallelism
 void transport_event_based();
 
+//! Simulate all particles using event-based parallelism, with a shared
+//! secondary bank
+void transport_event_based_shared_secondary();
+
 } // namespace openmc
 
 #endif // OPENMC_SIMULATION_H