added setting to control the maximum number of in flight particles. Also simplified enqueing procedure after XS lookup events to drop atomic usage there.

include/openmc/event.h

@@ -60,7 +60,6 @@ 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
 
@@ -72,7 +71,7 @@ void init_event_queues(int64_t n_particles);
 void free_event_queues(void);
 void dispatch_xs_event(int64_t i);
 void process_init_events(int64_t n_particles, int64_t source_offset);
-void process_calculate_xs_events(QueueItem* queue, int64_t n);
+void process_calculate_xs_events(QueueItem* queue, int64_t n_particles);
 void process_advance_particle_events();
 void process_surface_crossing_events();
 void process_collision_events();

include/openmc/settings.h

@@ -67,6 +67,8 @@ extern "C" int32_t n_inactive;     //!< number of inactive batches
 extern "C" int32_t gen_per_batch;  //!< number of generations per batch
 extern "C" int64_t n_particles;    //!< number of particles per generation
 
+extern "C" int64_t max_in_flight_particles; //!< Max num. event-based particles in flight
+
 extern int electron_treatment;       //!< how to treat secondary electrons
 extern std::array<double, 4> energy_cutoff;  //!< Energy cutoff in [eV] for each particle type
 extern int legendre_to_tabular_points; //!< number of points to convert Legendres

src/event.cpp

@@ -25,8 +25,6 @@ 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
 
 //==============================================================================
@@ -54,15 +52,11 @@ void free_event_queues(void)
 }
 
 void enqueue_particle(QueueItem* queue, int64_t& length, Particle* p,
-    int64_t buffer_idx, bool use_atomic)
+    int64_t buffer_idx)
 {
   int64_t idx;
-  if (use_atomic) {
-    #pragma omp atomic capture
-    idx = length++;
-  } else {
-    idx = length++;
-  }
+  #pragma omp atomic capture
+  idx = length++;
 
   queue[idx].idx = buffer_idx;
   queue[idx].E = p->E_;
@@ -77,10 +71,10 @@ void dispatch_xs_event(int64_t i)
   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, i, true);
+        simulation::calculate_nonfuel_xs_queue_length, p, i);
   } else {
       enqueue_particle(simulation::calculate_fuel_xs_queue.get(),
-          simulation::calculate_fuel_xs_queue_length, p, i, true);
+          simulation::calculate_fuel_xs_queue_length, p, i);
   }
 }
 
@@ -95,7 +89,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)
+void process_calculate_xs_events(QueueItem* queue, int64_t n_particles)
 {
   simulation::time_event_calculate_xs.start();
 
@@ -105,13 +99,19 @@ void process_calculate_xs_events(QueueItem* queue, int64_t n)
   //std::sort(queue, queue+n);
 
   #pragma omp parallel for schedule(runtime)
-  for (int64_t i = 0; i < n; i++) {
+  for (int64_t i = 0; i < n_particles; i++) {
     Particle* p = &simulation::particles[queue[i].idx]; 
     p->event_calculate_xs();
-    enqueue_particle(simulation::advance_particle_queue.get(),
-        simulation::advance_particle_queue_length, p, queue[i].idx, true);
+    
+    // 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_length += n_particles;
+
   simulation::time_event_calculate_xs.stop();
 }
 
@@ -126,10 +126,10 @@ void process_advance_particle_events()
     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, true);
+          simulation::surface_crossing_queue_length, p, buffer_idx);
     } else {
       enqueue_particle(simulation::collision_queue.get(),
-          simulation::collision_queue_length, p, buffer_idx, true);
+          simulation::collision_queue_length, p, buffer_idx);
     }
   }
 

src/settings.cpp

@@ -81,6 +81,8 @@ int32_t n_inactive {0};
 int32_t gen_per_batch {1};
 int64_t n_particles {-1};
 
+int64_t max_in_flight_particles {100000};
+
 int electron_treatment {ELECTRON_TTB};
 std::array<double, 4> energy_cutoff {0.0, 1000.0, 0.0, 0.0};
 int legendre_to_tabular_points {C_NONE};
@@ -127,6 +129,11 @@ void get_run_parameters(pugi::xml_node node_base)
   if (n_particles == -1) {
     n_particles = std::stoll(get_node_value(node_base, "particles"));
   }
+  
+  // Get maximum number of in flight particles for event-based mode
+  if (check_for_node(node_base, "max_in_flight_particles")) {
+    max_in_flight_particles = std::stoll(get_node_value(node_base, "max_in_flight_particles"));
+  }
 
   // Get number of basic batches
   if (check_for_node(node_base, "batches")) {

src/simulation.cpp

@@ -77,7 +77,7 @@ int openmc_simulation_init()
   // and event queues
   #ifdef EVENT_BASED
   int64_t event_buffer_length = std::min(simulation::work_per_rank,
-      simulation::max_particles_in_flight);
+      settings::max_in_flight_particles);
   init_event_queues(event_buffer_length);
   #endif
 
@@ -617,14 +617,14 @@ void transport_event_based()
   // loop is executed multiple times 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);
+    int64_t n_particles = std::min(remaining_work, settings::max_in_flight_particles);
 
     // Initialize all particle histories for this subiteration
     process_init_events(n_particles, source_offset);
 
     // Event-based transport loop
     while (true) {
-      // Determine which event kernel has the most particles in its queue
+      // 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,