converted event based queues and particle buffer into SharedArrays

include/openmc/event.h

@@ -5,6 +5,7 @@
 //! \brief Event-based data structures and methods
 
 #include "openmc/particle.h"
+#include "openmc/shared_array.h"
 
 
 namespace openmc {
@@ -49,23 +50,14 @@ namespace simulation {
 // coordinated with atomics. This means that normal vector methods (e.g.,
 // push_back(), size()) would cause undefined or unintended behavior. Rather,
 // adding particles to queues will be done via the enqueue_particle() function.
-extern std::unique_ptr<QueueItem[]> calculate_fuel_xs_queue;
-extern std::unique_ptr<QueueItem[]> calculate_nonfuel_xs_queue;
-extern std::unique_ptr<QueueItem[]> advance_particle_queue;
-extern std::unique_ptr<QueueItem[]> surface_crossing_queue;
-extern std::unique_ptr<QueueItem[]> collision_queue;
+extern SharedArray<QueueItem> calculate_fuel_xs_queue;
+extern SharedArray<QueueItem> calculate_nonfuel_xs_queue;
+extern SharedArray<QueueItem> advance_particle_queue;
+extern SharedArray<QueueItem> surface_crossing_queue;
+extern SharedArray<QueueItem> collision_queue;
 
-// Event queue lengths
-extern int64_t calculate_fuel_xs_queue_length;
-extern int64_t calculate_nonfuel_xs_queue_length;
-extern int64_t advance_particle_queue_length;
-extern int64_t surface_crossing_queue_length;
-extern int64_t collision_queue_length;
-
-// Particle buffer. This is an allocated pointer rather than a vector as it
-// will be shared between threads, so most vector methods would result in
-// undefined or unintended behavior.
-extern std::unique_ptr<Particle[]>  particles;
+// Particle buffer
+extern std::vector<Particle>  particles;
 
 } // namespace simulation
 
@@ -82,11 +74,9 @@ void free_event_queues(void);
 
 //! Atomically adds a particle to the specified queue
 //! \param queue The queue to append the particle to
-//! \param length A reference to the length variable for the queue
 //! \param p A pointer to the particle
 //! \param buffer_idx The particle's actual index in the particle buffer
-void enqueue_particle(QueueItem* queue, int64_t& length, const Particle* p,
-    int64_t buffer_idx);
+void enqueue_particle(SharedArray<QueueItem>& queue, const Particle* p, int64_t buffer_idx);
 
 //! Enqueues a particle based on if it is in fuel or a non-fuel material
 //! \param buffer_idx The particle's actual index in the particle buffer
@@ -98,9 +88,8 @@ void dispatch_xs_event(int64_t buffer_idx);
 void process_init_events(int64_t n_particles, int64_t source_offset);
 
 //! Executes the calculate XS event for all particles in this event's buffer
-//! \param queue The XS lookup queue to use
-//! \param n_particles The number of particles in this queue
-void process_calculate_xs_events(QueueItem* queue, int64_t n_particles);
+//! \param queue A reference to the desired XS lookup queue
+void process_calculate_xs_events(SharedArray<QueueItem>& queue);
 
 //! Executes the advance particle event for all particles in this event's buffer
 void process_advance_particle_events();

include/openmc/shared_array.h

@@ -29,7 +29,6 @@ public:
   SharedArray(int64_t capacity) : capacity_(capacity)
   {
     data_ = std::make_unique<T[]>(capacity);
-    size_ = 0;
   }
   
   //==========================================================================
@@ -45,7 +44,7 @@ public:
 
   //! Increases the size of the SharedArray by one and returns an index to the
   //! last element of the array.
-  int64_t append()
+  int64_t thread_safe_append()
   {
     // Atomically capture the index we want to write to
     int64_t idx;
@@ -70,9 +69,13 @@ public:
   }
 
   int64_t size() {return size_;}
+  
+  int64_t resize(int64_t size) {size_ = size;}
 
   int64_t capacity() {return capacity_;}
 
+  T* data() {return data_.get();}
+
 }; 
 
 } // namespace openmc

src/event.cpp

@@ -12,19 +12,13 @@ namespace openmc {
 
 namespace simulation {
 
-std::unique_ptr<QueueItem[]> calculate_fuel_xs_queue;
-std::unique_ptr<QueueItem[]> calculate_nonfuel_xs_queue;
-std::unique_ptr<QueueItem[]> advance_particle_queue;
-std::unique_ptr<QueueItem[]> surface_crossing_queue;
-std::unique_ptr<QueueItem[]> collision_queue;
+SharedArray<QueueItem> calculate_fuel_xs_queue;
+SharedArray<QueueItem> calculate_nonfuel_xs_queue;
+SharedArray<QueueItem> advance_particle_queue;
+SharedArray<QueueItem> surface_crossing_queue;
+SharedArray<QueueItem> collision_queue;
 
-int64_t calculate_fuel_xs_queue_length {0};
-int64_t calculate_nonfuel_xs_queue_length {0};
-int64_t advance_particle_queue_length {0};
-int64_t surface_crossing_queue_length {0};
-int64_t collision_queue_length {0};
-
-std::unique_ptr<Particle[]>  particles;
+std::vector<Particle>  particles;
 
 } // namespace simulation
 
@@ -34,30 +28,29 @@ std::unique_ptr<Particle[]>  particles;
 
 void init_event_queues(int64_t n_particles)
 {
-  simulation::calculate_fuel_xs_queue =    std::make_unique<QueueItem[]>(n_particles);
-  simulation::calculate_nonfuel_xs_queue = std::make_unique<QueueItem[]>(n_particles);
-  simulation::advance_particle_queue =     std::make_unique<QueueItem[]>(n_particles);
-  simulation::surface_crossing_queue =     std::make_unique<QueueItem[]>(n_particles);
-  simulation::collision_queue =            std::make_unique<QueueItem[]>(n_particles);
-  simulation::particles =                  std::make_unique<Particle[] >(n_particles);
+  simulation::calculate_fuel_xs_queue.reserve(n_particles);
+  simulation::calculate_nonfuel_xs_queue.reserve(n_particles);
+  simulation::advance_particle_queue.reserve(n_particles);
+  simulation::surface_crossing_queue.reserve(n_particles);
+  simulation::collision_queue.reserve(n_particles);
+
+  simulation::particles.resize(n_particles);
 }
 
 void free_event_queues(void)
 {
-  simulation::calculate_fuel_xs_queue.reset();
-  simulation::calculate_nonfuel_xs_queue.reset();
-  simulation::advance_particle_queue.reset();
-  simulation::surface_crossing_queue.reset();
-  simulation::collision_queue.reset();
-  simulation::particles.reset();
+  simulation::calculate_fuel_xs_queue.clear();
+  simulation::calculate_nonfuel_xs_queue.clear();
+  simulation::advance_particle_queue.clear();
+  simulation::surface_crossing_queue.clear();
+  simulation::collision_queue.clear();
+
+  simulation::particles.clear();
 }
 
-void enqueue_particle(QueueItem* queue, int64_t& length, const Particle* p,
-    int64_t buffer_idx)
+void enqueue_particle(SharedArray<QueueItem>& queue, const Particle* p, int64_t buffer_idx)
 {
-  int64_t idx;
-  #pragma omp atomic capture
-  idx = length++;
+  int64_t idx = queue.thread_safe_append();
 
   queue[idx].idx = buffer_idx;
   queue[idx].E = p->E_;
@@ -68,13 +61,10 @@ void enqueue_particle(QueueItem* queue, int64_t& length, const Particle* p,
 void dispatch_xs_event(int64_t buffer_idx)
 {
   Particle* p = &simulation::particles[buffer_idx];
-  if (p->material_ == MATERIAL_VOID ||
-      !model::materials[p->material_]->fissionable_) {
-    enqueue_particle(simulation::calculate_nonfuel_xs_queue.get(),
-        simulation::calculate_nonfuel_xs_queue_length, p, buffer_idx);
+  if (p->material_ == MATERIAL_VOID || !model::materials[p->material_]->fissionable_) {
+    enqueue_particle(simulation::calculate_nonfuel_xs_queue, p, buffer_idx);
   } else {
-    enqueue_particle(simulation::calculate_fuel_xs_queue.get(),
-        simulation::calculate_fuel_xs_queue_length, p, buffer_idx);
+    enqueue_particle(simulation::calculate_fuel_xs_queue, p, buffer_idx);
   }
 }
 
@@ -89,7 +79,7 @@ void process_init_events(int64_t n_particles, int64_t source_offset)
   simulation::time_event_init.stop();
 }
 
-void process_calculate_xs_events(QueueItem* queue, int64_t n_particles)
+void process_calculate_xs_events(SharedArray<QueueItem>& queue)
 {
   simulation::time_event_calculate_xs.start();
 
@@ -99,21 +89,24 @@ void process_calculate_xs_events(QueueItem* queue, int64_t n_particles)
   // to C++17, std::sort is a serial only operation, which in this case
   // makes it too slow to be practical for most test problems. 
   //
-  // std::sort(std::execution::par_unseq, queue, queue+n);
+  // std::sort(std::execution::par_unseq, queue.data(), queue.data() + n);
+  
+  int64_t offset = simulation::advance_particle_queue.size();;
 
   #pragma omp parallel for schedule(runtime)
-  for (int64_t i = 0; i < n_particles; i++) {
+  for (int64_t i = 0; i < queue.size(); i++) {
     Particle* p = &simulation::particles[queue[i].idx]; 
     p->event_calculate_xs();
     
     // After executing a calculate_xs event, particles will
     // always require an advance event. Therefore, we don't need to use
     // the protected enqueuing function.
-    int64_t offset = simulation::advance_particle_queue_length + i;
-    simulation::advance_particle_queue[offset] = queue[i];
+    simulation::advance_particle_queue[offset + i] = queue[i];
   }
 
-  simulation::advance_particle_queue_length += n_particles;
+  simulation::advance_particle_queue.resize(offset + queue.size());
+
+  queue.resize(0);
 
   simulation::time_event_calculate_xs.stop();
 }
@@ -123,19 +116,19 @@ void process_advance_particle_events()
   simulation::time_event_advance_particle.start();
 
   #pragma omp parallel for schedule(runtime)
-  for (int64_t i = 0; i < simulation::advance_particle_queue_length; i++) {
+  for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) {
     int64_t buffer_idx = simulation::advance_particle_queue[i].idx;
     Particle* p = &simulation::particles[buffer_idx];
     p->event_advance();
     if (p->collision_distance_ > p->boundary_.distance) {
-      enqueue_particle(simulation::surface_crossing_queue.get(),
-          simulation::surface_crossing_queue_length, p, buffer_idx);
+      enqueue_particle(simulation::surface_crossing_queue, p, buffer_idx);
     } else {
-      enqueue_particle(simulation::collision_queue.get(),
-          simulation::collision_queue_length, p, buffer_idx);
+      enqueue_particle(simulation::collision_queue, p, buffer_idx);
     }
   }
 
+  simulation::advance_particle_queue.resize(0);
+
   simulation::time_event_advance_particle.stop();
 }
 
@@ -144,7 +137,7 @@ void process_surface_crossing_events()
   simulation::time_event_surface_crossing.start();
 
   #pragma omp parallel for schedule(runtime)
-  for (int64_t i = 0; i < simulation::surface_crossing_queue_length; i++) {
+  for (int64_t i = 0; i < simulation::surface_crossing_queue.size(); i++) {
     int64_t buffer_index = simulation::surface_crossing_queue[i].idx;
     Particle* p = &simulation::particles[buffer_index];
     p->event_cross_surface();
@@ -152,6 +145,8 @@ void process_surface_crossing_events()
     if (p->alive_)
       dispatch_xs_event(buffer_index);
   }
+
+  simulation::surface_crossing_queue.resize(0);
   
   simulation::time_event_surface_crossing.stop();
 }
@@ -161,7 +156,7 @@ void process_collision_events()
   simulation::time_event_collision.start();
 
   #pragma omp parallel for schedule(runtime)
-  for (int64_t i = 0; i < simulation::collision_queue_length; i++) {
+  for (int64_t i = 0; i < simulation::collision_queue.size(); i++) {
     int64_t buffer_index = simulation::collision_queue[i].idx;
     Particle* p = &simulation::particles[buffer_index];
     p->event_collide();
@@ -170,6 +165,8 @@ void process_collision_events()
       dispatch_xs_event(buffer_index);
   }
 
+  simulation::collision_queue.resize(0);
+
   simulation::time_event_collision.stop();
 }
 

src/physics.cpp

@@ -16,6 +16,7 @@
 #include "openmc/secondary_uncorrelated.h"
 #include "openmc/search.h"
 #include "openmc/settings.h"
+#include "openmc/shared_array.h"
 #include "openmc/simulation.h"
 #include "openmc/string_utils.h"
 #include "openmc/thermal.h"

src/simulation.cpp

@@ -626,32 +626,25 @@ void transport_event_based()
     while (true) {
       // Determine which event kernel has the longest queue
       int64_t max = std::max({
-        simulation::calculate_fuel_xs_queue_length,
-        simulation::calculate_nonfuel_xs_queue_length,
-        simulation::advance_particle_queue_length,
-        simulation::surface_crossing_queue_length,
-        simulation::collision_queue_length});
+        simulation::calculate_fuel_xs_queue.size(),
+        simulation::calculate_nonfuel_xs_queue.size(),
+        simulation::advance_particle_queue.size(),
+        simulation::surface_crossing_queue.size(),
+        simulation::collision_queue.size()});
 
       // Execute event with the longest queue
       if (max == 0) {
         break;
-      } else if (max == simulation::calculate_fuel_xs_queue_length) {
-        process_calculate_xs_events(simulation::calculate_fuel_xs_queue.get(),
-            simulation::calculate_fuel_xs_queue_length);
-        simulation::calculate_fuel_xs_queue_length = 0;
-      } else if (max == simulation::calculate_nonfuel_xs_queue_length) {
-        process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue.get(),
-            simulation::calculate_nonfuel_xs_queue_length);
-        simulation::calculate_nonfuel_xs_queue_length = 0;
-      } else if (max == simulation::advance_particle_queue_length) {
+      } else if (max == simulation::calculate_fuel_xs_queue.size()) {
+        process_calculate_xs_events(simulation::calculate_fuel_xs_queue);
+      } else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
+        process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
+      } else if (max == simulation::advance_particle_queue.size()) {
         process_advance_particle_events();
-        simulation::advance_particle_queue_length = 0;
-      } else if (max == simulation::surface_crossing_queue_length) {
+      } else if (max == simulation::surface_crossing_queue.size()) {
         process_surface_crossing_events();
-        simulation::surface_crossing_queue_length = 0;
-      } else if (max == simulation::collision_queue_length) {
+      } else if (max == simulation::collision_queue.size()) {
         process_collision_events();
-        simulation::collision_queue_length = 0;
       }
     }