From 0fb047f56cabe7b95f7ee2ae5a7d36496f21ab42 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 15 Nov 2019 21:28:24 +0000 Subject: [PATCH] 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