From 8c884a17de91c0cbf23eab30fa3c9128ca86e38b Mon Sep 17 00:00:00 2001 From: John Tramm Date: Sun, 19 Jan 2020 19:12:32 +0000 Subject: [PATCH] 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 | 6 ++ include/openmc/event.h | 84 +++++++++++++++++ include/openmc/simulation.h | 3 + src/event.cpp | 176 ++++++++++++++++++++++++++++++++++++ src/finalize.cpp | 4 + src/simulation.cpp | 82 +++++++++++++++++ 6 files changed, 355 insertions(+) create mode 100644 include/openmc/event.h create mode 100644 src/event.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 4e8211677..95b334ce4 100644 --- a/CMakeLists.txt +++ b/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 #=============================================================================== diff --git a/include/openmc/event.h b/include/openmc/event.h new file mode 100644 index 000000000..7a654b6af --- /dev/null +++ b/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 + + +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 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; +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 diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 5bfba138c..d02308fe1 100644 --- a/include/openmc/simulation.h +++ b/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 diff --git a/src/event.cpp b/src/event.cpp new file mode 100644 index 000000000..2dfbb3ed5 --- /dev/null +++ b/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(n_particles); + simulation::calculate_nonfuel_xs_queue = std::make_unique(n_particles); + simulation::advance_particle_queue = std::make_unique(n_particles); + simulation::surface_crossing_queue = std::make_unique(n_particles); + simulation::collision_queue = std::make_unique(n_particles); + simulation::particles = std::make_unique(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 calculate_fuel_xs_queue; +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}; +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 diff --git a/src/finalize.cpp b/src/finalize.cpp index 5bda587ac..16e768733 100644 --- a/src/finalize.cpp +++ b/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 } } diff --git a/src/simulation.cpp b/src/simulation.cpp index d4a0438e9..fef56dee8 100644 --- a/src/simulation.cpp +++ b/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