Shared Secondary Particle Bank

docs/source/io_formats/settings.rst

@@ -1711,6 +1711,17 @@ following sub-elements/attributes:
   The ``weight_windows_file`` element has no attributes and contains the path to
   a weight windows HDF5 file to load during simulation initialization.
 
+----------------------------------------
+``<shared_secondary_bank>`` Element
+----------------------------------------
+
+  The ``shared_secondary_bank`` element indicates whether to use a shared
+  secondary particle bank. When enabled, secondary particles are collected into
+  a global bank, sorted for reproducibility, and load-balanced across MPI ranks
+  between generations. If not specified, the shared secondary bank is enabled
+  automatically for fixed-source simulations with weight windows active, and
+  disabled otherwise.
+
 -------------------------------
 ``<weight_windows_on>`` Element
 -------------------------------

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,
+  SharedArray<SourceSite>& shared_secondary_bank);
+
 } // namespace openmc
 
 #endif // OPENMC_EVENT_H

include/openmc/particle.h

@@ -69,7 +69,8 @@ class Particle : public ParticleData {
   void event_advance();
   void event_cross_surface();
   void event_collide();
-  void event_revive_from_secondary();
+  void event_revive_from_secondary(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 @@ class ParticleData : public GeometryState {
   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 @@ class ParticleData : public GeometryState {
 
   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 @@ class ParticleData : public GeometryState {
   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 @@ class ParticleData : public GeometryState {
   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

@@ -200,6 +200,8 @@ extern int trigger_batch_interval; //!< Batch interval for triggers
 extern "C" int verbosity;          //!< How verbose to make output
 extern double weight_cutoff;       //!< Weight cutoff for Russian roulette
 extern double weight_survive;      //!< Survival weight after Russian roulette
+extern bool
+  use_shared_secondary_bank; //!< Use shared bank for secondary particles
 
 } // namespace settings
 

include/openmc/shared_array.h

@@ -30,14 +30,13 @@ class SharedArray {
   //! Default constructor.
   SharedArray() = default;
 
-  //! Construct a zero size container with space to hold capacity number of
-  //! elements.
+  //! Construct a container with size elements, with space to hold size number
+  //! of elements.
   //
-  //! \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 +96,28 @@ class SharedArray {
     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 ? 1 : 2 * capacity_;
+      unique_ptr<T[]> new_data = make_unique<T[]>(new_capacity);
+      for (int64_t i = 0; i < size_; i++) {
+        new_data[i] = std::move(data_[i]);
+      }
+      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();
@@ -71,7 +74,7 @@ void initialize_batch();
 void initialize_generation();
 
 //! Full initialization of a particle history
-void initialize_history(Particle& p, int64_t index_source);
+void initialize_history(Particle& p, int64_t index_source, bool is_secondary);
 
 //! Finalize a batch
 //!
@@ -92,16 +95,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

openmc/lib/core.py

@@ -31,7 +31,10 @@ class _SourceSite(Structure):
                 ('particle', c_int32),
                 ('parent_nuclide', c_int),
                 ('parent_id', c_int64),
-                ('progeny_id', c_int64)]
+                ('progeny_id', c_int64),
+                ('wgt_born', c_double),
+                ('wgt_ww_born', c_double),
+                ('n_split', c_int64)]
 
 # Define input type for numpy arrays that will be passed into C++ functions
 # Must be an int or double array, with single dimension that is contiguous

openmc/settings.py

@@ -371,6 +371,15 @@ class Settings:
         Whether delayed neutrons are created in fission.
 
         .. versionadded:: 0.13.3
+    shared_secondary_bank : bool
+        Whether to use a shared secondary particle bank. When enabled,
+        secondary particles are collected into a global bank, sorted for
+        reproducibility, and load-balanced across MPI ranks between
+        generations. If not specified, the shared secondary bank is
+        enabled automatically for fixed-source simulations with weight
+        windows active, and disabled otherwise.
+
+        .. versionadded:: 0.15.1
     weight_windows_on : bool
         Whether weight windows are enabled
 
@@ -475,6 +484,7 @@ def __init__(self, **kwargs):
         self._weight_window_generators = cv.CheckedList(
             WeightWindowGenerator, 'weight window generators')
         self._weight_windows_on = None
+        self._shared_secondary_bank = None
         self._weight_windows_file = None
         self._weight_window_checkpoints = {}
         self._max_history_splits = None
@@ -1287,6 +1297,15 @@ def weight_windows_on(self, value: bool):
         cv.check_type('weight windows on', value, bool)
         self._weight_windows_on = value
 
+    @property
+    def shared_secondary_bank(self) -> bool:
+        return self._shared_secondary_bank
+
+    @shared_secondary_bank.setter
+    def shared_secondary_bank(self, value: bool):
+        cv.check_type('shared secondary bank', value, bool)
+        self._shared_secondary_bank = value
+
     @property
     def weight_window_checkpoints(self) -> dict:
         return self._weight_window_checkpoints
@@ -1891,6 +1910,11 @@ def _create_weight_windows_on_subelement(self, root):
             elem = ET.SubElement(root, "weight_windows_on")
             elem.text = str(self._weight_windows_on).lower()
 
+    def _create_shared_secondary_bank_subelement(self, root):
+        if self._shared_secondary_bank is not None:
+            elem = ET.SubElement(root, "shared_secondary_bank")
+            elem.text = str(self._shared_secondary_bank).lower()
+
     def _create_weight_window_generators_subelement(self, root, mesh_memo=None):
         if not self.weight_window_generators:
             return
@@ -2379,6 +2403,11 @@ def _weight_windows_on_from_xml_element(self, root):
         if text is not None:
             self.weight_windows_on = text in ('true', '1')
 
+    def _shared_secondary_bank_from_xml_element(self, root):
+        text = get_text(root, 'shared_secondary_bank')
+        if text is not None:
+            self.shared_secondary_bank = text in ('true', '1')
+
     def _weight_windows_file_from_xml_element(self, root):
         text = get_text(root, 'weight_windows_file')
         if text is not None:
@@ -2538,6 +2567,7 @@ def to_xml_element(self, mesh_memo=None):
         self._create_write_initial_source_subelement(element)
         self._create_weight_windows_subelement(element, mesh_memo)
         self._create_weight_windows_on_subelement(element)
+        self._create_shared_secondary_bank_subelement(element)
         self._create_weight_window_generators_subelement(element, mesh_memo)
         self._create_weight_windows_file_element(element)
         self._create_weight_window_checkpoints_subelement(element)
@@ -2654,6 +2684,7 @@ def from_xml_element(cls, elem, meshes=None):
         settings._write_initial_source_from_xml_element(elem)
         settings._weight_windows_from_xml_element(elem, meshes)
         settings._weight_windows_on_from_xml_element(elem)
+        settings._shared_secondary_bank_from_xml_element(elem)
         settings._weight_windows_file_from_xml_element(elem)
         settings._weight_window_generators_from_xml_element(elem, meshes)
         settings._weight_window_checkpoints_from_xml_element(elem)