diff --git a/include/openmc/event.h b/include/openmc/event.h index 9a39dec883..eeb657ec0a 100644 --- a/include/openmc/event.h +++ b/include/openmc/event.h @@ -1,6 +1,9 @@ #ifndef OPENMC_EVENT_H #define OPENMC_EVENT_H +//! \file event.h +//! \brief Event-based data structures and methods + #include "openmc/particle.h" #include "openmc/tallies/filter.h" @@ -11,11 +14,22 @@ namespace openmc { // Structs //============================================================================== +// In the event-based model, instead of moving or sorting the particles +// themselves based on which event they need, a queue is used to store the +// index (and other useful info) for each event type. +// The QueueItem struct holds the relevant information about a particle needed +// for sorting the queue. For very high particle counts, a sorted queue has the +// potential to result in greatly improved cache efficiency. However, sorting +// will introduce some overhead due to the sorting process itself, and may not +// result in any benefits if not enough particles are present for them to achieve +// consistent locality improvements. struct QueueItem{ - int64_t idx; // particle index in event-based buffer - double E; // particle energy - int64_t material; // material that particle is in - Particle::Type type; // particle type + int64_t idx; //!< particle index in event-based particle buffer + double E; //!< particle energy + int64_t material; //!< material that particle is in + Particle::Type type; //!< particle type + + // Comparator sorts by particle type, then by material type, then by energy bool operator<(const QueueItem& rhs) const { // First, compare by particle type @@ -30,15 +44,14 @@ struct QueueItem{ // Next, compare by material: // TODO: Currently, material IDs are not usually unique to material // types. When unique material type IDs are available, we can alter the - // material field in this struct to contain the type ID, and then sort - // by that, as below: - //if( material < rhs.material) - // return true; - //if( material > rhs.material) - // return false; + // material field in this struct to contain the 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 + // Now, compare by energy. return (E < rhs.E); } }; @@ -49,32 +62,74 @@ struct QueueItem{ namespace simulation { +// Event queues. These are allocated pointer variables rather than vectors, +// because they are shared between threads and writing to them must be +// 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 calculate_fuel_xs_queue; extern std::unique_ptr calculate_nonfuel_xs_queue; extern std::unique_ptr advance_particle_queue; extern std::unique_ptr surface_crossing_queue; extern std::unique_ptr collision_queue; -extern std::unique_ptr particles; + +// 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::vector particles; + } // namespace simulation //============================================================================== // Functions //============================================================================== +//! Allocates space for the event queues and particle buffer void init_event_queues(int64_t n_particles); + +//! Frees the event queues and particle buffer void free_event_queues(void); -void dispatch_xs_event(int64_t i); + +//! 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, 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 +void dispatch_xs_event(int64_t buffer_idx); + +//! Executes the initialization event for all particles +//! \param n_particles The number of particles in the particle buffer +//! \param source_offset The offset index in the source bank to use 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); + +//! Executes the advance particle event for all particles in this event's buffer void process_advance_particle_events(); + +//! Executes the surface crossing event for all particles in this event's buffer void process_surface_crossing_events(); + +//! Executes the collision event for all particles in this event's buffer void process_collision_events(); + +//! Executes the death event for all particles +//! \param n_particles The number of particles in the particle buffer void process_death_events(int64_t n_particles); } // namespace openmc diff --git a/src/event.cpp b/src/event.cpp index f56d8ec04a..2d94234e7a 100644 --- a/src/event.cpp +++ b/src/event.cpp @@ -17,7 +17,6 @@ std::unique_ptr calculate_nonfuel_xs_queue; std::unique_ptr advance_particle_queue; std::unique_ptr surface_crossing_queue; std::unique_ptr collision_queue; -std::unique_ptr particles; int64_t calculate_fuel_xs_queue_length {0}; int64_t calculate_nonfuel_xs_queue_length {0}; @@ -25,6 +24,8 @@ int64_t advance_particle_queue_length {0}; int64_t surface_crossing_queue_length {0}; int64_t collision_queue_length {0}; +std::unique_ptr particles; + } // namespace simulation //============================================================================== @@ -64,16 +65,16 @@ void enqueue_particle(QueueItem* queue, int64_t& length, Particle* p, queue[idx].type = p->type_; } -void dispatch_xs_event(int64_t i) +void dispatch_xs_event(int64_t buffer_idx) { - Particle* p = &simulation::particles[i]; + 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, i); + simulation::calculate_nonfuel_xs_queue_length, p, buffer_idx); } else { enqueue_particle(simulation::calculate_fuel_xs_queue.get(), - simulation::calculate_fuel_xs_queue_length, p, i); + simulation::calculate_fuel_xs_queue_length, p, buffer_idx); } } @@ -94,8 +95,11 @@ void process_calculate_xs_events(QueueItem* queue, int64_t n_particles) // TODO: If using C++17, perform a parallel sort of the queue // by particle type, material type, and then energy, in order to - // improve cache locality and reduce thread divergence on GPU. - //std::sort(queue, queue+n); + // improve cache locality and reduce thread divergence on GPU. Prior + // 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); #pragma omp parallel for schedule(runtime) for (int64_t i = 0; i < n_particles; i++) {