From 0fb047f56cabe7b95f7ee2ae5a7d36496f21ab42 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 15 Nov 2019 21:28:24 +0000 Subject: [PATCH 01/40] Brings things up to level with the @paulromano event-based code --- include/openmc/geometry.h | 2 + include/openmc/particle.h | 24 ++ include/openmc/simulation.h | 6 +- src/material.cpp | 6 + src/particle.cpp | 14 +- src/simulation.cpp | 469 ++++++++++++++++++++++++++++++++++++ 6 files changed, 516 insertions(+), 5 deletions(-) diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index e386348c4..a6e1e4063 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -28,12 +28,14 @@ extern std::vector overlap_check_count; // Information about nearest boundary crossing //============================================================================== +/* struct BoundaryInfo { double distance {INFINITY}; //!< distance to nearest boundary int surface_index {0}; //!< if boundary is surface, index in surfaces vector int coord_level; //!< coordinate level after crossing boundary std::array lattice_translation {}; //!< which way lattice indices will change }; +*/ //============================================================================== //! Check two distances by coincidence tolerance diff --git a/include/openmc/particle.h b/include/openmc/particle.h index a2b7ce557..38e12574c 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -128,8 +128,21 @@ struct MacroXS { double incoherent; //!< macroscopic incoherent xs double photoelectric; //!< macroscopic photoelectric xs double pair_production; //!< macroscopic pair production xs + int i_grid; }; +//============================================================================== +// Information about nearest boundary crossing +//============================================================================== + +struct BoundaryInfo { + double distance {INFINITY}; //!< distance to nearest boundary + int surface_index {0}; //!< if boundary is surface, index in surfaces vector + int coord_level; //!< coordinate level after crossing boundary + std::array lattice_translation {}; //!< which way lattice indices will change +}; + + //============================================================================ //! State of a particle being transported through geometry //============================================================================ @@ -272,6 +285,7 @@ public: // Indices for various arrays int surface_ {0}; //!< index for surface particle is on + BoundaryInfo boundary_; int cell_born_ {-1}; //!< index for cell particle was born in int material_ {-1}; //!< index for current material int material_last_ {-1}; //!< index for last material @@ -289,4 +303,14 @@ public: } // namespace openmc +extern template class std::vector; + +namespace openmc{ + +extern std::vector particle_bank; +#pragma omp threadprivate(particle_bank) + +} + + #endif // OPENMC_PARTICLE_H diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 8efac4c7f..5a652b35b 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -38,7 +38,9 @@ extern "C" int restart_batch; //!< batch at which a restart job resumed extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied? extern "C" int total_gen; //!< total number of generations simulated extern double total_weight; //!< Total source weight in a batch -extern int64_t work_per_rank; //!< number of particles per MPI rank +extern int64_t thread_work_index; +extern int64_t work_per_rank; //!< number of particles per MPI rank +extern int64_t work_per_thread; //!< number of particles on a given tread extern const RegularMesh* entropy_mesh; extern const RegularMesh* ufs_mesh; @@ -49,7 +51,7 @@ extern std::vector work_index; // Threadprivate variables extern "C" bool trace; //!< flag to show debug information -#pragma omp threadprivate(current_work, trace) +#pragma omp threadprivate(current_work, work_per_thread, trace) } // namespace simulation diff --git a/src/material.cpp b/src/material.cpp index a5bc3dc13..16b9275a0 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -751,6 +751,7 @@ void Material::calculate_xs(Particle& p) const void Material::calculate_neutron_xs(Particle& p) const { + /* // Find energy index on energy grid int neutron = static_cast(Particle::Type::neutron); int i_grid = std::log(p.E_/data::energy_min[neutron])/simulation::log_spacing; @@ -760,9 +761,11 @@ void Material::calculate_neutron_xs(Particle& p) const // Initialize position in i_sab_nuclides int j = 0; + */ // Add contribution from each nuclide in material for (int i = 0; i < nuclide_.size(); ++i) { + /* // ====================================================================== // CHECK FOR S(A,B) TABLE @@ -792,18 +795,21 @@ void Material::calculate_neutron_xs(Particle& p) const // ====================================================================== // CALCULATE MICROSCOPIC CROSS SECTION + */ // Determine microscopic cross sections for this nuclide int i_nuclide = nuclide_[i]; // Calculate microscopic cross section for this nuclide const auto& micro {p.neutron_xs_[i_nuclide]}; + /* if (p.E_ != micro.last_E || p.sqrtkT_ != micro.last_sqrtkT || i_sab != micro.index_sab || sab_frac != micro.sab_frac) { data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, p); } + */ // ====================================================================== // ADD TO MACROSCOPIC CROSS SECTION diff --git a/src/particle.cpp b/src/particle.cpp index 030c8dc63..a4642ccac 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -27,8 +27,14 @@ #include "openmc/tallies/tally_scoring.h" #include "openmc/track_output.h" +// Explicit template instantiation definition +template class std::vector; + + namespace openmc { +std::vector particle_bank; + //============================================================================== // LocalCoord implementation //============================================================================== @@ -132,14 +138,15 @@ void Particle::transport() { // Display message if high verbosity or trace is on - if (settings::verbosity >= 9 || simulation::trace) { - write_message("Simulating Particle " + std::to_string(id_)); - } + //if (settings::verbosity >= 9 || simulation::trace) { + // write_message("Simulating Particle " + std::to_string(id_)); + //} // Initialize number of events to zero int n_event = 0; // Add paricle's starting weight to count for normalizing tallies later + /* #pragma omp atomic simulation::total_weight += wgt_; @@ -153,6 +160,7 @@ Particle::transport() // Every particle starts with no accumulated flux derivative. if (!model::active_tallies.empty()) zero_flux_derivs(); + */ while (true) { // Set the random number stream diff --git a/src/simulation.cpp b/src/simulation.cpp index 8de1e7fbc..1b566f518 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -2,23 +2,32 @@ #include "openmc/bank.h" #include "openmc/capi.h" +#include "openmc/cell.h" #include "openmc/container_util.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" +#include "openmc/geometry.h" #include "openmc/material.h" #include "openmc/message_passing.h" +#include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" #include "openmc/output.h" #include "openmc/particle.h" #include "openmc/photon.h" +#include "openmc/physics.h" + #include "openmc/physics_mg.h" #include "openmc/random_lcg.h" #include "openmc/settings.h" #include "openmc/source.h" #include "openmc/state_point.h" +#include "openmc/thermal.h" #include "openmc/timer.h" +#include "openmc/tallies/derivative.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/tally.h" +#include "openmc/tallies/tally_scoring.h" #include "openmc/tallies/trigger.h" +#include "openmc/track_output.h" #ifdef _OPENMP #include @@ -28,6 +37,408 @@ #include #include +namespace openmc { + +extern std::vector calculate_fuel_xs_queue; +extern std::vector calculate_nonfuel_xs_queue; +extern std::vector advance_particle_queue; +extern std::vector surface_crossing_queue; +extern std::vector collision_queue; +#pragma omp threadprivate(calculate_fuel_xs_queue, calculate_nonfuel_xs_queue, advance_particle_queue, surface_crossing_queue, collision_queue) + +std::vector calculate_fuel_xs_queue; +std::vector calculate_nonfuel_xs_queue; +std::vector advance_particle_queue; +std::vector surface_crossing_queue; +std::vector collision_queue; + +constexpr size_t MAX_PARTICLES_PER_THREAD {100}; + +void initialize_histories(int& index_source, + size_t& remaining_work_per_thread) +{ + int work = std::min(remaining_work_per_thread, MAX_PARTICLES_PER_THREAD); + for (int i = 0; i < work; ++i) { + particle_bank.emplace_back(); + auto& p {particle_bank.back()}; + + initialize_history(&p, index_source); + ++index_source; + } + + remaining_work_per_thread -= work; +} + +void revive_particle_from_secondary(Particle* p) +{ + p->from_source(&simulation::secondary_bank.back()); + simulation::secondary_bank.pop_back(); + // n_event = 0; + + // Enter new particle in particle track file + if (p->write_track_) add_particle_track(); +} + +void dispatch_xs_event(Particle* p) +{ + if (p->material_ == MATERIAL_VOID) { + calculate_nonfuel_xs_queue.push_back(p); + } else { + if (model::materials[p->material_]->fissionable_) { + calculate_fuel_xs_queue.push_back(p); + } else { + calculate_nonfuel_xs_queue.push_back(p); + } + } +} + +void process_calculate_xs_events(std::vector& queue) +{ + // Save last_ members, find grid index + for (auto& p : queue) { + // Set the random number stream + if (p->type_ == Particle::Type::neutron) { + prn_set_stream(STREAM_TRACKING); + } else { + prn_set_stream(STREAM_PHOTON); + } + + // Store pre-collision particle properties + p->wgt_last_ = p->wgt_; + p->E_last_ = p->E_; + p->u_last_ = p->u(); + p->r_last_ = p->r(); + + // If the cell hasn't been determined based on the particle's location, + // initiate a search for the current cell. This generally happens at the + // beginning of the history and again for any secondary particles + if (p->coord_[p->n_coord_ - 1].cell == C_NONE) { + if (!find_cell(p, false)) { + p->mark_as_lost("Could not find the cell containing particle " + + std::to_string(p->id_)); + return; + } + + // set birth cell attribute + if (p->cell_born_ == C_NONE) p->cell_born_ = p->coord_[p->n_coord_ - 1].cell; + } + + // Write particle track. + if (p->write_track_) write_particle_track(*p); + + if (settings::check_overlaps) check_cell_overlap(p); + + if (settings::run_CE) { + if (p->material_ == p->material_last_ && p->sqrtkT_ != p->sqrtkT_last_) { + // Remove particle from queue + } + } + + // Find energy index on energy grid + // TODO: Calculate this separately? + int neutron = static_cast(Particle::Type::neutron); + p->macro_xs_.i_grid = std::log(p->E_/data::energy_min[neutron]) / simulation::log_spacing; + } + + // Calculate nuclide micros + for (int i = 0; i < data::nuclides.size(); ++i) { + for (auto& p : queue) { + if (p->material_ == MATERIAL_VOID) continue; + + // If material doesn't have this nuclide, skip it + const auto& mat {model::materials[p->material_]}; + if (mat->mat_nuclide_index_[i] == -1) continue; + + // ====================================================================== + // CHECK FOR S(A,B) TABLE + + // Check if this nuclide matches one of the S(a,b) tables specified. + // This relies on thermal_tables_ being sorted by .index_nuclide + int i_sab = C_NONE; + double sab_frac = 0.0; + for (const auto& sab : mat->thermal_tables_) { + if (i == sab.index_nuclide) { + // Get index in sab_tables + i_sab = sab.index_table; + sab_frac = sab.fraction; + + // If particle energy is greater than the highest energy for the + // S(a,b) table, then don't use the S(a,b) table + //if (p->E_ > data::thermal_scatt[i_sab]->threshold()) i_sab = C_NONE; + if (p->E_ > data::thermal_scatt[i_sab]->energy_max_) i_sab = C_NONE; + } + } + + // ====================================================================== + // CALCULATE MICROSCOPIC CROSS SECTION + + // Calculate microscopic cross section for this nuclide + const auto& micro {p->neutron_xs_[i]}; + if (p->E_ != micro.last_E + || p->sqrtkT_ != micro.last_sqrtkT + || i_sab != micro.index_sab + || sab_frac != micro.sab_frac) { + data::nuclides[i]->calculate_xs(i_sab, p->macro_xs_.i_grid, sab_frac, *p); + } + } + } + + for (auto& p : queue) { + // Calculate microscopic and macroscopic cross sections + if (p->material_ != MATERIAL_VOID) { + if (settings::run_CE) { + // If the material is the same as the last material and the + // temperature hasn't changed, we don't need to lookup cross + // sections again. + model::materials[p->material_]->calculate_xs(*p); + } else { + // Get the MG data + calculate_xs_c(p->material_, p->g_, p->sqrtkT_, p->u_local(), + p->macro_xs_.total, p->macro_xs_.absorption, p->macro_xs_.nu_fission); + + // Finally, update the particle group while we have already checked + // for if multi-group + p->g_last_ = p->g_; + } + } else { + p->macro_xs_.total = 0.0; + p->macro_xs_.absorption = 0.0; + p->macro_xs_.fission = 0.0; + p->macro_xs_.nu_fission = 0.0; + } + + advance_particle_queue.push_back(p); + } + + queue.clear(); +} + +void process_advance_particle_events() +{ + for (auto& p : advance_particle_queue) { + simulation::trace == (p->id_ == 0); + + // Sample a distance to collision + double d_collision; + if (p->type_ == Particle::Type::electron || + p->type_ == Particle::Type::positron) { + d_collision = 0.0; + } else if (p->macro_xs_.total == 0.0) { + d_collision = INFINITY; + } else { + d_collision = -std::log(prn()) / p->macro_xs_.total; + } + + // -------------- break here? ------------------- + + // Find the distance to the nearest boundary + p->boundary_ = distance_to_boundary(p); + + // Select smaller of the two distances + double distance; + if (p->boundary_.distance < d_collision) { + surface_crossing_queue.push_back(p); + distance = p->boundary_.distance; + } else { + collision_queue.push_back(p); + distance = d_collision; + } + + // -------------- break here? ------------------- + + // Advance particle + for (int j = 0; j < p->n_coord_; ++j) { + p->coord_[j].r += distance * p->coord_[j].u; + } + + // -------------- break here? ------------------- + + // Score track-length tallies + if (!model::active_tracklength_tallies.empty()) { + score_tracklength_tally(p, distance); + } + + // Score track-length estimate of k-eff + if (settings::run_mode == RUN_MODE_EIGENVALUE && + p->type_ == Particle::Type::neutron) { + global_tally_tracklength += p->wgt_ * distance * p->macro_xs_.nu_fission; + } + + // Score flux derivative accumulators for differential tallies. + if (!model::active_tallies.empty()) { + score_track_derivative(p, distance); + } + } + + advance_particle_queue.clear(); +} + +void process_surface_crossing_events() +{ + for (auto& p : surface_crossing_queue) { + // Set surface that particle is on and adjust coordinate levels + p->surface_ = p->boundary_.surface_index; + p->n_coord_ = p->boundary_.coord_level; + + // Saving previous cell data + for (int j = 0; j < p->n_coord_; ++j) { + p->cell_last_[j] = p->coord_[j].cell; + } + p->n_coord_last_ = p->n_coord_; + + if (p->boundary_.lattice_translation[0] != 0 || + p->boundary_.lattice_translation[1] != 0 || + p->boundary_.lattice_translation[2] != 0) { + // Particle crosses lattice boundary + cross_lattice(p, p->boundary_); + p->event_ = EVENT_LATTICE; + } else { + // Particle crosses surface + p->cross_surface(); + p->event_ = EVENT_SURFACE; + } + // Score cell to cell partial currents + if (!model::active_surface_tallies.empty()) { + score_surface_tally(p, model::active_surface_tallies); + } + + if (!p->alive_ && !simulation::secondary_bank.empty()) { + revive_particle_from_secondary(p); + } + + if (p->alive_) dispatch_xs_event(p); + } + + surface_crossing_queue.clear(); +} + +void process_collision_events() +{ + for (auto& p : collision_queue) { + // Score collision estimate of keff + if (settings::run_mode == RUN_MODE_EIGENVALUE && + p->type_ == Particle::Type::neutron) { + global_tally_collision += p->wgt_ * p->macro_xs_.nu_fission + / p->macro_xs_.total; + } + + // Score surface current tallies -- this has to be done before the collision + // since the direction of the particle will change and we need to use the + // pre-collision direction to figure out what mesh surfaces were crossed + + if (!model::active_meshsurf_tallies.empty()) + score_surface_tally(p, model::active_meshsurf_tallies); + + // Clear surface component + p->surface_ = 0; + + if (settings::run_CE) { + collision(p); + } else { + collision_mg(p); + } + + // Score collision estimator tallies -- this is done after a collision + // has occurred rather than before because we need information on the + // outgoing energy for any tallies with an outgoing energy filter + if (!model::active_collision_tallies.empty()) score_collision_tally(p); + if (!model::active_analog_tallies.empty()) { + if (settings::run_CE) { + score_analog_tally_ce(p); + } else { + score_analog_tally_mg(p); + } + } + + // Reset banked weight during collision + p->n_bank_ = 0; + p->wgt_bank_ = 0.0; + for (int& v : p->n_delayed_bank_) v = 0; + + // Reset fission logical + p->fission_ = false; + + // Save coordinates for tallying purposes + p->r_last_current_ = p->r(); + + // Set last material to none since cross sections will need to be + // re-evaluated + p->material_last_ = C_NONE; + + // Set all directions to base level -- right now, after a collision, only + // the base level directions are changed + for (int j = 0; j < p->n_coord_ - 1; ++j) { + if (p->coord_[j + 1].rotated) { + // If next level is rotated, apply rotation matrix + const auto& m {model::cells[p->coord_[j].cell]->rotation_}; + const auto& u {p->coord_[j].u}; + p->coord_[j + 1].u.x = m[3]*u.x + m[4]*u.y + m[5]*u.z; + p->coord_[j + 1].u.y = m[6]*u.x + m[7]*u.y + m[8]*u.z; + p->coord_[j + 1].u.z = m[9]*u.x + m[10]*u.y + m[11]*u.z; + } else { + // Otherwise, copy this level's direction + p->coord_[j+1].u = p->coord_[j].u; + } + } + + // Score flux derivative accumulators for differential tallies. + if (!model::active_tallies.empty()) score_collision_derivative(p); + + if (!p->alive_ && !simulation::secondary_bank.empty()) { + revive_particle_from_secondary(p); + } + + if (p->alive_) dispatch_xs_event(p); + } + + collision_queue.clear(); +} + +void transport() +{ + int index_source = simulation::thread_work_index; + size_t remaining_work_per_thread = simulation::work_per_thread; + + while (remaining_work_per_thread > 0) { + // Initialize all histories + initialize_histories(index_source, remaining_work_per_thread); + + // Add all particles to advance particle queue + for (auto& p : particle_bank) { + dispatch_xs_event(&p); + } + + while (true) { + // Determine size of each queue + int n_fuel_xs = calculate_fuel_xs_queue.size(); + int n_nonfuel_xs = calculate_nonfuel_xs_queue.size(); + int n_advance = advance_particle_queue.size(); + int n_surface = surface_crossing_queue.size(); + int n_collision = collision_queue.size(); + //std::cout << n_xs << " " << n_advance << " " << n_surface << " " << n_collision << '\n'; + + int max = std::max({n_fuel_xs, n_nonfuel_xs, n_advance, n_surface, n_collision}); + if (max == 0) { + break; + } else if (max == n_fuel_xs) { + process_calculate_xs_events(calculate_fuel_xs_queue); + } else if (max == n_nonfuel_xs) { + process_calculate_xs_events(calculate_nonfuel_xs_queue); + } else if (max == n_advance) { + process_advance_particle_events(); + } else if (max == n_surface) { + process_surface_crossing_events(); + } else if (max == n_collision) { + process_collision_events(); + } + } + + particle_bank.clear(); + } +} + +} // namespace openmc + //============================================================================== // C API functions //============================================================================== @@ -188,6 +599,14 @@ int openmc_next_batch(int* status) // ==================================================================== // LOOP OVER PARTICLES + simulation::current_work = 1; + + #pragma omp parallel + { + transport(); + } + +/* #pragma omp parallel for schedule(runtime) for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) { simulation::current_work = i_work; @@ -199,6 +618,7 @@ int openmc_next_batch(int* status) // transport particle p.transport(); } +*/ // Accumulate time for transport simulation::time_transport.stop(); @@ -255,7 +675,9 @@ int restart_batch; bool satisfy_triggers {false}; int total_gen {0}; double total_weight; +int64_t thread_work_index; int64_t work_per_rank; +int64_t work_per_thread; const RegularMesh* entropy_mesh {nullptr}; const RegularMesh* ufs_mesh {nullptr}; @@ -498,6 +920,29 @@ void initialize_history(Particle* p, int64_t index_source) } } } + +// Display message if high verbosity or trace is on + if (settings::verbosity >= 9 || simulation::trace) { + write_message("Simulating Particle " + std::to_string(p->id_)); + } + + // // Initialize number of events to zero + // int n_event = 0; + + // Add paricle's starting weight to count for normalizing tallies later + #pragma omp atomic + simulation::total_weight += p->wgt_; + + // Force calculation of cross-sections by setting last energy to zero + if (settings::run_CE) { + for (auto& micro : p->neutron_xs_) micro.last_E = 0.0; + } + + // Prepare to write out particle track. + if (p->write_track_) add_particle_track(); + + // Every particle starts with no accumulated flux derivative. + if (!model::active_tallies.empty()) zero_flux_derivs(); } int overall_generation() @@ -528,6 +973,30 @@ void calculate_work() i_bank += work_i; simulation::work_index[i + 1] = i_bank; } + #ifdef _OPENMP + // Determine work per thread + int remainder_thread = simulation::work_per_rank % omp_get_max_threads(); + + #pragma omp parallel + { + simulation::work_per_thread = simulation::work_per_rank / omp_get_num_threads(); + if (omp_get_thread_num() < remainder_thread) { + ++simulation::work_per_thread; + } + } + + int64_t work_i = 0; + #pragma omp parallel for ordered + for (int i = 0; i < omp_get_num_threads(); ++i) { + #pragma omp ordered + { + simulation::thread_work_index = work_i; + work_i += simulation::work_per_thread; + std::cout << "Thread " << omp_get_thread_num() << ": " << simulation::work_per_thread + << std::endl; + } + } + #endif } #ifdef OPENMC_MPI From 633339293617cc63001cc413d90a3f6d627ea41b Mon Sep 17 00:00:00 2001 From: John Tramm Date: Mon, 18 Nov 2019 15:08:41 +0000 Subject: [PATCH 02/40] fixed bug by partially removing multigroup support --- src/simulation.cpp | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/src/simulation.cpp b/src/simulation.cpp index 1b566f518..d341b1c12 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -186,21 +186,16 @@ void process_calculate_xs_events(std::vector& queue) for (auto& p : queue) { // Calculate microscopic and macroscopic cross sections if (p->material_ != MATERIAL_VOID) { - if (settings::run_CE) { + // Only works for CE, no MG support + //if (settings::run_CE) { // If the material is the same as the last material and the // temperature hasn't changed, we don't need to lookup cross // sections again. model::materials[p->material_]->calculate_xs(*p); - } else { - // Get the MG data - calculate_xs_c(p->material_, p->g_, p->sqrtkT_, p->u_local(), - p->macro_xs_.total, p->macro_xs_.absorption, p->macro_xs_.nu_fission); - - // Finally, update the particle group while we have already checked - // for if multi-group - p->g_last_ = p->g_; } - } else { + //else { + //} + else { p->macro_xs_.total = 0.0; p->macro_xs_.absorption = 0.0; p->macro_xs_.fission = 0.0; From 2cc5f751b0d0120004c1091d083dd8e24a585bfa Mon Sep 17 00:00:00 2001 From: John Tramm Date: Mon, 18 Nov 2019 20:19:18 +0000 Subject: [PATCH 03/40] initial switch to SHM queues. Needs debugging. --- src/simulation.cpp | 241 +++++++++++++++++++++++++++++------------- src/tallies/tally.cpp | 1 + 2 files changed, 171 insertions(+), 71 deletions(-) diff --git a/src/simulation.cpp b/src/simulation.cpp index d341b1c12..39e6181a7 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -39,6 +39,7 @@ namespace openmc { + /* extern std::vector calculate_fuel_xs_queue; extern std::vector calculate_nonfuel_xs_queue; extern std::vector advance_particle_queue; @@ -51,24 +52,47 @@ std::vector calculate_nonfuel_xs_queue; std::vector advance_particle_queue; std::vector surface_crossing_queue; std::vector collision_queue; +*/ + +int * calculate_fuel_xs_queue; +int * calculate_nonfuel_xs_queue; +int * advance_particle_queue; +int * surface_crossing_queue; +int * collision_queue; +Particle * particles; + +int calculate_fuel_xs_queue_length = 0; +int calculate_nonfuel_xs_queue_length = 0; +int advance_particle_queue_length = 0; +int surface_crossing_queue_length = 0; +int collision_queue_length = 0; + +const int MAX_PARTICLES_IN_FLIGHT = 100; + +void init_event_queues(void) +{ + int n_particles = MAX_PARTICLES_IN_FLIGHT; + calculate_fuel_xs_queue = new int[n_particles]; + calculate_nonfuel_xs_queue = new int[n_particles]; + advance_particle_queue = new int[n_particles]; + surface_crossing_queue = new int[n_particles]; + collision_queue = new int[n_particles]; + particles = new Particle[n_particles]; +} + +void free_event_queues(void) +{ + free( calculate_fuel_xs_queue ); + free( calculate_nonfuel_xs_queue ); + free( advance_particle_queue ); + free( surface_crossing_queue ); + free( collision_queue ); + delete[] particles; +} constexpr size_t MAX_PARTICLES_PER_THREAD {100}; -void initialize_histories(int& index_source, - size_t& remaining_work_per_thread) -{ - int work = std::min(remaining_work_per_thread, MAX_PARTICLES_PER_THREAD); - for (int i = 0; i < work; ++i) { - particle_bank.emplace_back(); - auto& p {particle_bank.back()}; - - initialize_history(&p, index_source); - ++index_source; - } - - remaining_work_per_thread -= work; -} - +// TODO: What is going on here? void revive_particle_from_secondary(Particle* p) { p->from_source(&simulation::secondary_bank.back()); @@ -79,24 +103,39 @@ void revive_particle_from_secondary(Particle* p) if (p->write_track_) add_particle_track(); } -void dispatch_xs_event(Particle* p) +void dispatch_xs_event(int i) { + Particle * p = particles + i; + int idx; if (p->material_ == MATERIAL_VOID) { - calculate_nonfuel_xs_queue.push_back(p); + + #pragma omp atomic capture + idx = calculate_nonfuel_xs_queue_length++; + //std::cout << "Dispatching particle to non Fuel XS queue idx = " << idx << std::endl; + calculate_nonfuel_xs_queue[idx] = i; } else { if (model::materials[p->material_]->fissionable_) { - calculate_fuel_xs_queue.push_back(p); + #pragma omp atomic capture + idx = calculate_fuel_xs_queue_length++; + //std::cout << "Dispatching particle to Fuel XS queue idx = " << idx << std::endl; + calculate_fuel_xs_queue[idx] = i; } else { - calculate_nonfuel_xs_queue.push_back(p); + #pragma omp atomic capture + idx = calculate_nonfuel_xs_queue_length++; + //std::cout << "Dispatching particle to non Fuel XS queue idx = " << idx << std::endl; + calculate_nonfuel_xs_queue[idx] = i; } } } -void process_calculate_xs_events(std::vector& queue) +void process_calculate_xs_events(int * queue, int n) { // Save last_ members, find grid index - for (auto& p : queue) { + for (int i = 0; i < n; i++) { + Particle *p = particles + queue[i]; + //std::cout << "particle offset = " << queue[i] << std::endl; // Set the random number stream + // TODO: Move RNG seeds to particle storage if (p->type_ == Particle::Type::neutron) { prn_set_stream(STREAM_TRACKING); } else { @@ -142,7 +181,9 @@ void process_calculate_xs_events(std::vector& queue) // Calculate nuclide micros for (int i = 0; i < data::nuclides.size(); ++i) { - for (auto& p : queue) { + // loop over particles + for (int j = 0; j < n; j++) { + Particle * p = particles + queue[i]; if (p->material_ == MATERIAL_VOID) continue; // If material doesn't have this nuclide, skip it @@ -183,7 +224,8 @@ void process_calculate_xs_events(std::vector& queue) } } - for (auto& p : queue) { + for (int i = 0; i < n; i++) { + Particle * p = particles + queue[i]; // Calculate microscopic and macroscopic cross sections if (p->material_ != MATERIAL_VOID) { // Only works for CE, no MG support @@ -202,15 +244,18 @@ void process_calculate_xs_events(std::vector& queue) p->macro_xs_.nu_fission = 0.0; } - advance_particle_queue.push_back(p); + int idx; + #pragma omp atomic capture + idx = advance_particle_queue_length++; + advance_particle_queue[idx] = queue[i]; } - - queue.clear(); } void process_advance_particle_events() { - for (auto& p : advance_particle_queue) { + //for (auto& p : advance_particle_queue) { + for (int i = 0; i < advance_particle_queue_length; i++) { + Particle * p = particles + advance_particle_queue[i]; simulation::trace == (p->id_ == 0); // Sample a distance to collision @@ -231,11 +276,16 @@ void process_advance_particle_events() // Select smaller of the two distances double distance; + int idx; if (p->boundary_.distance < d_collision) { - surface_crossing_queue.push_back(p); + #pragma omp atomic capture + idx = surface_crossing_queue_length++; + surface_crossing_queue[idx] = advance_particle_queue[i]; distance = p->boundary_.distance; } else { - collision_queue.push_back(p); + #pragma omp atomic capture + idx = collision_queue_length++; + collision_queue[idx] = advance_particle_queue[i]; distance = d_collision; } @@ -265,12 +315,14 @@ void process_advance_particle_events() } } - advance_particle_queue.clear(); + advance_particle_queue_length = 0; } void process_surface_crossing_events() { - for (auto& p : surface_crossing_queue) { + //for (auto& p : surface_crossing_queue) { + for (int i = 0; i < surface_crossing_queue_length; i++) { + Particle * p = particles + surface_crossing_queue[i]; // Set surface that particle is on and adjust coordinate levels p->surface_ = p->boundary_.surface_index; p->n_coord_ = p->boundary_.coord_level; @@ -297,19 +349,27 @@ void process_surface_crossing_events() score_surface_tally(p, model::active_surface_tallies); } + // TODO: Add back in support for secondary particles + /* if (!p->alive_ && !simulation::secondary_bank.empty()) { revive_particle_from_secondary(p); } + */ - if (p->alive_) dispatch_xs_event(p); + if (p->alive_) + { + dispatch_xs_event(surface_crossing_queue[i]); + } } - surface_crossing_queue.clear(); + surface_crossing_queue_length = 0; } void process_collision_events() { - for (auto& p : collision_queue) { + //for (auto& p : collision_queue) { + for (int i = 0; i < collision_queue_length; i++) { + Particle * p = particles + collision_queue[i]; // Score collision estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE && p->type_ == Particle::Type::neutron) { @@ -379,57 +439,92 @@ void process_collision_events() // Score flux derivative accumulators for differential tallies. if (!model::active_tallies.empty()) score_collision_derivative(p); + // TODO: Add back in secondary particle support + /* if (!p->alive_ && !simulation::secondary_bank.empty()) { revive_particle_from_secondary(p); } + */ - if (p->alive_) dispatch_xs_event(p); + if (p->alive_) + { + dispatch_xs_event(collision_queue[i]); + } } - collision_queue.clear(); + collision_queue_length = 0; } void transport() { - int index_source = simulation::thread_work_index; - size_t remaining_work_per_thread = simulation::work_per_thread; + int remaining_work = simulation::work_per_rank; + int source_offset = 0; + // Subiterations to complete sets of particles + while (remaining_work > 0) { - while (remaining_work_per_thread > 0) { - // Initialize all histories - initialize_histories(index_source, remaining_work_per_thread); + // TODO: Do this only once per PI + init_event_queues(); - // Add all particles to advance particle queue - for (auto& p : particle_bank) { - dispatch_xs_event(&p); - } + // Figure out work for this subiteration + int n_particles = MAX_PARTICLES_IN_FLIGHT; + if( n_particles > remaining_work) + n_particles = remaining_work; - while (true) { - // Determine size of each queue - int n_fuel_xs = calculate_fuel_xs_queue.size(); - int n_nonfuel_xs = calculate_nonfuel_xs_queue.size(); - int n_advance = advance_particle_queue.size(); - int n_surface = surface_crossing_queue.size(); - int n_collision = collision_queue.size(); - //std::cout << n_xs << " " << n_advance << " " << n_surface << " " << n_collision << '\n'; + //std::cout << "Initializing particle histories..." << std::endl; + // Initialize all histories + // TODO: Parallelize + for (int i = 0; i < n_particles; i++) { + initialize_history(particles + i, source_offset + i); + } - int max = std::max({n_fuel_xs, n_nonfuel_xs, n_advance, n_surface, n_collision}); - if (max == 0) { - break; - } else if (max == n_fuel_xs) { - process_calculate_xs_events(calculate_fuel_xs_queue); - } else if (max == n_nonfuel_xs) { - process_calculate_xs_events(calculate_nonfuel_xs_queue); - } else if (max == n_advance) { - process_advance_particle_events(); - } else if (max == n_surface) { - process_surface_crossing_events(); - } else if (max == n_collision) { - process_collision_events(); - } - } + //std::cout << "Enqueing particles for XS Lookups..." << std::endl; + // Add all particles to advance particle queue + // TODO: Parallelize + for (int i = 0; i < n_particles; i++) { + dispatch_xs_event(i); + } - particle_bank.clear(); - } + while (true) { + std::cout << "Fuel XS Lookups = " << calculate_fuel_xs_queue_length << std::endl; + std::cout << "Non Fuel XS Lookups = " << calculate_nonfuel_xs_queue_length << std::endl; + std::cout << "Advance Particles = " << advance_particle_queue_length << std::endl; + std::cout << "Surface Crossings = " << surface_crossing_queue_length << std::endl; + std::cout << "Collisions = " << collision_queue_length << std::endl; + int max = std::max({calculate_fuel_xs_queue_length, calculate_nonfuel_xs_queue_length, advance_particle_queue_length, surface_crossing_queue_length, collision_queue_length}); + if (max == 0) { + break; + } else if (max == calculate_fuel_xs_queue_length) { + //std::cout << "Performing Fuel XS Lookups..." << std::endl; + process_calculate_xs_events(calculate_fuel_xs_queue, calculate_fuel_xs_queue_length); + //std::cout << "Done." << std::endl; + calculate_fuel_xs_queue_length = 0; + } else if (max == calculate_nonfuel_xs_queue_length) { + process_calculate_xs_events(calculate_nonfuel_xs_queue, calculate_nonfuel_xs_queue_length); + calculate_nonfuel_xs_queue_length = 0; + } else if (max == advance_particle_queue_length) { + process_advance_particle_events(); + } else if (max == surface_crossing_queue_length) { + process_surface_crossing_events(); + } else if (max == collision_queue_length) { + process_collision_events(); + } + } + + remaining_work -= n_particles; + source_offset += n_particles; + + // Should all be zero + /* + calculate_fuel_xs_queue_length = 0; + calculate_nonfuel_xs_queue_length = 0; + advance_particle_queue_length = 0; + surface_crossing_queue_length = 0; + collision_queue_length = 0; + */ + + // TODO: Do this only once per PI + free_event_queues(); + } } } // namespace openmc @@ -477,9 +572,11 @@ int openmc_simulation_init() allocate_banks(); // Allocate tally results arrays if they're not allocated yet + std::cout << "Trying to Initialize Results..." << std::endl; for (auto& t : model::tallies) { t->init_results(); } + std::cout << "Success in Initializing Results!" << std::endl; // Set up material nuclide index mapping for (auto& mat : model::materials) { @@ -596,7 +693,7 @@ int openmc_next_batch(int* status) simulation::current_work = 1; - #pragma omp parallel + //#pragma omp parallel { transport(); } @@ -987,8 +1084,10 @@ void calculate_work() { simulation::thread_work_index = work_i; work_i += simulation::work_per_thread; + /* std::cout << "Thread " << omp_get_thread_num() << ": " << simulation::work_per_thread << std::endl; + */ } } #endif diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 852918a0c..bfd91df1d 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -804,6 +804,7 @@ Tally::init_triggers(pugi::xml_node node) void Tally::init_results() { int n_scores = scores_.size() * nuclides_.size(); + std::cout << "TRYING TO ALLOCATE n filter bins, n scores = " << n_filter_bins_ << " " << n_scores << std::endl; results_ = xt::empty({n_filter_bins_, n_scores, 3}); } From 3ef3ee9fcbc4ef4f75252bac63c213605a9a8d8d Mon Sep 17 00:00:00 2001 From: John Tramm Date: Mon, 18 Nov 2019 22:04:11 +0000 Subject: [PATCH 04/40] fixed the bug -- was using wrong index in one of the event loops (should have been j, but was using i). Now to remove all the debugging print statements --- src/distribution_angle.cpp | 1 + src/physics.cpp | 18 +++++++ src/simulation.cpp | 97 ++++++++++++++++++++++++++++++++++++-- 3 files changed, 113 insertions(+), 3 deletions(-) diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp index 615f84f68..f12c00286 100644 --- a/src/distribution_angle.cpp +++ b/src/distribution_angle.cpp @@ -85,6 +85,7 @@ double AngleDistribution::sample(double E) const // Sample between the ith and (i+1)th bin if (r > prn()) ++i; + //assert(i >= 0 ); // Sample i-th distribution double mu = distribution_[i]->sample(); diff --git a/src/physics.cpp b/src/physics.cpp index 0eb8c494e..7dae2431c 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -79,6 +79,7 @@ void collision(Particle* p) void sample_neutron_reaction(Particle* p) { + std::cout << "particle energy e = " << p->E_ << std::endl; // Sample a nuclide within the material int i_nuclide = sample_nuclide(p); @@ -127,9 +128,12 @@ void sample_neutron_reaction(Particle* p) } if (!p->alive_) return; + std::cout << "before scatter particle energy e = " << p->E_ << " with nuclide = " << i_nuclide << std::endl; // Sample a scattering reaction and determine the secondary energy of the // exiting neutron scatter(p, i_nuclide); + + std::cout << "after scatter particle energy e = " << p->E_ << std::endl; // Advance URR seed stream 'N' times after energy changes if (p->E_ != p->E_last_) { @@ -598,6 +602,7 @@ void scatter(Particle* p, int i_nuclide) // Determine temperature double kT = nuc->multipole_ ? p->sqrtkT_*p->sqrtkT_ : nuc->kTs_[i_temp]; + std::cout << "we are doing an elastic scatter" << std::endl; // Perform collision physics for elastic scattering elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p); @@ -610,6 +615,7 @@ void scatter(Particle* p, int i_nuclide) // ======================================================================= // S(A,B) SCATTERING + std::cout << "we are doing an SAB scatter" << std::endl; sab_scatter(i_nuclide, micro.index_sab, p); p->event_mt_ = ELASTIC; @@ -617,6 +623,7 @@ void scatter(Particle* p, int i_nuclide) } if (!sampled) { + std::cout << "we are doing an Inelastic scatter" << std::endl; // ======================================================================= // INELASTIC SCATTERING @@ -643,7 +650,9 @@ void scatter(Particle* p, int i_nuclide) // Perform collision physics for inelastic scattering const auto& rx {nuc->reactions_[i]}; + std::cout << "Energy before Inelastic scatter = " << p->E_ << std::endl; inelastic_scatter(nuc.get(), rx.get(), p); + std::cout << "Energy after Inelastic scatter = " << p->E_ << std::endl; p->event_mt_ = rx->mt_; } @@ -1050,28 +1059,37 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p) // sample outgoing energy and scattering cosine double E; double mu; + // This gives a negative sometimes for some unknown reason rx->products_[0].sample(E_in, E, mu); +std::cout << "E_in = " << E_in << "E = " << E << " mu = " << mu << std::endl; // if scattering system is in center-of-mass, transfer cosine of scattering // angle and outgoing energy from CM to LAB if (rx->scatter_in_cm_) { double E_cm = E; + std::cout << "E_cm = " << E_cm << std::endl; // determine outgoing energy in lab double A = nuc->awr_; + std::cout << "A = " << nuc->awr_ << std::endl; E = E_cm + (E_in + 2.0*mu*(A + 1.0) * std::sqrt(E_in*E_cm)) / ((A + 1.0)*(A + 1.0)); + // here's where it goes wrong + std::cout << "E = " << E << std::endl; // determine outgoing angle in lab mu = mu*std::sqrt(E_cm/E) + 1.0/(A+1.0) * std::sqrt(E_in/E); + std::cout << "mu = " << mu << std::endl; } // Because of floating-point roundoff, it may be possible for mu to be // outside of the range [-1,1). In these cases, we just set mu to exactly -1 // or 1 if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu); + std::cout << "mu = " << mu << std::endl; // Set outgoing energy and scattering angle + std::cout << "FINAL E = " << E << " mu = " << mu << std::endl; p->E_ = E; p->mu_ = mu; diff --git a/src/simulation.cpp b/src/simulation.cpp index 39e6181a7..f56dfb2ef 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -183,7 +183,8 @@ void process_calculate_xs_events(int * queue, int n) for (int i = 0; i < data::nuclides.size(); ++i) { // loop over particles for (int j = 0; j < n; j++) { - Particle * p = particles + queue[i]; + //Particle * p = particles + queue[i]; + Particle * p = particles + queue[j]; if (p->material_ == MATERIAL_VOID) continue; // If material doesn't have this nuclide, skip it @@ -370,6 +371,7 @@ void process_collision_events() //for (auto& p : collision_queue) { for (int i = 0; i < collision_queue_length; i++) { Particle * p = particles + collision_queue[i]; + std::cout << "Beginning collision of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl; // Score collision estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE && p->type_ == Particle::Type::neutron) { @@ -383,6 +385,8 @@ void process_collision_events() if (!model::active_meshsurf_tallies.empty()) score_surface_tally(p, model::active_meshsurf_tallies); + + std::cout << "After surface tally of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl; // Clear surface component p->surface_ = 0; @@ -393,6 +397,8 @@ void process_collision_events() collision_mg(p); } + std::cout << "After collision() of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl; + // Score collision estimator tallies -- this is done after a collision // has occurred rather than before because we need information on the // outgoing energy for any tallies with an outgoing energy filter @@ -404,6 +410,7 @@ void process_collision_events() score_analog_tally_mg(p); } } + std::cout << "After analog tally of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl; // Reset banked weight during collision p->n_bank_ = 0; @@ -449,12 +456,82 @@ void process_collision_events() if (p->alive_) { dispatch_xs_event(collision_queue[i]); + std::cout << "Ended collision of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl; + assert(std::isfinite(p->E_) ); } } collision_queue_length = 0; } +void check_energies(void) +{ + int * Q; + int n; + + + Q = calculate_fuel_xs_queue; + n = calculate_fuel_xs_queue_length; + for( int i = 0; i < n; i++ ) + { + Particle * p = particles + Q[i]; + Particle p_true = *p; + if( !std::isfinite(particles[Q[i]].E_ ) ) + { + std::cout << "NAN energy particle found at index xs FUEL " << Q[i] << std::endl; + assert(0); + } + } + Q = calculate_nonfuel_xs_queue; + n = calculate_nonfuel_xs_queue_length ; + for( int i = 0; i < n; i++ ) + { + Particle * p = particles + Q[i]; + Particle p_true = *p; + if( !std::isfinite(particles[Q[i]].E_ ) ) + { + std::cout << "NAN energy particle found at index xs Non fuel " << Q[i] << std::endl; + assert(0); + } + } + Q = advance_particle_queue; + n = advance_particle_queue_length ; + for( int i = 0; i < n; i++ ) + { + Particle * p = particles + Q[i]; + Particle p_true = *p; + if( !std::isfinite(particles[Q[i]].E_ ) ) + { + std::cout << "NAN energy particle found at index advance particle " << Q[i] << std::endl; + assert(0); + } + } + Q = surface_crossing_queue; + n = surface_crossing_queue_length ; + for( int i = 0; i < n; i++ ) + { + Particle * p = particles + Q[i]; + Particle p_true = *p; + if( !std::isfinite(particles[Q[i]].E_ ) ) + { + std::cout << "NAN energy particle found at index surface crossing " << Q[i] << std::endl; + assert(0); + } + } + Q = collision_queue; + n = collision_queue_length ; + for( int i = 0; i < n; i++ ) + { + Particle * p = particles + Q[i]; + Particle p_true = *p; + if( !std::isfinite(particles[Q[i]].E_ ) ) + { + std::cout << "NAN energy particle found at index collision " << Q[i] << std::endl; + assert(0); + } + } +} + void transport() { int remaining_work = simulation::work_per_rank; @@ -491,21 +568,35 @@ void transport() std::cout << "Surface Crossings = " << surface_crossing_queue_length << std::endl; std::cout << "Collisions = " << collision_queue_length << std::endl; int max = std::max({calculate_fuel_xs_queue_length, calculate_nonfuel_xs_queue_length, advance_particle_queue_length, surface_crossing_queue_length, collision_queue_length}); + check_energies(); if (max == 0) { break; } else if (max == calculate_fuel_xs_queue_length) { - //std::cout << "Performing Fuel XS Lookups..." << std::endl; + std::cout << "pre fuel XS check..." << std::endl; + //check_energies(calculate_fuel_xs_queue, calculate_fuel_xs_queue_length); + std::cout << "Performing Fuel XS Lookups..." << std::endl; process_calculate_xs_events(calculate_fuel_xs_queue, calculate_fuel_xs_queue_length); - //std::cout << "Done." << std::endl; calculate_fuel_xs_queue_length = 0; } else if (max == calculate_nonfuel_xs_queue_length) { + std::cout << "pre non fuel XS check..." << std::endl; + //check_energies(calculate_nonfuel_xs_queue, calculate_nonfuel_xs_queue_length); + std::cout << "Performing Non Fuel XS Lookups..." << std::endl; process_calculate_xs_events(calculate_nonfuel_xs_queue, calculate_nonfuel_xs_queue_length); calculate_nonfuel_xs_queue_length = 0; } else if (max == advance_particle_queue_length) { + std::cout << "pre advancing check..." << std::endl; + //check_energies(advance_particle_queue, advance_particle_queue_length); + std::cout << "Advancing Particles..." << std::endl; process_advance_particle_events(); } else if (max == surface_crossing_queue_length) { + std::cout << "pre surface crossing check..." << std::endl; + //check_energies(surface_crossing_queue, surface_crossing_queue_length); + std::cout << "Surface Crossings..." << std::endl; process_surface_crossing_events(); } else if (max == collision_queue_length) { + std::cout << "pre Colliding check..." << std::endl; + //check_energies(collision_queue, collision_queue_length); + std::cout << "Colliding..." << std::endl; process_collision_events(); } } From 2f9dba6a5813505a5de60dbd0b735628632815ef Mon Sep 17 00:00:00 2001 From: John Tramm Date: Mon, 18 Nov 2019 22:28:19 +0000 Subject: [PATCH 05/40] removed print statements. Works in serial as expected. --- examples/xml/pincell/settings.xml | 4 ++-- src/physics.cpp | 30 +++++++++++++------------- src/simulation.cpp | 35 ++++++++++++++++++------------- 3 files changed, 37 insertions(+), 32 deletions(-) diff --git a/examples/xml/pincell/settings.xml b/examples/xml/pincell/settings.xml index 5845898c8..7b0a985f1 100644 --- a/examples/xml/pincell/settings.xml +++ b/examples/xml/pincell/settings.xml @@ -3,8 +3,8 @@ eigenvalue - 100 - 10 + 10 + 5 1000 eigenvalue - 10 - 5 + 100 + 10 1000