General import of event-based model from old branch. Seems to work, but need to test then work on cleaning up some more.

CMakeLists.txt

@@ -19,6 +19,7 @@ option(debug    "Compile with debug flags"                       OFF)
 option(optimize "Turn on all compiler optimization flags"        OFF)
 option(coverage "Compile with coverage analysis flags"           OFF)
 option(dagmc    "Enable support for DAGMC (CAD) geometry"        OFF)
+option(event    "Enable event-based transport mode"              OFF)
 
 #===============================================================================
 # MPI for distributed-memory parallelism
@@ -223,6 +224,7 @@ list(APPEND libopenmc_SOURCES
   src/eigenvalue.cpp
   src/endf.cpp
   src/error.cpp
+  src/event.cpp
   src/initialize.cpp
   src/finalize.cpp
   src/geometry.cpp
@@ -352,6 +354,10 @@ if(dagmc)
   target_include_directories(libopenmc PRIVATE ${DAGMC_INCLUDE_DIRS})
 endif()
 
+if(event)
+  target_compile_definitions(libopenmc PRIVATE EVENT_BASED)
+endif()
+
 #===============================================================================
 # openmc executable
 #===============================================================================

include/openmc/event.h

@@ -0,0 +1,84 @@
+#ifndef OPENMC_EVENT_H
+#define OPENMC_EVENT_H
+
+#include "openmc/particle.h"
+#include "openmc/tallies/filter.h"
+
+#include <vector>
+
+
+namespace openmc {
+
+//==============================================================================
+// Structs
+//==============================================================================
+
+struct QueueItem{
+  int64_t idx;      // particle index in event-based buffer
+  double E;     // particle energy
+  int material; // material that particle is in
+  Particle::Type type;
+  bool operator<(const QueueItem& rhs) const
+  {
+    // First, compare by type
+    if (type < rhs.type) {
+      return true;
+    } else if (type > rhs.type) {
+      return false;
+    }
+
+    // At this point, we have the same particle types.
+    // Now, compare by material
+
+    // TODO: Temporarily disabled as SMR problem has different material IDs for every pin
+    // Need to sort by material type instead...
+    /*
+       if( material < rhs.material)
+       return true;
+       if( material > rhs.material)
+       return false;
+       */
+
+    // At this point, we have the same particle type, in the same material.
+    // Now, compare by energy
+    return (E < rhs.E);
+  }
+};
+
+//==============================================================================
+// Global variable declarations
+//==============================================================================
+//
+namespace simulation {
+
+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 std::unique_ptr<Particle[]>  particles;
+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;
+extern int64_t max_particles_in_flight;
+
+} // namespace simulation
+
+//==============================================================================
+// Functions
+//==============================================================================
+
+void init_event_queues(int64_t n_particles);
+void free_event_queues(void);
+void dispatch_xs_event(int64_t i);
+void process_calculate_xs_events(QueueItem* queue, int64_t n);
+void process_advance_particle_events();
+void process_surface_crossing_events();
+void process_collision_events();
+void process_death_events(int64_t n_particles);
+
+} // namespace openmc
+
+#endif // OPENMC_EVENT_H

include/openmc/simulation.h

@@ -95,6 +95,9 @@ void transport_history_based_single_particle(Particle& p);
 //! Simulate all particle histories using history-based parallelism
 void transport_history_based();
 
+//! Simulate all particle histories using event-based parallelism
+void transport_event_based();
+
 } // namespace openmc
 
 #endif // OPENMC_SIMULATION_H

src/event.cpp

@@ -0,0 +1,176 @@
+#include "openmc/event.h"
+#include "openmc/material.h"
+
+namespace openmc {
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+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);
+}
+
+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();
+}
+
+void dispatch_xs_event(int64_t i)
+{
+  Particle* p = &simulation::particles[i];
+  int64_t idx;
+  if (p->material_ == MATERIAL_VOID) {
+    #pragma omp atomic capture
+    idx = simulation::calculate_nonfuel_xs_queue_length++;
+    simulation::calculate_nonfuel_xs_queue[idx].idx = i;
+    simulation::calculate_nonfuel_xs_queue[idx].E = p->E_;
+    simulation::calculate_nonfuel_xs_queue[idx].material = p->material_;
+    simulation::calculate_nonfuel_xs_queue[idx].type = p->type_;
+  }
+  else
+  {
+    if (model::materials[p->material_]->fissionable_) {
+      #pragma omp atomic capture
+      idx = simulation::calculate_fuel_xs_queue_length++;
+      simulation::calculate_fuel_xs_queue[idx].idx = i;
+      simulation::calculate_fuel_xs_queue[idx].E = p->E_;
+      simulation::calculate_fuel_xs_queue[idx].material = p->material_;
+      simulation::calculate_fuel_xs_queue[idx].type = p->type_;
+    }
+    else
+    {
+      #pragma omp atomic capture
+      idx = simulation::calculate_nonfuel_xs_queue_length++;
+      simulation::calculate_nonfuel_xs_queue[idx].idx = i;
+      simulation::calculate_nonfuel_xs_queue[idx].E = p->E_;
+      simulation::calculate_nonfuel_xs_queue[idx].material = p->material_;
+      simulation::calculate_nonfuel_xs_queue[idx].type = p->type_;
+    }
+  }
+}
+
+void process_calculate_xs_events(QueueItem* queue, int64_t n)
+{
+  // Sort queue by energy
+  std::sort(queue, queue+n);
+
+  // Save last_ members, find grid index
+  #pragma omp parallel for schedule(runtime)
+  for (auto i = 0; i < n; i++) {
+    Particle* p = &simulation::particles[queue[i].idx]; 
+    p->event_calculate_xs();
+  }
+
+  int64_t start = simulation::advance_particle_queue_length;
+  int64_t end = start + n;
+  int64_t j = 0;
+  for (auto i = start; i < end; i++) {
+    simulation::advance_particle_queue[i].idx = queue[j].idx;
+    simulation::advance_particle_queue[i].E = simulation::particles[queue[j].idx].E_;
+    simulation::advance_particle_queue[i].material = simulation::particles[queue[j].idx].material_;
+    simulation::advance_particle_queue[i].type = simulation::particles[queue[j].idx].type_;
+    j++;
+  }
+  simulation::advance_particle_queue_length += n;
+}
+
+void process_advance_particle_events()
+{
+  #pragma omp parallel for schedule(runtime)
+  for (auto i = 0; i < simulation::advance_particle_queue_length; i++) {
+    Particle* p = &simulation::particles[simulation::advance_particle_queue[i].idx];
+    p->event_advance();
+    if( p->collision_distance_ > p->boundary_.distance ) 
+    {
+      int64_t idx;
+      #pragma omp atomic capture
+      idx = simulation::surface_crossing_queue_length++;
+      simulation::surface_crossing_queue[idx].idx = simulation::advance_particle_queue[i].idx;
+      simulation::surface_crossing_queue[idx].E = p->E_;
+      simulation::surface_crossing_queue[idx].material = p->material_;
+      simulation::surface_crossing_queue[idx].type = p->type_;
+    }
+    else
+    {
+      int64_t idx;
+      #pragma omp atomic capture
+      idx = simulation::collision_queue_length++;
+      simulation::collision_queue[idx].idx = simulation::advance_particle_queue[i].idx;
+      simulation::collision_queue[idx].E = p->E_;
+      simulation::collision_queue[idx].material = p->material_;
+      simulation::collision_queue[idx].type = p->type_;
+    }
+  }
+}
+
+void process_surface_crossing_events()
+{
+  #pragma omp parallel for schedule(runtime)
+  for (auto i = 0; i < simulation::surface_crossing_queue_length; i++) {
+    Particle* p = &simulation::particles[simulation::surface_crossing_queue[i].idx];
+    p->event_cross_surface();
+    p->event_revive_from_secondary();
+    if (p->alive_)
+      dispatch_xs_event(simulation::surface_crossing_queue[i].idx);
+  }
+
+}
+
+void process_collision_events()
+{
+  #pragma omp parallel for schedule(runtime)
+  for (auto i = 0; i < simulation::collision_queue_length; i++) {
+    Particle* p = &simulation::particles[simulation::collision_queue[i].idx];
+    p->event_collide();
+    p->event_revive_from_secondary();
+    if (p->alive_)
+      dispatch_xs_event(simulation::collision_queue[i].idx);
+  }
+
+}
+
+void process_death_events(int64_t n_particles)
+{
+  #pragma omp parallel for schedule(runtime)
+  for (auto i = 0; i < n_particles; i++) {
+    Particle& p = simulation::particles[i];
+    p.event_death();
+  }
+}
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+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;
+std::unique_ptr<Particle[]>  particles;
+
+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};
+
+int64_t max_particles_in_flight {100000};
+
+} // namespace simulation
+
+} // namespace openmc

src/finalize.cpp

@@ -7,6 +7,7 @@
 #include "openmc/cross_sections.h"
 #include "openmc/dagmc.h"
 #include "openmc/eigenvalue.h"
+#include "openmc/event.h"
 #include "openmc/geometry.h"
 #include "openmc/geometry_aux.h"
 #include "openmc/material.h"
@@ -50,6 +51,9 @@ void free_memory()
 #ifdef DAGMC
   free_memory_dagmc();
 #endif
+#ifdef EVENT_BASED
+  free_event_queues();
+#endif
 }
 
 }

src/simulation.cpp

@@ -5,6 +5,7 @@
 #include "openmc/container_util.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
+#include "openmc/event.h"
 #include "openmc/geometry_aux.h"
 #include "openmc/material.h"
 #include "openmc/message_passing.h"
@@ -72,6 +73,14 @@ int openmc_simulation_init()
   // fission bank
   allocate_banks();
 
+  // If doing an event-based simulation, intialize the particle buffer
+  // and event queues
+  #ifdef EVENT_BASED
+  int64_t event_buffer_length = std::min(simulation::work_per_rank,
+      simulation::max_particles_in_flight);
+  init_event_queues(event_buffer_length);
+  #endif
+
   // Allocate tally results arrays if they're not allocated yet
   for (auto& t : model::tallies) {
     t->init_results();
@@ -183,7 +192,11 @@ int openmc_next_batch(int* status)
     simulation::time_transport.start();
 
     // Transport loop
+    #ifdef EVENT_BASED
+    transport_event_based();
+    #else
     transport_history_based();
+    #endif
 
     // Accumulate time for transport
     simulation::time_transport.stop();
@@ -593,4 +606,73 @@ void transport_history_based()
   }
 }
 
+void transport_event_based()
+{
+  int64_t remaining_work = simulation::work_per_rank;
+  int64_t source_offset = 0;
+
+  // To cap the total amount of memory used to store particle object data, the
+  // number of particles in flight at any point in time can set. In the case
+  // that the maximum in flight particle count is lower than the total number
+  // of particles that need to be run this iteration, the event-based transport
+  // loop is executed using subierations until all particles have been completed.
+  while (remaining_work > 0) {
+    // Figure out # of particles to run for this subiteration
+    int64_t n_particles = std::min(remaining_work, simulation::max_particles_in_flight);
+
+    // Initialize all particle histories for this subiteration
+    #pragma omp parallel for schedule(runtime)
+    for (auto i = 0; i < n_particles; i++) {
+      initialize_history(&simulation::particles[i], source_offset + i + 1);
+    }
+
+    // Add all particles to the XS lookup queue
+    #pragma omp parallel for schedule(runtime)
+    for (auto i = 0; i < n_particles; i++) {
+      dispatch_xs_event(i);
+    }
+
+    int events_retired = 0;
+
+    // Event-based transport loop
+    while (true) {
+      // Determine which event kernel has the most particles in its 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});
+
+      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) {
+        process_advance_particle_events();
+        simulation::advance_particle_queue_length = 0;
+      } else if (max == simulation::surface_crossing_queue_length) {
+        process_surface_crossing_events();
+        simulation::surface_crossing_queue_length = 0;
+      } else if (max == simulation::collision_queue_length) {
+        process_collision_events();
+        simulation::collision_queue_length = 0;
+      }
+
+      events_retired++;
+    }
+    
+    process_death_events(n_particles);
+
+    remaining_work -= n_particles;
+    source_offset += n_particles;
+  }
+}
+
 } // namespace openmc