From 6f80abfe2bdc430914a7933d0b364255e3378a3f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 9 Oct 2018 22:34:08 -0500 Subject: [PATCH 01/23] Move variables in simulation.h into simulation namespace --- include/openmc/capi.h | 3 --- include/openmc/simulation.h | 25 +++++++++++++++++-------- openmc/capi/core.py | 2 +- src/cmfd_execute.cpp | 3 ++- src/eigenvalue.cpp | 6 +++--- src/geometry.cpp | 4 ++-- src/initialize.cpp | 7 ++++--- src/particle.cpp | 16 ++++++++-------- src/settings.cpp | 17 +++++++++-------- src/simulation.cpp | 12 ++++++++---- src/simulation_header.F90 | 19 +++++++++---------- src/source.cpp | 16 ++++++++-------- src/state_point.cpp | 9 +++++---- 13 files changed, 76 insertions(+), 63 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 57d817d23..3bf1d83ce 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -148,7 +148,6 @@ extern "C" { extern int32_t n_surfaces; extern int32_t n_tallies; extern int32_t n_universes; - extern bool openmc_simulation_initialized; extern double global_tally_absorption; #pragma omp threadprivate(global_tally_absorption) @@ -156,9 +155,7 @@ extern "C" { // later) extern bool openmc_master; extern int openmc_n_procs; - extern int openmc_n_threads; extern int openmc_rank; - extern int64_t openmc_work; #ifdef __cplusplus } diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 9c0779072..939b26c14 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -10,19 +10,28 @@ namespace openmc { //============================================================================== -// Global variables +// Global variable declarations //============================================================================== -extern "C" int openmc_current_batch; -extern "C" int openmc_current_gen; -extern "C" int64_t openmc_current_work; -extern "C" int openmc_n_lost_particles; -extern "C" int openmc_total_gen; -extern "C" bool openmc_trace; +namespace simulation { + +extern "C" int current_batch; +extern "C" int current_gen; +extern "C" int64_t current_work; +extern "C" int n_lost_particles; +extern "C" int total_gen; +extern "C" bool trace; +extern "C" int64_t work; + +#ifdef _OPENMP +extern "C" int n_threads; +#endif extern std::vector work_index; -#pragma omp threadprivate(openmc_current_work, openmc_trace) +#pragma omp threadprivate(current_work, trace) + +} // namespace simulation //============================================================================== // Functions diff --git a/openmc/capi/core.py b/openmc/capi/core.py index fa6517488..205123ff8 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -71,7 +71,7 @@ def current_batch(): Current batch of the simulation """ - return c_int.in_dll(_dll, 'openmc_current_batch').value + return c_int.in_dll(_dll, 'current_batch').value def finalize(): diff --git a/src/cmfd_execute.cpp b/src/cmfd_execute.cpp index 740c1f51f..00d162069 100644 --- a/src/cmfd_execute.cpp +++ b/src/cmfd_execute.cpp @@ -7,6 +7,7 @@ #include "openmc/capi.h" #include "openmc/mesh.h" +#include "openmc/simulation.h" namespace openmc { @@ -23,7 +24,7 @@ cmfd_populate_sourcecounts(int n_energy, const double* energies, // Get source counts in each mesh bin / energy bin auto& m = meshes.at(index_cmfd_mesh); - xt::xarray counts = m->count_sites(openmc_work, source_bank, n_energy, energies, outside); + xt::xarray counts = m->count_sites(simulation::work, source_bank, n_energy, energies, outside); // Copy data from the xarray into the source counts array std::copy(counts.begin(), counts.end(), source_counts); diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 454fd57e6..f03820492 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -66,7 +66,7 @@ void ufs_count_sites() { auto &m = meshes[settings::index_ufs_mesh]; - if (openmc_current_batch == 1 && openmc_current_gen == 1) { + if (simulation::current_batch == 1 && simulation::current_gen == 1) { // On the first generation, just assume that the source is already evenly // distributed so that effectively the production of fission sites is not // biased @@ -82,7 +82,7 @@ void ufs_count_sites() // count number of source sites in each ufs mesh cell bool sites_outside; - source_frac = m->count_sites(openmc_work, source_bank, 0, nullptr, + source_frac = m->count_sites(simulation::work, source_bank, 0, nullptr, &sites_outside); // Check for sites outside of the mesh @@ -102,7 +102,7 @@ void ufs_count_sites() // Since the total starting weight is not equal to n_particles, we need to // renormalize the weight of the source sites - for (int i = 0; i < openmc_work; ++i) { + for (int i = 0; i < simulation::work; ++i) { source_bank[i].wgt *= settings::n_particles / total; } } diff --git a/src/geometry.cpp b/src/geometry.cpp index 857dbd32e..28885a5db 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -91,7 +91,7 @@ find_cell(Particle* p, int search_surf) { if (cells[i_cell]->contains(r, u, surf)) { p->coord[p->n_coord-1].cell = i_cell; - if (settings::verbosity >= 10 || openmc_trace) { + if (settings::verbosity >= 10 || simulation::trace) { std::stringstream msg; msg << " Entering cell " << cells[i_cell]->id_; write_message(msg, 1); @@ -255,7 +255,7 @@ cross_lattice(Particle* p, int lattice_translation[3]) { Lattice& lat {*lattices[p->coord[p->n_coord-1].lattice-1]}; - if (settings::verbosity >= 10 || openmc_trace) { + if (settings::verbosity >= 10 || simulation::trace) { std::stringstream msg; msg << " Crossing lattice " << lat.id_ << ". Current position (" << p->coord[p->n_coord-1].lattice_x << "," diff --git a/src/initialize.cpp b/src/initialize.cpp index 0c26135dc..ac7a5ef15 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -15,6 +15,7 @@ #include "openmc/hdf5_interface.h" #include "openmc/message_passing.h" #include "openmc/settings.h" +#include "openmc/simulation.h" #include "openmc/string_utils.h" // data/functions from Fortran side @@ -172,13 +173,13 @@ parse_command_line(int argc, char* argv[]) #ifdef _OPENMP // Read and set number of OpenMP threads - openmc_n_threads = std::stoi(argv[i]); - if (openmc_n_threads < 1) { + simulation::n_threads = std::stoi(argv[i]); + if (simulation::n_threads < 1) { std::string msg {"Number of threads must be positive."}; strcpy(openmc_err_msg, msg.c_str()); return OPENMC_E_INVALID_ARGUMENT; } - omp_set_num_threads(openmc_n_threads); + omp_set_num_threads(simulation::n_threads); #else if (openmc_master) warning("Ignoring number of threads specified on command line."); diff --git a/src/particle.cpp b/src/particle.cpp index ca57c2faa..d6dfa969c 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -127,15 +127,15 @@ Particle::mark_as_lost(const char* message) // Increment number of lost particles alive = false; #pragma omp atomic - openmc_n_lost_particles += 1; + simulation::n_lost_particles += 1; // Count the total number of simulated particles (on this processor) - auto n = openmc_current_batch * settings::gen_per_batch * openmc_work; + auto n = simulation::current_batch * settings::gen_per_batch * simulation::work; // Abort the simulation if the maximum number of lost particles has been // reached - if (openmc_n_lost_particles >= MAX_LOST_PARTICLES && - openmc_n_lost_particles >= REL_MAX_LOST_PARTICLES*n) { + if (simulation::n_lost_particles >= MAX_LOST_PARTICLES && + simulation::n_lost_particles >= REL_MAX_LOST_PARTICLES*n) { fatal_error("Maximum number of lost particles has been reached."); } } @@ -148,7 +148,7 @@ Particle::write_restart() const // Set up file name std::stringstream filename; - filename << settings::path_output << "particle_" << openmc_current_batch + filename << settings::path_output << "particle_" << simulation::current_batch << '_' << id << ".h5"; #pragma omp critical (WriteParticleRestart) @@ -165,9 +165,9 @@ Particle::write_restart() const #endif // Write data to file - write_dataset(file_id, "current_batch", openmc_current_batch); + write_dataset(file_id, "current_batch", simulation::current_batch); write_dataset(file_id, "generations_per_batch", settings::gen_per_batch); - write_dataset(file_id, "current_generation", openmc_current_gen); + write_dataset(file_id, "current_generation", simulation::current_gen); write_dataset(file_id, "n_particles", settings::n_particles); switch (settings::run_mode) { case RUN_MODE_FIXEDSOURCE: @@ -188,7 +188,7 @@ Particle::write_restart() const int64_t n; openmc_source_bank(&src, &n); - int64_t i = openmc_current_work; + int64_t i = simulation::current_work; write_dataset(file_id, "weight", src[i-1].wgt); write_dataset(file_id, "energy", src[i-1].E); hsize_t dims[] {3}; diff --git a/src/settings.cpp b/src/settings.cpp index b5f2bf5c6..bf97e44b8 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -17,6 +17,7 @@ #include "openmc/mesh.h" #include "openmc/output.h" #include "openmc/random_lcg.h" +#include "openmc/simulation.h" #include "openmc/source.h" #include "openmc/string_utils.h" #include "openmc/xml_interface.h" @@ -57,7 +58,7 @@ bool urr_ptables_on {true}; bool write_all_tracks {false}; bool write_initial_source {false}; bool dagmc {false}; - + std::string path_cross_sections; std::string path_input; std::string path_multipole; @@ -218,7 +219,7 @@ void read_settings_xml() fatal_error("DAGMC mode unsupported for this build of OpenMC"); } #endif - + // To this point, we haven't displayed any output since we didn't know what // the verbosity is. Now that we checked for it, show the title if necessary if (openmc_master) { @@ -393,14 +394,14 @@ void read_settings_xml() // Number of OpenMP threads if (check_for_node(root, "threads")) { #ifdef _OPENMP - if (openmc_n_threads == 0) { - openmc_n_threads = std::stoi(get_node_value(root, "threads")); - if (openmc_n_threads < 1) { + if (simulation::n_threads == 0) { + simulation::n_threads = std::stoi(get_node_value(root, "threads")); + if (simulation::n_threads < 1) { std::stringstream msg; - msg << "Invalid number of threads: " << openmc_n_threads; + msg << "Invalid number of threads: " << simulation::n_threads; fatal_error(msg); } - omp_set_num_threads(openmc_n_threads); + omp_set_num_threads(simulation::n_threads); } #else if (openmc_master) warning("OpenMC was not compiled with OpenMP support; " @@ -414,7 +415,7 @@ void read_settings_xml() omp_set_num_threads(1); } #endif - + // ========================================================================== // EXTERNAL SOURCE diff --git a/src/simulation.cpp b/src/simulation.cpp index 74a3c14fc..5c2a13678 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -27,8 +27,12 @@ namespace openmc { // Global variables //============================================================================== +namespace simulation { + std::vector work_index; +} // namespace simulation + //============================================================================== // Functions //============================================================================== @@ -48,18 +52,18 @@ void calculate_work() int64_t remainder = settings::n_particles % mpi::n_procs; int64_t i_bank = 0; - work_index.reserve(mpi::n_procs); - work_index.push_back(0); + simulation::work_index.reserve(mpi::n_procs); + simulation::work_index.push_back(0); for (int i = 0; i < mpi::n_procs; ++i) { // Number of particles for rank i int64_t work_i = i < remainder ? min_work + 1 : min_work; // Set number of particles - if (mpi::rank == i) openmc_work = work_i; + if (mpi::rank == i) simulation::work = work_i; // Set index into source bank for rank i i_bank += work_i; - work_index.push_back(i_bank); + simulation::work_index.push_back(i_bank); } } diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index ad5513c61..b96caf257 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -13,18 +13,17 @@ module simulation_header ! GEOMETRY-RELATED VARIABLES ! Number of lost particles - integer(C_INT), bind(C, name='openmc_n_lost_particles') :: n_lost_particles = 0 + integer(C_INT), bind(C) :: n_lost_particles = 0 real(8) :: log_spacing ! spacing on logarithmic grid ! ============================================================================ ! SIMULATION VARIABLES - integer(C_INT), bind(C, name='openmc_current_batch') :: current_batch ! current batch - integer(C_INT), bind(C, name='openmc_current_gen') :: current_gen ! current generation within a batch - integer(C_INT), bind(C, name='openmc_total_gen') :: total_gen = 0 ! total number of generations simulated - logical(C_BOOL), bind(C, name='openmc_simulation_initialized') :: & - simulation_initialized = .false. + integer(C_INT), bind(C) :: current_batch ! current batch + integer(C_INT), bind(C) :: current_gen ! current generation within a batch + integer(C_INT), bind(C) :: total_gen = 0 ! total number of generations simulated + logical(C_BOOL), bind(C) :: simulation_initialized = .false. logical :: need_depletion_rx ! need to calculate depletion reaction rx? ! ============================================================================ @@ -32,9 +31,9 @@ module simulation_header logical :: satisfy_triggers = .false. ! whether triggers are satisfied - integer(C_INT64_T), bind(C, name='openmc_work') :: work ! number of particles per processor + integer(C_INT64_T), bind(C) :: work ! number of particles per processor integer(C_INT64_T), allocatable :: work_index(:) ! starting index in source bank for each process - integer(C_INT64_T), bind(C, name='openmc_current_work') :: current_work ! index in source bank of current history simulated + integer(C_INT64_T), bind(C) :: current_work ! index in source bank of current history simulated ! ============================================================================ ! K-EIGENVALUE SIMULATION VARIABLES @@ -51,7 +50,7 @@ module simulation_header ! PARALLEL PROCESSING VARIABLES #ifdef _OPENMP - integer(C_INT), bind(C, name='openmc_n_threads') :: n_threads = NONE ! number of OpenMP threads + integer(C_INT), bind(C) :: n_threads = NONE ! number of OpenMP threads integer :: thread_id ! ID of a given thread #endif @@ -60,7 +59,7 @@ module simulation_header integer :: restart_batch - logical(C_BOOL), bind(C, name='openmc_trace') :: trace + logical(C_BOOL), bind(C) :: trace !$omp threadprivate(trace, thread_id, current_work) diff --git a/src/source.cpp b/src/source.cpp index 666bbb783..1290c5f0d 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -264,17 +264,17 @@ void initialize_source() } // Read in the source bank - read_source_bank(file_id, work_index.data(), source_bank); + read_source_bank(file_id, simulation::work_index.data(), source_bank); // Close file file_close(file_id); } else { // Generation source sites from specified distribution in user input - for (int64_t i = 0; i < openmc_work; ++i) { + for (int64_t i = 0; i < simulation::work; ++i) { // initialize random number seed - int64_t id = openmc_total_gen*settings::n_particles + - work_index[openmc::mpi::rank] + i + 1; + int64_t id = simulation::total_gen*settings::n_particles + + simulation::work_index[mpi::rank] + i + 1; set_particle_seed(id); // sample external source distribution @@ -287,7 +287,7 @@ void initialize_source() write_message("Writing out initial source...", 5); std::string filename = settings::path_output + "initial_source.h5"; hid_t file_id = file_open(filename, 'w', true); - write_source_bank(file_id, work_index.data(), source_bank); + write_source_bank(file_id, simulation::work_index.data(), source_bank); file_close(file_id); } } @@ -364,10 +364,10 @@ extern "C" void fill_source_bank_fixedsource() int64_t n; openmc_source_bank(&source_bank, &n); - for (int64_t i = 0; i < openmc_work; ++i) { + for (int64_t i = 0; i < simulation::work; ++i) { // initialize random number seed - int64_t id = (openmc_total_gen + overall_generation()) * - settings::n_particles + work_index[openmc::mpi::rank] + i + 1; + int64_t id = (simulation::total_gen + overall_generation()) * + settings::n_particles + simulation::work_index[mpi::rank] + i + 1; set_particle_seed(id); // sample external source distribution diff --git a/src/state_point.cpp b/src/state_point.cpp index dba0fe1b8..e24993b6c 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -11,6 +11,7 @@ #include "openmc/error.h" #include "openmc/message_passing.h" #include "openmc/settings.h" +#include "openmc/simulation.h" namespace openmc { @@ -46,7 +47,7 @@ write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); // Create another data space but for each proc individually - hsize_t count[] {static_cast(openmc_work)}; + hsize_t count[] {static_cast(simulation::work)}; hid_t memspace = H5Screate_simple(1, count, nullptr); // Select hyperslab for this dataspace @@ -77,7 +78,7 @@ write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) // Save source bank sites since the souce_bank array is overwritten below #ifdef OPENMC_MPI - std::vector temp_source {source_bank, source_bank + openmc_work}; + std::vector temp_source {source_bank, source_bank + simulation::work}; #endif for (int i = 0; i < openmc::mpi::n_procs; ++i) { @@ -113,7 +114,7 @@ write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) #endif } else { #ifdef OPENMC_MPI - MPI_Send(source_bank, openmc_work, openmc::mpi::bank, 0, openmc::mpi::rank, + MPI_Send(source_bank, simulation::work, openmc::mpi::bank, 0, openmc::mpi::rank, openmc::mpi::intracomm); #endif } @@ -131,7 +132,7 @@ void read_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) hid_t dset = H5Dopen(group_id, "source_bank", H5P_DEFAULT); // Create another data space but for each proc individually - hsize_t dims[] {static_cast(openmc_work)}; + hsize_t dims[] {static_cast(simulation::work)}; hid_t memspace = H5Screate_simple(1, dims, nullptr); // Make sure source bank is big enough From 388c687fd2c08a6754fb1b38ed74ff0ab849e762 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 9 Oct 2018 22:42:21 -0500 Subject: [PATCH 02/23] Make remaining variables in simulation_header.F90 C interoperable --- src/simulation_header.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index b96caf257..2e35749b9 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -15,7 +15,7 @@ module simulation_header ! Number of lost particles integer(C_INT), bind(C) :: n_lost_particles = 0 - real(8) :: log_spacing ! spacing on logarithmic grid + real(C_DOUBLE), bind(C) :: log_spacing ! spacing on logarithmic grid ! ============================================================================ ! SIMULATION VARIABLES @@ -24,12 +24,12 @@ module simulation_header integer(C_INT), bind(C) :: current_gen ! current generation within a batch integer(C_INT), bind(C) :: total_gen = 0 ! total number of generations simulated logical(C_BOOL), bind(C) :: simulation_initialized = .false. - logical :: need_depletion_rx ! need to calculate depletion reaction rx? + logical(C_BOOL), bind(C) :: need_depletion_rx ! need to calculate depletion reaction rx? ! ============================================================================ ! TALLY PRECISION TRIGGER VARIABLES - logical :: satisfy_triggers = .false. ! whether triggers are satisfied + logical(C_BOOL), bind(C) :: satisfy_triggers = .false. ! whether triggers are satisfied integer(C_INT64_T), bind(C) :: work ! number of particles per processor integer(C_INT64_T), allocatable :: work_index(:) ! starting index in source bank for each process @@ -42,22 +42,22 @@ module simulation_header type(VectorReal) :: k_generation ! single-generation estimates of k real(C_DOUBLE), bind(C, name='openmc_keff') :: keff = ONE ! average k over active batches real(C_DOUBLE), bind(C, name='openmc_keff_std') :: keff_std ! standard deviation of average k - real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption - real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength - real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength + real(C_DOUBLE), bind(C) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption + real(C_DOUBLE), bind(C) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength + real(C_DOUBLE), bind(C) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength ! ============================================================================ ! PARALLEL PROCESSING VARIABLES #ifdef _OPENMP integer(C_INT), bind(C) :: n_threads = NONE ! number of OpenMP threads - integer :: thread_id ! ID of a given thread + integer(C_INT), bind(C) :: thread_id ! ID of a given thread #endif ! ============================================================================ ! MISCELLANEOUS VARIABLES - integer :: restart_batch + integer(C_INT), bind(C) :: restart_batch logical(C_BOOL), bind(C) :: trace From 3b85a657111275132a3fb8393999a2682ef99011 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 10 Oct 2018 07:02:07 -0500 Subject: [PATCH 03/23] Add remaining variables in simulation.h --- include/openmc/capi.h | 2 -- include/openmc/simulation.h | 36 ++++++++++++++++++++++++------------ openmc/capi/core.py | 4 ++-- src/physics_mg.cpp | 4 ++-- src/simulation_header.F90 | 4 ++-- 5 files changed, 30 insertions(+), 20 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 3bf1d83ce..8e9edc4a0 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -134,8 +134,6 @@ extern "C" { // Global variables extern char openmc_err_msg[256]; - extern double openmc_keff; - extern double openmc_keff_std; extern int32_t n_cells; extern int32_t n_filters; extern int32_t n_lattices; diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 939b26c14..38b0ad314 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -15,21 +15,33 @@ namespace openmc { namespace simulation { -extern "C" int current_batch; -extern "C" int current_gen; -extern "C" int64_t current_work; -extern "C" int n_lost_particles; -extern "C" int total_gen; -extern "C" bool trace; -extern "C" int64_t work; - -#ifdef _OPENMP -extern "C" int n_threads; -#endif +extern "C" int current_batch; //!< current batch +extern "C" int current_gen; //!< current fission generation +extern "C" int64_t current_work; //!< index in source back of current particle +extern "C" double keff; //!< average k over batches +extern "C" double keff_std; //!< standard deviation of average k +extern "C" double k_col_abs; //!< sum over batches of k_collision * k_absorption +extern "C" double k_col_tra; //!< sum over batches of k_collision * k_tracklength +extern "C" double k_abs_tra; //!< sum over batches of k_absorption * k_tracklength +extern "C" double log_spacing; //!< lethargy spacing for energy grid searches +extern "C" int n_lost_particles; //!< cumulative number of lost particles +extern "C" bool need_depletion_rx; //!< need to calculate depletion rx? +extern "C" bool restart_batch; //!< batch at which a restart job resumed +extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied? +extern "C" bool simulation_initialized; //!< has simulation been initialized? +extern "C" int total_gen; //!< total number of generations simulated +extern "C" int64_t work; //!< number of particles per process extern std::vector work_index; -#pragma omp threadprivate(current_work, trace) +// Threadprivate variables +extern "C" bool trace; //!< flag to show debug information +#ifdef _OPENMP +extern "C" int n_threads; //!< number of OpenMP threads +extern "C" int thread_id; //!< ID of a given thread +#endif + +#pragma omp threadprivate(current_work, thread_id, trace) } // namespace simulation diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 205123ff8..5be31ce86 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -220,8 +220,8 @@ def keff(): return tuple(k) else: # Otherwise, return the tracklength estimator - mean = c_double.in_dll(_dll, 'openmc_keff').value - std_dev = c_double.in_dll(_dll, 'openmc_keff_std').value \ + mean = c_double.in_dll(_dll, 'keff').value + std_dev = c_double.in_dll(_dll, 'keff_std').value \ if n > 1 else np.inf return (mean, std_dev) diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 38eb9a5b9..18301c3ab 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -30,7 +30,7 @@ collision_mg(Particle* p, const double* energy_bin_avg, sample_reaction(p, energy_bin_avg, material_xs); // Display information about collision - if ((settings::verbosity >= 10) || (openmc_trace)) { + if ((settings::verbosity >= 10) || (simulation::trace)) { std::stringstream msg; msg << " Energy Group = " << p->g; write_message(msg, 1); @@ -115,7 +115,7 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; // Determine the expected number of neutrons produced - double nu_t = p->wgt / openmc_keff * weight * material_xs->nu_fission / + double nu_t = p->wgt / simulation::keff * weight * material_xs->nu_fission / material_xs->total; // Sample the number of neutrons produced diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 2e35749b9..8328b13b1 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -40,8 +40,8 @@ module simulation_header ! Temporary k-effective values type(VectorReal) :: k_generation ! single-generation estimates of k - real(C_DOUBLE), bind(C, name='openmc_keff') :: keff = ONE ! average k over active batches - real(C_DOUBLE), bind(C, name='openmc_keff_std') :: keff_std ! standard deviation of average k + real(C_DOUBLE), bind(C) :: keff = ONE ! average k over active batches + real(C_DOUBLE), bind(C) :: keff_std ! standard deviation of average k real(C_DOUBLE), bind(C) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption real(C_DOUBLE), bind(C) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength real(C_DOUBLE), bind(C) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength From d3b9d65fc7a6b9c6a8af1949b3eb7eea6eecf676 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 10 Oct 2018 07:17:33 -0500 Subject: [PATCH 04/23] Add ifdef around variable in tracking.F90 --- src/tracking.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/tracking.F90 b/src/tracking.F90 index 94c009be0..dc0c008d6 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -322,7 +322,9 @@ contains real(8) :: norm ! "norm" of surface normal real(8) :: xyz(3) ! Saved global coordinate integer :: i_surface ! index in surfaces +#ifdef DAGMC integer :: i_cell ! index of new cell +#endif logical :: rotational ! if rotational periodic BC applied logical :: found ! particle found in universe? class(Surface), pointer :: surf From 3ba886cd6ddc159d4d071c74480c8a140adbd931 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 10 Oct 2018 07:17:52 -0500 Subject: [PATCH 05/23] Move variable definitions from simulation_header.F90 to simulation.cpp --- include/openmc/simulation.h | 6 +++++- src/simulation.cpp | 31 +++++++++++++++++++++++++++++++ src/simulation_header.F90 | 37 ++++++++++++++++--------------------- 3 files changed, 52 insertions(+), 22 deletions(-) diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 38b0ad314..564516e60 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -26,12 +26,13 @@ extern "C" double k_abs_tra; //!< sum over batches of k_absorption * k_track extern "C" double log_spacing; //!< lethargy spacing for energy grid searches extern "C" int n_lost_particles; //!< cumulative number of lost particles extern "C" bool need_depletion_rx; //!< need to calculate depletion rx? -extern "C" bool restart_batch; //!< batch at which a restart job resumed +extern "C" int restart_batch; //!< batch at which a restart job resumed extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied? extern "C" bool simulation_initialized; //!< has simulation been initialized? extern "C" int total_gen; //!< total number of generations simulated extern "C" int64_t work; //!< number of particles per process +extern std::vector k_generation; extern std::vector work_index; // Threadprivate variables @@ -52,6 +53,9 @@ extern "C" int thread_id; //!< ID of a given thread //! Initialize simulation extern "C" void openmc_simulation_init_c(); +//! Determine overall generation number +extern "C" int overall_generation(); + //! Determine number of particles to transport per process void calculate_work(); diff --git a/src/simulation.cpp b/src/simulation.cpp index 5c2a13678..98c3d3677 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -29,8 +29,33 @@ namespace openmc { namespace simulation { +int current_batch; +int current_gen; +int64_t current_work; +double keff {1.0}; +double keff_std; +double k_col_abs {0.0}; +double k_col_tra {0.0}; +double k_abs_tra {0.0}; +double log_spacing; +int n_lost_particles {0}; +bool need_depletion_rx {false}; +int restart_batch; +bool satisfy_triggers {false}; +bool simulation_initialized {false}; +int total_gen {0}; +int64_t work; + +std::vector k_generation; std::vector work_index; +// Threadprivate variables +bool trace; //!< flag to show debug information +#ifdef _OPENMP +int n_threads {-1}; //!< number of OpenMP threads +int thread_id; //!< ID of a given thread +#endif + } // namespace simulation //============================================================================== @@ -43,6 +68,12 @@ void openmc_simulation_init_c() calculate_work(); } +int overall_generation() +{ + using namespace simulation; + return settings::gen_per_batch*(current_batch - 1) + current_gen; +} + void calculate_work() { // Determine minimum amount of particles to simulate on each processor diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 8328b13b1..29ee34ccf 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -13,7 +13,7 @@ module simulation_header ! GEOMETRY-RELATED VARIABLES ! Number of lost particles - integer(C_INT), bind(C) :: n_lost_particles = 0 + integer(C_INT), bind(C) :: n_lost_particles real(C_DOUBLE), bind(C) :: log_spacing ! spacing on logarithmic grid @@ -22,14 +22,14 @@ module simulation_header integer(C_INT), bind(C) :: current_batch ! current batch integer(C_INT), bind(C) :: current_gen ! current generation within a batch - integer(C_INT), bind(C) :: total_gen = 0 ! total number of generations simulated - logical(C_BOOL), bind(C) :: simulation_initialized = .false. + integer(C_INT), bind(C) :: total_gen ! total number of generations simulated + logical(C_BOOL), bind(C) :: simulation_initialized logical(C_BOOL), bind(C) :: need_depletion_rx ! need to calculate depletion reaction rx? ! ============================================================================ ! TALLY PRECISION TRIGGER VARIABLES - logical(C_BOOL), bind(C) :: satisfy_triggers = .false. ! whether triggers are satisfied + logical(C_BOOL), bind(C) :: satisfy_triggers ! whether triggers are satisfied integer(C_INT64_T), bind(C) :: work ! number of particles per processor integer(C_INT64_T), allocatable :: work_index(:) ! starting index in source bank for each process @@ -40,18 +40,18 @@ module simulation_header ! Temporary k-effective values type(VectorReal) :: k_generation ! single-generation estimates of k - real(C_DOUBLE), bind(C) :: keff = ONE ! average k over active batches - real(C_DOUBLE), bind(C) :: keff_std ! standard deviation of average k - real(C_DOUBLE), bind(C) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption - real(C_DOUBLE), bind(C) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength - real(C_DOUBLE), bind(C) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength + real(C_DOUBLE), bind(C) :: keff ! average k over active batches + real(C_DOUBLE), bind(C) :: keff_std ! standard deviation of average k + real(C_DOUBLE), bind(C) :: k_col_abs ! sum over batches of k_collision * k_absorption + real(C_DOUBLE), bind(C) :: k_col_tra ! sum over batches of k_collision * k_tracklength + real(C_DOUBLE), bind(C) :: k_abs_tra ! sum over batches of k_absorption * k_tracklength ! ============================================================================ ! PARALLEL PROCESSING VARIABLES #ifdef _OPENMP - integer(C_INT), bind(C) :: n_threads = NONE ! number of OpenMP threads - integer(C_INT), bind(C) :: thread_id ! ID of a given thread + integer(C_INT), bind(C) :: n_threads ! number of OpenMP threads + integer(C_INT), bind(C) :: thread_id ! ID of a given thread #endif ! ============================================================================ @@ -66,20 +66,15 @@ module simulation_header interface subroutine entropy_clear() bind(C) end subroutine + + pure function overall_generation() result(gen) bind(C) + import C_INT + integer(C_INT) :: gen + end function overall_generation end interface - contains -!=============================================================================== -! OVERALL_GENERATION determines the overall generation number -!=============================================================================== - - pure function overall_generation() result(gen) bind(C) - integer(C_INT) :: gen - gen = gen_per_batch*(current_batch - 1) + current_gen - end function overall_generation - !=============================================================================== ! FREE_MEMORY_SIMULATION deallocates global arrays defined in this module !=============================================================================== From f2748f64a91894aa70ef07471928b67832af14f6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 26 Apr 2018 22:26:29 -0500 Subject: [PATCH 06/23] Perform end of simulation broadcasts from C++ --- include/openmc/capi.h | 1 + include/openmc/simulation.h | 9 +++++-- src/simulation.F90 | 43 ++++++++------------------------- src/simulation.cpp | 48 +++++++++++++++++++++++++++++++++++++ 4 files changed, 66 insertions(+), 35 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 8e9edc4a0..5aba98453 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -48,6 +48,7 @@ extern "C" { int64_t openmc_get_seed(); int openmc_get_tally_index(int32_t id, int32_t* index); void openmc_get_tally_next_id(int32_t* id); + int openmc_global_tallies(double** ptr); int openmc_hard_reset(); int openmc_init(int argc, char* argv[], const void* intracomm); int openmc_init_f(const int* intracomm); diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 564516e60..ceb5eb92a 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -50,14 +50,19 @@ extern "C" int thread_id; //!< ID of a given thread // Functions //============================================================================== +//! Determine number of particles to transport per process +void calculate_work(); + //! Initialize simulation extern "C" void openmc_simulation_init_c(); //! Determine overall generation number extern "C" int overall_generation(); -//! Determine number of particles to transport per process -void calculate_work(); +#ifdef OPENMC_MPI +extern "C" void broadcast_results(); +extern "C" void broadcast_triggers(); +#endif } // namespace openmc diff --git a/src/simulation.F90 b/src/simulation.F90 index b1087e21d..c4158062e 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -330,6 +330,11 @@ contains #endif character(MAX_FILE_LEN) :: filename + interface + subroutine broadcast_triggers() bind(C) + end subroutine broadcast_triggers + end interface + ! Reduce tallies onto master process and accumulate call time_tallies % start() call accumulate_tallies() @@ -351,8 +356,7 @@ contains ! Check_triggers if (master) call check_triggers() #ifdef OPENMC_MPI - call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, & - mpi_intracomm, mpi_err) + call broadcast_triggers() #endif if (satisfy_triggers .or. & (trigger_on .and. current_batch == n_max_batches)) then @@ -488,6 +492,9 @@ contains interface subroutine print_overlap_check() bind(C) end subroutine print_overlap_check + + subroutine broadcast_results() bind(C) + end subroutine broadcast_results end interface err = 0 @@ -515,37 +522,7 @@ contains total_gen = total_gen + current_batch*gen_per_batch #ifdef OPENMC_MPI - ! Broadcast tally results so that each process has access to results - if (allocated(tallies)) then - do i = 1, size(tallies) - associate (results => tallies(i) % obj % results) - ! Create a new datatype that consists of all values for a given filter - ! bin and then use that to broadcast. This is done to minimize the - ! chance of the 'count' argument of MPI_BCAST exceeding 2**31 - n = size(results, 3) - count_per_filter = size(results, 1) * size(results, 2) - call MPI_TYPE_CONTIGUOUS(count_per_filter, MPI_DOUBLE, & - result_block, mpi_err) - call MPI_TYPE_COMMIT(result_block, mpi_err) - call MPI_BCAST(results, n, result_block, 0, mpi_intracomm, mpi_err) - call MPI_TYPE_FREE(result_block, mpi_err) - end associate - end do - end if - - ! Also broadcast global tally results - n = size(global_tallies) - call MPI_BCAST(global_tallies, n, MPI_DOUBLE, 0, mpi_intracomm, mpi_err) - - ! These guys are needed so that non-master processes can calculate the - ! combined estimate of k-effective - tempr(1) = k_col_abs - tempr(2) = k_col_tra - tempr(3) = k_abs_tra - call MPI_BCAST(tempr, 3, MPI_REAL8, 0, mpi_intracomm, mpi_err) - k_col_abs = tempr(1) - k_col_tra = tempr(2) - k_abs_tra = tempr(3) + call broadcast_results() #endif ! Write tally results to tallies.out diff --git a/src/simulation.cpp b/src/simulation.cpp index 98c3d3677..249dd090a 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -4,6 +4,12 @@ #include "openmc/message_passing.h" #include "openmc/settings.h" +#ifdef OPENMC_MPI +#include +#endif + +#include + // OPENMC_RUN encompasses all the main logic where iterations are performed // over the batches, generations, and histories in a fixed source or k-eigenvalue // calculation. @@ -98,4 +104,46 @@ void calculate_work() } } +#ifdef OPENMC_MPI +void broadcast_results() { + // Broadcast tally results so that each process has access to results + double* buffer; + if (n_tallies > 0) { + for (int i=1; i <= n_tallies; ++i) { + // Create a new datatype that consists of all values for a given filter + // bin and then use that to broadcast. This is done to minimize the + // chance of the 'count' argument of MPI_BCAST exceeding 2**31 + int shape[3]; + openmc_tally_results(i, &buffer, shape); + + int n = shape[0]; + int count_per_filter = shape[1] * shape[2]; + MPI_Datatype result_block; + MPI_Type_contiguous(count_per_filter, MPI_DOUBLE, &result_block); + MPI_Type_commit(&result_block); + MPI_Bcast(buffer, n, result_block, 0, mpi::intracomm); + MPI_Type_free(&result_block); + } + } + + // Also broadcast global tally results + openmc_global_tallies(&buffer); + MPI_Bcast(buffer, 4, MPI_DOUBLE, 0, mpi::intracomm); + + // These guys are needed so that non-master processes can calculate the + // combined estimate of k-effective + double temp[] {simulation::k_col_abs, simulation::k_col_tra, + simulation::k_abs_tra}; + MPI_Bcast(temp, 3, MPI_DOUBLE, 0, mpi::intracomm); + simulation::k_col_abs = temp[0]; + simulation::k_col_tra = temp[1]; + simulation::k_abs_tra = temp[2]; +} + +void broadcast_triggers() +{ + MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm); +} +#endif + } // namespace openmc From 6ba4d68cb9cfa3132c32f08848781bffc846c8d5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 10 Oct 2018 07:50:30 -0500 Subject: [PATCH 07/23] Move call to MPI_Abort to C++ side --- CMakeLists.txt | 1 + include/openmc/error.h | 5 ++++- include/openmc/message_passing.h | 2 +- src/error.F90 | 12 +++++++++--- src/error.cpp | 14 ++++++++++++++ 5 files changed, 29 insertions(+), 5 deletions(-) create mode 100644 src/error.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index b72bb57eb..b7517a130 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -407,6 +407,7 @@ add_library(libopenmc SHARED src/distribution_spatial.cpp src/eigenvalue.cpp src/endf.cpp + src/error.cpp src/initialize.cpp src/finalize.cpp src/geometry.cpp diff --git a/include/openmc/error.h b/include/openmc/error.h index 17b2f392e..b02c0a9f8 100644 --- a/include/openmc/error.h +++ b/include/openmc/error.h @@ -9,7 +9,6 @@ namespace openmc { - extern "C" void fatal_error_from_c(const char* message, int message_len); extern "C" void warning_from_c(const char* message, int message_len); extern "C" void write_message_from_c(const char* message, int message_len, @@ -81,5 +80,9 @@ void write_message(const std::stringstream& message, int level) write_message(message.str(), level); } +#ifdef OPENMC_MPI +extern "C" void abort_mpi(int code); +#endif + } // namespace openmc #endif // OPENMC_ERROR_H diff --git a/include/openmc/message_passing.h b/include/openmc/message_passing.h index c910210db..c730674ab 100644 --- a/include/openmc/message_passing.h +++ b/include/openmc/message_passing.h @@ -2,7 +2,7 @@ #define OPENMC_MESSAGE_PASSING_H #ifdef OPENMC_MPI -#include "mpi.h" +#include #endif namespace openmc { diff --git a/src/error.F90 b/src/error.F90 index c00006c8d..ad07987eb 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -130,14 +130,20 @@ contains character(*) :: message integer, optional :: error_code ! error code - integer :: code ! error code + integer(C_INT) :: code ! error code integer :: i_start ! starting position integer :: i_end ! ending position integer :: line_wrap ! length of line integer :: length ! length of message integer :: indent ! length of indentation + #ifdef OPENMC_MPI - integer :: mpi_err + interface + subroutine abort_mpi(code) bind(C) + import C_INT + integer(C_INT), value :: code + end subroutine + end interface #endif @@ -190,7 +196,7 @@ contains #ifdef OPENMC_MPI ! Abort MPI - call MPI_ABORT(mpi_intracomm, code, mpi_err) + call abort_mpi(code) #endif ! Abort program diff --git a/src/error.cpp b/src/error.cpp new file mode 100644 index 000000000..d17219f31 --- /dev/null +++ b/src/error.cpp @@ -0,0 +1,14 @@ +#include "openmc/error.h" + +#include "openmc/message_passing.h" + +namespace openmc { + +#ifdef OPENMC_MPI +void abort_mpi(int code) +{ + MPI_Abort(mpi::intracomm, code); +} +#endif + +} // namespace openmc From ab4167452719095e726b557ec5048bfce9825629 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 10 Oct 2018 07:56:34 -0500 Subject: [PATCH 08/23] Remove unused variables --- src/simulation.F90 | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/src/simulation.F90 b/src/simulation.F90 index c4158062e..767ec5244 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -325,9 +325,6 @@ contains subroutine finalize_batch() integer(C_INT) :: err -#ifdef OPENMC_MPI - integer :: mpi_err ! MPI error code -#endif character(MAX_FILE_LEN) :: filename interface @@ -477,17 +474,6 @@ contains integer(C_INT) :: err integer :: i ! loop index -#ifdef OPENMC_MPI - integer :: n ! size of arrays - integer :: mpi_err ! MPI error code - integer :: count_per_filter ! number of result values for one filter bin - real(8) :: tempr(3) ! temporary array for communication -#ifdef OPENMC_MPIF08 - type(MPI_Datatype) :: result_block -#else - integer :: result_block -#endif -#endif interface subroutine print_overlap_check() bind(C) From 6b6dbb73c571847057763f429e2d7629b98b3bd9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 10 Oct 2018 10:00:42 -0500 Subject: [PATCH 09/23] Reverse shape for openmc_tally_results (makes sense for C/C++) --- openmc/capi/tally.py | 2 +- src/tallies/tally_header.F90 | 9 +++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index 44bf89415..56c4799aa 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -298,7 +298,7 @@ class Tally(_FortranObjectWithID): data = POINTER(c_double)() shape = (c_int*3)() _dll.openmc_tally_results(self._index, data, shape) - return as_array(data, tuple(shape[::-1])) + return as_array(data, tuple(shape)) @property def scores(self): diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index c6238a88d..85cc5133b 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -153,7 +153,7 @@ module tally_header ! Normalization for statistics integer(C_INT32_T), public, bind(C) :: n_realizations = 0 ! # of independent realizations - real(8), public :: total_weight ! total starting particle weight in realization + real(C_DOUBLE), public, bind(C) :: total_weight ! total starting particle weight in realization contains @@ -688,7 +688,12 @@ contains associate (t => tallies(index) % obj) if (allocated(t % results)) then ptr = C_LOC(t % results(1,1,1)) - shape_(:) = shape(t % results) + + ! Note that shape is reversed since it is assumed to be used from + ! C/C++ code + shape_(1) = size(t % results, 3) + shape_(2) = size(t % results, 2) + shape_(3) = size(t % results, 1) err = 0 else err = E_ALLOCATE From 1fd4135f8fda64ddaf751bcc399d13602d014e67 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 10 Oct 2018 11:43:53 -0500 Subject: [PATCH 10/23] Move reduce_tally_results to C++ --- CMakeLists.txt | 3 +- include/openmc/constants.h | 7 +-- include/openmc/tallies/tally.h | 23 ++++++++ src/tallies/tally.F90 | 66 ++--------------------- src/tallies/tally.cpp | 96 ++++++++++++++++++++++++++++++++++ 5 files changed, 130 insertions(+), 65 deletions(-) create mode 100644 include/openmc/tallies/tally.h create mode 100644 src/tallies/tally.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index b7517a130..99d4e4a5e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -445,7 +445,8 @@ add_library(libopenmc SHARED src/surface.cpp src/thermal.cpp src/xml_interface.cpp - src/xsdata.cpp) + src/xsdata.cpp + src/tallies/tally.cpp) set_target_properties(libopenmc PROPERTIES OUTPUT_NAME openmc diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 9c96258fb..6ead514d7 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -314,9 +314,9 @@ constexpr int MG_GET_XS_CHI_DELAYED {14}; // TALLY-RELATED CONSTANTS // Tally result entries -constexpr int RESULT_VALUE {1}; -constexpr int RESULT_SUM {2}; -constexpr int RESULT_SUM_SQ {3}; +constexpr int RESULT_VALUE {0}; +constexpr int RESULT_SUM {1}; +constexpr int RESULT_SUM_SQ {2}; // Tally type // TODO: Convert to enum @@ -407,6 +407,7 @@ constexpr int RELATIVE_ERROR {2}; constexpr int STANDARD_DEVIATION {3}; // Global tally parameters +constexpr int N_GLOBAL_TALLIES {4}; constexpr int K_COLLISION {1}; constexpr int K_ABSORPTION {2}; constexpr int K_TRACKLENGTH {3}; diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h new file mode 100644 index 000000000..a68278bb1 --- /dev/null +++ b/include/openmc/tallies/tally.h @@ -0,0 +1,23 @@ +#ifndef OPENMC_TALLIES_TALLY_H +#define OPENMC_TALLIES_TALLY_H + +#include "openmc/constants.h" + +#include "xtensor/xtensor.hpp" + +namespace openmc { + +extern "C" double total_weight; + +xt::xtensor tally_results(int idx); + +xt::xtensor global_tallies(); + +#ifdef OPENMC_MPI +//! Collect all tally results onto master process +extern "C" void reduce_tally_results(); +#endif + +} // namespace openmc + +#endif // OPENMC_TALLIES_TALLY_H diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index a618f351e..ddba2dc07 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -3778,6 +3778,11 @@ contains real(C_DOUBLE) :: val #ifdef OPENMC_MPI + interface + subroutine reduce_tally_results() bind(C) + end subroutine + end interface + ! Combine tally results onto master process if (reduce_tallies) call reduce_tally_results() #endif @@ -3821,67 +3826,6 @@ contains end subroutine accumulate_tallies -!=============================================================================== -! REDUCE_TALLY_RESULTS collects all the results from tallies onto one processor -!=============================================================================== - -#ifdef OPENMC_MPI - subroutine reduce_tally_results() - - integer :: i - integer :: n ! number of filter bins - integer :: m ! number of score bins - integer :: n_bins ! total number of bins - integer :: mpi_err ! MPI error code - real(C_DOUBLE), allocatable :: tally_temp(:,:) ! contiguous array of results - real(C_DOUBLE), allocatable :: tally_temp2(:,:) ! reduced contiguous results - real(C_DOUBLE) :: temp(N_GLOBAL_TALLIES), temp2(N_GLOBAL_TALLIES) - - do i = 1, active_tallies % size() - associate (t => tallies(active_tallies % data(i)) % obj) - - m = size(t % results, 2) - n = size(t % results, 3) - n_bins = m*n - - allocate(tally_temp(m,n), tally_temp2(m,n)) - - ! Reduce contiguous set of tally results - tally_temp = t % results(RESULT_VALUE,:,:) - call MPI_REDUCE(tally_temp, tally_temp2, n_bins, MPI_DOUBLE, & - MPI_SUM, 0, mpi_intracomm, mpi_err) - - if (master) then - ! Transfer values to value on master - t % results(RESULT_VALUE,:,:) = tally_temp2 - else - ! Reset value on other processors - t % results(RESULT_VALUE,:,:) = ZERO - end if - - deallocate(tally_temp, tally_temp2) - end associate - end do - - ! Reduce global tallies onto master - temp = global_tallies(RESULT_VALUE, :) - call MPI_REDUCE(temp, temp2, N_GLOBAL_TALLIES, MPI_DOUBLE, MPI_SUM, & - 0, mpi_intracomm, mpi_err) - if (master) then - global_tallies(RESULT_VALUE, :) = temp2 - else - global_tallies(RESULT_VALUE, :) = ZERO - end if - - ! We also need to determine the total starting weight of particles from the - ! last realization - temp(1) = total_weight - call MPI_REDUCE(temp, total_weight, 1, MPI_REAL8, MPI_SUM, & - 0, mpi_intracomm, mpi_err) - - end subroutine reduce_tally_results -#endif - !=============================================================================== ! SETUP_ACTIVE_TALLIES !=============================================================================== diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp new file mode 100644 index 000000000..c98442b2c --- /dev/null +++ b/src/tallies/tally.cpp @@ -0,0 +1,96 @@ +#include "openmc/tallies/tally.h" + +#include "openmc/capi.h" +#include "openmc/constants.h" +#include "openmc/message_passing.h" + +#include "xtensor/xadapt.hpp" +#include "xtensor/xbuilder.hpp" // for empty_like +#include "xtensor/xview.hpp" + +#include + +namespace openmc { + +xt::xtensor global_tallies() +{ + // Get pointer to global tallies + double* buffer; + openmc_global_tallies(&buffer); + + // Adapt into xtensor + std::vector shape {N_GLOBAL_TALLIES, 3}; + int size {3*N_GLOBAL_TALLIES}; + return xt::adapt(buffer, size, xt::no_ownership(), shape); +} + +xt::xtensor tally_results(int idx) +{ + // Get pointer to tally results + double* results; + std::vector shape(3); + openmc_tally_results(idx, &results, shape.data()); + + // Adapt array into xtensor with no ownership + int size {shape[0] * shape[1] * shape[2]}; + return xt::adapt(results, size, xt::no_ownership(), shape); +} + +#ifdef OPENMC_MPI +void reduce_tally_results() +{ + for (int i = 1; i <= n_tallies; ++i) { + // Skip any tallies that are not active + bool active; + openmc_tally_get_active(i, &active); + if (!active) continue; + + // Get view of accumulated tally values + auto results = tally_results(i); + auto values_view = xt::view(results, xt::all(), xt::all(), RESULT_VALUE); + + // Make copy of tally values in contiguous array + xt::xtensor values = values_view; + xt::xtensor values_reduced = xt::empty_like(values); + + // Reduce contiguous set of tally results + MPI_Reduce(values.data(), values_reduced.data(), values.size(), + MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); + + // Transfer values on master and reset on other ranks + if (mpi::master) { + values_view = mpi::master; + } else { + values_view = 0.0; + } + } + + // Get view of global tally values + auto gt = global_tallies(); + auto gt_values_view = xt::view(gt, xt::all(), RESULT_VALUE); + + // Make copy of values in contiguous array + xt::xtensor gt_values = gt_values_view; + xt::xtensor gt_values_reduced = xt::empty_like(gt_values); + + // Reduce contiguous data + MPI_Reduce(gt_values.data(), gt_values_reduced.data(), N_GLOBAL_TALLIES, + MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); + + // Transfer values on master and reset on other ranks + if (mpi::master) { + gt_values_view = gt_values_reduced; + } else { + gt_values_view = 0.0; + } + + // We also need to determine the total starting weight of particles from the + // last realization + double weight_reduced; + MPI_Reduce(&total_weight, &weight_reduced, 1, MPI_DOUBLE, MPI_SUM, + 0, mpi::intracomm); + if (mpi::master) total_weight = weight_reduced; +} +#endif + +} // namespace openmc From 2ccb7e3b5877d17c953dda895ffef7abb8b6c0cd Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 10 Oct 2018 13:23:16 -0500 Subject: [PATCH 11/23] Move write_tally_results_nr to C++ --- include/openmc/hdf5_interface.h | 9 +++ include/openmc/state_point.h | 2 + src/state_point.F90 | 127 +------------------------------- src/state_point.cpp | 123 +++++++++++++++++++++++++++++-- 4 files changed, 131 insertions(+), 130 deletions(-) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index fa6460551..cf7fb732f 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -368,6 +368,15 @@ write_dataset(hid_t obj_id, const char* name, const xt::xarray& arr) arr.data(), false); } +template inline void +write_dataset(hid_t obj_id, const char* name, const xt::xtensor& arr) +{ + auto s = arr.shape(); + std::vector dims {s.cbegin(), s.cend()}; + write_dataset(obj_id, dims.size(), dims.data(), name, H5TypeMap::type_id, + arr.data(), false); +} + inline void write_dataset(hid_t obj_id, const char* name, Position r) { diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index 64713ae96..6541a67fc 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -14,5 +14,7 @@ extern "C" void write_source_bank(hid_t group_id, int64_t* work_index, extern "C" void read_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank); +extern "C" void write_tally_results_nr(hid_t file_id); + } // namespace openmc #endif // OPENMC_STATE_POINT_H diff --git a/src/state_point.F90 b/src/state_point.F90 index 91d0f1090..11e4ad621 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -87,6 +87,10 @@ contains import HID_T integer(HID_T), value :: group end subroutine + subroutine write_tally_results_nr(file_id) bind(C) + import HID_T + integer(HID_T), value :: file_id + end subroutine end interface err = 0 @@ -490,129 +494,6 @@ contains end subroutine write_source_point -!=============================================================================== -! WRITE_TALLY_RESULTS_NR -!=============================================================================== - - subroutine write_tally_results_nr(file_id) - integer(HID_T), intent(in) :: file_id - - integer :: i ! loop index - integer :: n ! number of filter bins - integer :: m ! number of score bins - integer :: n_bins ! total number of bins - integer(HID_T) :: tallies_group, tally_group - real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of results - real(8), target :: global_temp(3,N_GLOBAL_TALLIES) -#ifdef OPENMC_MPI - integer :: mpi_err ! MPI error code - real(8) :: dummy ! temporary receive buffer for non-root reduces -#endif - type(TallyObject) :: dummy_tally - - ! ========================================================================== - ! COLLECT AND WRITE GLOBAL TALLIES - - if (master) then - ! Write number of realizations - call write_dataset(file_id, "n_realizations", n_realizations) - - ! Write number of global tallies - call write_dataset(file_id, "n_global_tallies", N_GLOBAL_TALLIES) - - tallies_group = open_group(file_id, "tallies") - end if - - -#ifdef OPENMC_MPI - ! Reduce global tallies - n_bins = size(global_tallies) - call MPI_REDUCE(global_tallies, global_temp, n_bins, MPI_REAL8, MPI_SUM, & - 0, mpi_intracomm, mpi_err) -#endif - - if (master) then - ! Transfer values to value on master - if (current_batch == n_max_batches .or. satisfy_triggers) then - global_tallies(:,:) = global_temp(:,:) - end if - - ! Write out global tallies sum and sum_sq - call write_dataset(file_id, "global_tallies", global_temp) - end if - - if (active_tallies % size() > 0) then - ! Indicate that tallies are on - if (master) then - call write_attribute(file_id, "tallies_present", 1) - end if - - ! Write all tally results - TALLY_RESULTS: do i = 1, n_tallies - associate (t => tallies(i) % obj) - ! Determine size of tally results array - m = size(t % results, 2) - n = size(t % results, 3) - n_bins = m*n*2 - - ! Allocate array for storing sums and sums of squares, but - ! contiguously in memory for each - allocate(tally_temp(2,m,n)) - tally_temp(1,:,:) = t % results(RESULT_SUM,:,:) - tally_temp(2,:,:) = t % results(RESULT_SUM_SQ,:,:) - - if (master) then - tally_group = open_group(tallies_group, "tally " // & - trim(to_str(t % id))) - - ! The MPI_IN_PLACE specifier allows the master to copy values into - ! a receive buffer without having a temporary variable -#ifdef OPENMC_MPI - call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, & - MPI_SUM, 0, mpi_intracomm, mpi_err) -#endif - - ! At the end of the simulation, store the results back in the - ! regular TallyResults array - if (current_batch == n_max_batches .or. satisfy_triggers) then - t % results(RESULT_SUM,:,:) = tally_temp(1,:,:) - t % results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:) - end if - - ! Put in temporary tally result - allocate(dummy_tally % results(3,m,n)) - dummy_tally % results(RESULT_SUM,:,:) = tally_temp(1,:,:) - dummy_tally % results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:) - - ! Write reduced tally results to file - call dummy_tally % write_results_hdf5(tally_group) - - ! Deallocate temporary tally result - deallocate(dummy_tally % results) - else - ! Receive buffer not significant at other processors -#ifdef OPENMC_MPI - call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, & - 0, mpi_intracomm, mpi_err) -#endif - end if - - ! Deallocate temporary copy of tally results - deallocate(tally_temp) - - if (master) call close_group(tally_group) - end associate - end do TALLY_RESULTS - - if (master) call close_group(tallies_group) - else - ! Indicate that tallies are off - if (master) call write_dataset(file_id, "tallies_present", 0) - end if - - - end subroutine write_tally_results_nr - !=============================================================================== ! LOAD_STATE_POINT !=============================================================================== diff --git a/src/state_point.cpp b/src/state_point.cpp index e24993b6c..4bd91e121 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -1,17 +1,20 @@ #include "openmc/state_point.h" -#include -#include - -#ifdef OPENMC_MPI -#include "mpi.h" -#endif - #include "openmc/capi.h" +#include "openmc/constants.h" #include "openmc/error.h" +#include "openmc/hdf5_interface.h" #include "openmc/message_passing.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include "openmc/tallies/tally.h" + +#include "xtensor/xbuilder.hpp" // for empty_like +#include "xtensor/xview.hpp" + +#include +#include +#include namespace openmc { @@ -165,4 +168,110 @@ void read_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) H5Tclose(banktype); } +void write_tally_results_nr(hid_t file_id) +{ + // ========================================================================== + // COLLECT AND WRITE GLOBAL TALLIES + + hid_t tallies_group; + if (mpi::master) { + // Write number of realizations + write_dataset(file_id, "n_realizations", n_realizations); + + // Write number of global tallies + write_dataset(file_id, "n_global_tallies", N_GLOBAL_TALLIES); + + tallies_group = open_group(file_id, "tallies"); + } + + // Get pointer to global tallies + auto gt = global_tallies(); + +#ifdef OPENMC_MPI + // Reduce global tallies + xt::xtensor gt_reduced = xt::empty_like(gt); + MPI_Reduce(gt.data(), gt_reduced.data(), gt.size(), MPI_DOUBLE, + MPI_SUM, 0, mpi::intracomm); + + // Transfer values to value on master + if (mpi::master) { + if (simulation::current_batch == settings::n_max_batches || + simulation::satisfy_triggers) { + std::copy(gt_reduced.begin(), gt_reduced.end(), gt.begin()); + } + } +#endif + + // Write out global tallies sum and sum_sq + if (mpi::master) { + write_dataset(file_id, "global_tallies", gt); + } + + for (int i = 1; i <= n_tallies; ++i) { + // Skip any tallies that are not active + bool active; + openmc_tally_get_active(i, &active); + if (!active) continue; + + if (mpi::master && !object_exists(file_id, "tallies_present")) { + write_attribute(file_id, "tallies_present", 1); + } + + // Get view of accumulated tally values + auto results = tally_results(i); + auto values_view = xt::view(results, xt::all(), xt::all(), + xt::range(RESULT_SUM, RESULT_SUM_SQ + 1)); + + // Make copy of tally values in contiguous array + xt::xtensor values = values_view; + + if (mpi::master) { + // Open group for tally + int id; + openmc_tally_get_id(i, &id); + std::string groupname {"tally " + std::to_string(id)}; + hid_t tally_group = open_group(tallies_group, groupname.c_str()); + + // The MPI_IN_PLACE specifier allows the master to copy values into + // a receive buffer without having a temporary variable +#ifdef OPENMC_MPI + MPI_Reduce(MPI_IN_PLACE, values.data(), values.size(), MPI_DOUBLE, + MPI_SUM, 0, mpi::intracomm); +#endif + + // At the end of the simulation, store the results back in the + // regular TallyResults array + if (simulation::current_batch == settings::n_max_batches || + simulation::satisfy_triggers) { + values_view = values; + } + + // Put in temporary tally result + xt::xtensor results_copy = xt::zeros_like(results); + auto copy_view = xt::view(results_copy, xt::all(), xt::all(), + xt::range(RESULT_SUM, RESULT_SUM_SQ + 1)); + copy_view = values; + + // Write reduced tally results to file + auto shape = results_copy.shape(); + write_tally_results(tally_group, shape[0], shape[1], results_copy.data()); + + close_group(tally_group); + } else { + // Receive buffer not significant at other processors +#ifdef OPENMC_MPI + MPI_Reduce(values.data(), nullptr, values.size(), MPI_REAL8, MPI_SUM, + 0, mpi::intracomm); +#endif + } + } + + if (mpi::master) { + if (!object_exists(file_id, "tallies_present")) { + // Indicate that tallies are off + write_dataset(file_id, "tallies_present", 0); + } + } +} + } // namespace openmc From 9c787f988e2a87995b43a02a49c3e2573b1d8243 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 11 Oct 2018 09:35:51 -0500 Subject: [PATCH 12/23] Fix tally reduction --- include/openmc/capi.h | 3 ++- include/openmc/hdf5_interface.h | 20 +++++++---------- include/openmc/tallies/tally.h | 6 +++-- openmc/capi/tally.py | 6 ++--- src/simulation.cpp | 39 ++++++++++++++------------------- src/tallies/tally.cpp | 16 ++++++++------ src/tallies/tally_header.F90 | 2 +- 7 files changed, 43 insertions(+), 49 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 5aba98453..43f5f9869 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -3,6 +3,7 @@ #include #include +#include #ifdef __cplusplus extern "C" { @@ -106,7 +107,7 @@ extern "C" { int openmc_tally_get_scores(int32_t index, int** scores, int* n); int openmc_tally_get_type(int32_t index, int32_t* type); int openmc_tally_reset(int32_t index); - int openmc_tally_results(int32_t index, double** ptr, int shape_[3]); + int openmc_tally_results(int32_t index, double** ptr, size_t shape_[3]); int openmc_tally_set_active(int32_t index, bool active); int openmc_tally_set_estimator(int32_t index, const char* estimator); int openmc_tally_set_filters(int32_t index, int n, const int32_t* indices); diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index cf7fb732f..cce7ca164 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -8,6 +8,7 @@ #include // for strlen #include #include +#include #include #include "hdf5.h" @@ -333,7 +334,9 @@ write_attribute(hid_t obj_id, const char* name, const std::vector& buffer) // Templates/overloads for write_dataset //============================================================================== -template inline void +// Template for scalars (ensured by SFINAE) +template inline +std::enable_if_t>::value> write_dataset(hid_t obj_id, const char* name, T buffer) { write_dataset(obj_id, 0, nullptr, name, H5TypeMap::type_id, &buffer, false); @@ -359,18 +362,11 @@ write_dataset(hid_t obj_id, const char* name, const std::vector& buffer) write_dataset(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data(), false); } -template inline void -write_dataset(hid_t obj_id, const char* name, const xt::xarray& arr) -{ - auto s = arr.shape(); - std::vector dims {s.cbegin(), s.cend()}; - write_dataset(obj_id, dims.size(), dims.data(), name, H5TypeMap::type_id, - arr.data(), false); -} - -template inline void -write_dataset(hid_t obj_id, const char* name, const xt::xtensor& arr) +// Template for xarray, xtensor, etc. +template inline void +write_dataset(hid_t obj_id, const char* name, const xt::xcontainer& arr) { + using T = typename D::value_type; auto s = arr.shape(); std::vector dims {s.cbegin(), s.cend()}; write_dataset(obj_id, dims.size(), dims.data(), name, H5TypeMap::type_id, diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index a68278bb1..a07cd6bd1 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -9,9 +9,11 @@ namespace openmc { extern "C" double total_weight; -xt::xtensor tally_results(int idx); +template +using adaptor_type = xt::xtensor_adaptor, N>; -xt::xtensor global_tallies(); +adaptor_type<2> global_tallies(); +adaptor_type<3> tally_results(int idx); #ifdef OPENMC_MPI //! Collect all tally results onto master process diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index 56c4799aa..a414ead45 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -1,5 +1,5 @@ from collections.abc import Mapping -from ctypes import c_int, c_int32, c_double, c_char_p, c_bool, POINTER +from ctypes import c_int, c_int32, c_size_t, c_double, c_char_p, c_bool, POINTER from weakref import WeakValueDictionary import numpy as np @@ -60,7 +60,7 @@ _dll.openmc_tally_reset.argtypes = [c_int32] _dll.openmc_tally_reset.restype = c_int _dll.openmc_tally_reset.errcheck = _error_handler _dll.openmc_tally_results.argtypes = [ - c_int32, POINTER(POINTER(c_double)), POINTER(c_int*3)] + c_int32, POINTER(POINTER(c_double)), POINTER(c_size_t*3)] _dll.openmc_tally_results.restype = c_int _dll.openmc_tally_results.errcheck = _error_handler _dll.openmc_tally_set_active.argtypes = [c_int32, c_bool] @@ -296,7 +296,7 @@ class Tally(_FortranObjectWithID): @property def results(self): data = POINTER(c_double)() - shape = (c_int*3)() + shape = (c_size_t*3)() _dll.openmc_tally_results(self._index, data, shape) return as_array(data, tuple(shape)) diff --git a/src/simulation.cpp b/src/simulation.cpp index 249dd090a..386e00510 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -3,10 +3,7 @@ #include "openmc/capi.h" #include "openmc/message_passing.h" #include "openmc/settings.h" - -#ifdef OPENMC_MPI -#include -#endif +#include "openmc/tallies/tally.h" #include @@ -106,29 +103,25 @@ void calculate_work() #ifdef OPENMC_MPI void broadcast_results() { - // Broadcast tally results so that each process has access to results - double* buffer; - if (n_tallies > 0) { - for (int i=1; i <= n_tallies; ++i) { - // Create a new datatype that consists of all values for a given filter - // bin and then use that to broadcast. This is done to minimize the - // chance of the 'count' argument of MPI_BCAST exceeding 2**31 - int shape[3]; - openmc_tally_results(i, &buffer, shape); + // Broadcast tally results so that each process has access to results + for (int i = 1; i <= n_tallies; ++i) { + // Create a new datatype that consists of all values for a given filter + // bin and then use that to broadcast. This is done to minimize the + // chance of the 'count' argument of MPI_BCAST exceeding 2**31 + auto results = tally_results(i); - int n = shape[0]; - int count_per_filter = shape[1] * shape[2]; - MPI_Datatype result_block; - MPI_Type_contiguous(count_per_filter, MPI_DOUBLE, &result_block); - MPI_Type_commit(&result_block); - MPI_Bcast(buffer, n, result_block, 0, mpi::intracomm); - MPI_Type_free(&result_block); - } + auto shape = results.shape(); + int count_per_filter = shape[1] * shape[2]; + MPI_Datatype result_block; + MPI_Type_contiguous(count_per_filter, MPI_DOUBLE, &result_block); + MPI_Type_commit(&result_block); + MPI_Bcast(results.data(), shape[0], result_block, 0, mpi::intracomm); + MPI_Type_free(&result_block); } // Also broadcast global tally results - openmc_global_tallies(&buffer); - MPI_Bcast(buffer, 4, MPI_DOUBLE, 0, mpi::intracomm); + auto gt = global_tallies(); + MPI_Bcast(gt.data(), gt.size(), MPI_DOUBLE, 0, mpi::intracomm); // These guys are needed so that non-master processes can calculate the // combined estimate of k-effective diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index c98442b2c..98783f449 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -8,31 +8,33 @@ #include "xtensor/xbuilder.hpp" // for empty_like #include "xtensor/xview.hpp" +#include #include namespace openmc { -xt::xtensor global_tallies() +adaptor_type<2> global_tallies() { // Get pointer to global tallies double* buffer; openmc_global_tallies(&buffer); // Adapt into xtensor - std::vector shape {N_GLOBAL_TALLIES, 3}; - int size {3*N_GLOBAL_TALLIES}; + std::array shape = {N_GLOBAL_TALLIES, 3}; + std::size_t size {3*N_GLOBAL_TALLIES}; + return xt::adapt(buffer, size, xt::no_ownership(), shape); } -xt::xtensor tally_results(int idx) +adaptor_type<3> tally_results(int idx) { // Get pointer to tally results double* results; - std::vector shape(3); + std::array shape; openmc_tally_results(idx, &results, shape.data()); // Adapt array into xtensor with no ownership - int size {shape[0] * shape[1] * shape[2]}; + std::size_t size {shape[0] * shape[1] * shape[2]}; return xt::adapt(results, size, xt::no_ownership(), shape); } @@ -59,7 +61,7 @@ void reduce_tally_results() // Transfer values on master and reset on other ranks if (mpi::master) { - values_view = mpi::master; + values_view = values_reduced; } else { values_view = 0.0; } diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 85cc5133b..06eb6d2a2 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -681,7 +681,7 @@ contains ! allows a user to obtain in-memory tally results from Python directly. integer(C_INT32_T), intent(in), value :: index type(C_PTR), intent(out) :: ptr - integer(C_INT), intent(out) :: shape_(3) + integer(C_SIZE_T), intent(out) :: shape_(3) integer(C_INT) :: err if (index >= 1 .and. index <= size(tallies)) then From 217f87c0852abf58ee8708d0a7857b96380e1a0e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 15 Oct 2018 11:52:03 -0500 Subject: [PATCH 13/23] Move global_tally_absorption to tally.h --- include/openmc/capi.h | 2 -- include/openmc/tallies/tally.h | 19 +++++++++++++++++++ src/physics_mg.cpp | 1 + src/tallies/tally.cpp | 4 ++++ 4 files changed, 24 insertions(+), 2 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 43f5f9869..f7d69290c 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -148,8 +148,6 @@ extern "C" { extern int32_t n_surfaces; extern int32_t n_tallies; extern int32_t n_universes; - extern double global_tally_absorption; -#pragma omp threadprivate(global_tally_absorption) // Variables that are shared by necessity (can be removed from public header // later) diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index a07cd6bd1..96b081615 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -7,12 +7,31 @@ namespace openmc { +//============================================================================== +// Global variable declarations +//============================================================================== + extern "C" double total_weight; +// Threadprivate variables + +extern "C" double global_tally_absorption; +#pragma omp threadprivate(global_tally_absorption) + +//============================================================================== +// Non-member functions +//============================================================================== + template using adaptor_type = xt::xtensor_adaptor, N>; +//! Get the global tallies as a multidimensional array +//! \return Global tallies array adaptor_type<2> global_tallies(); + +//! Get tally results as a multidimensional array +//! \param idx Index in tallies array +//! \return Tally results array adaptor_type<3> tally_results(int idx); #ifdef OPENMC_MPI diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 18301c3ab..40e5ba1ae 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -16,6 +16,7 @@ #include "openmc/random_lcg.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include "openmc/tallies/tally.h" namespace openmc { diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 98783f449..481f05edc 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -13,6 +13,10 @@ namespace openmc { +//============================================================================== +// Non-member functions +//============================================================================== + adaptor_type<2> global_tallies() { // Get pointer to global tallies From 02fbc1582d837e40e01c7aa4c82283440cfa1599 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 11 Oct 2018 12:55:33 -0500 Subject: [PATCH 14/23] Convert synchronize_bank to C++ --- include/openmc/eigenvalue.h | 3 + include/openmc/search.h | 10 +- src/api.F90 | 1 - src/eigenvalue.F90 | 271 ------------------------------------ src/eigenvalue.cpp | 244 ++++++++++++++++++++++++++++++++ src/initialize.F90 | 56 +------- src/initialize.cpp | 22 +-- src/message_passing.F90 | 2 - src/simulation.F90 | 5 +- src/simulation.cpp | 6 +- 10 files changed, 278 insertions(+), 342 deletions(-) diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index cc879461e..59b448e29 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -24,6 +24,9 @@ extern "C" int64_t n_bank; // Non-member functions //============================================================================== +//! Sample/redistribute source sites from accumulated fission sites +extern "C" void synchronize_bank(); + //! Calculates the Shannon entropy of the fission source distribution to assess //! source convergence extern "C" void shannon_entropy(); diff --git a/include/openmc/search.h b/include/openmc/search.h index c91ad18fe..446e2916a 100644 --- a/include/openmc/search.h +++ b/include/openmc/search.h @@ -4,7 +4,7 @@ #ifndef OPENMC_SEARCH_H #define OPENMC_SEARCH_H -#include // for lower_bound +#include // for lower_bound, upper_bound namespace openmc { @@ -19,6 +19,14 @@ lower_bound_index(It first, It last, const T& value) return (index == last) ? -1 : index - first; } +template +typename std::iterator_traits::difference_type +upper_bound_index(It first, It last, const T& value) +{ + It index = std::upper_bound(first, last, value) - 1; + return (index == last) ? -1 : index - first; +} + } // namespace openmc #endif // OPENMC_SEARCH_H diff --git a/src/api.F90 b/src/api.F90 index 38329ad80..d79aa2450 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -181,7 +181,6 @@ contains err = 0 #ifdef OPENMC_MPI ! Free all MPI types - call MPI_TYPE_FREE(MPI_BANK, err) call openmc_free_bank() #endif diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index bd2a20b41..df108f6d7 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -21,277 +21,6 @@ module eigenvalue contains -!=============================================================================== -! SYNCHRONIZE_BANK samples source sites from the fission sites that were -! accumulated during the generation. This routine is what allows this Monte -! Carlo to scale to large numbers of processors where other codes cannot. -!=============================================================================== - - subroutine synchronize_bank() - - integer :: i ! loop indices - integer :: j ! loop indices - integer(8) :: start ! starting index in global bank - integer(8) :: finish ! ending index in global bank - integer(8) :: total ! total sites in global fission bank - integer(8) :: index_temp ! index in temporary source bank - integer(8) :: sites_needed ! # of sites to be sampled - real(8) :: p_sample ! probability of sampling a site - type(Bank), save, allocatable :: & - & temp_sites(:) ! local array of extra sites on each node - -#ifdef OPENMC_MPI - integer :: mpi_err ! MPI error code - integer(8) :: n ! number of sites to send/recv - integer :: neighbor ! processor to send/recv data from -#ifdef OPENMC_MPIF08 - type(MPI_Request) :: request(20) -#else - integer :: request(20) ! communication request for send/recving sites -#endif - integer :: n_request ! number of communication requests - integer(8) :: index_local ! index in local source bank - integer(8), save, allocatable :: & - & bank_position(:) ! starting positions in global source bank -#endif - - ! In order to properly understand the fission bank algorithm, you need to - ! think of the fission and source bank as being one global array divided - ! over multiple processors. At the start, each processor has a random amount - ! of fission bank sites -- each processor needs to know the total number of - ! sites in order to figure out the probability for selecting - ! sites. Furthermore, each proc also needs to know where in the 'global' - ! fission bank its own sites starts in order to ensure reproducibility by - ! skipping ahead to the proper seed. - -#ifdef OPENMC_MPI - start = 0_8 - call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, & - mpi_intracomm, mpi_err) - - ! While we would expect the value of start on rank 0 to be 0, the MPI - ! standard says that the receive buffer on rank 0 is undefined and not - ! significant - if (rank == 0) start = 0_8 - - finish = start + n_bank - total = finish - call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, & - mpi_intracomm, mpi_err) - -#else - start = 0_8 - finish = n_bank - total = n_bank -#endif - - ! If there are not that many particles per generation, it's possible that no - ! fission sites were created at all on a single processor. Rather than add - ! extra logic to treat this circumstance, we really want to ensure the user - ! runs enough particles to avoid this in the first place. - - if (n_bank == 0) then - call fatal_error("No fission sites banked on processor " // to_str(rank)) - end if - - ! Make sure all processors start at the same point for random sampling. Then - ! skip ahead in the sequence using the starting index in the 'global' - ! fission bank for each processor. - - call set_particle_seed(int(total_gen + overall_generation(), 8)) - call advance_prn_seed(start) - - ! Determine how many fission sites we need to sample from the source bank - ! and the probability for selecting a site. - - if (total < n_particles) then - sites_needed = mod(n_particles,total) - else - sites_needed = n_particles - end if - p_sample = real(sites_needed,8)/real(total,8) - - call time_bank_sample % start() - - ! ========================================================================== - ! SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES - - ! Allocate temporary source bank - index_temp = 0_8 - if (.not. allocated(temp_sites)) allocate(temp_sites(3*work)) - - do i = 1, int(n_bank,4) - - ! If there are less than n_particles particles banked, automatically add - ! int(n_particles/total) sites to temp_sites. For example, if you need - ! 1000 and 300 were banked, this would add 3 source sites per banked site - ! and the remaining 100 would be randomly sampled. - if (total < n_particles) then - do j = 1, int(n_particles/total) - index_temp = index_temp + 1 - temp_sites(index_temp) = fission_bank(i) - end do - end if - - ! Randomly sample sites needed - if (prn() < p_sample) then - index_temp = index_temp + 1 - temp_sites(index_temp) = fission_bank(i) - end if - end do - - ! At this point, the sampling of source sites is done and now we need to - ! figure out where to send source sites. Since it is possible that one - ! processor's share of the source bank spans more than just the immediate - ! neighboring processors, we have to perform an ALLGATHER to determine the - ! indices for all processors - -#ifdef OPENMC_MPI - ! First do an exclusive scan to get the starting indices for - start = 0_8 - call MPI_EXSCAN(index_temp, start, 1, MPI_INTEGER8, MPI_SUM, & - mpi_intracomm, mpi_err) - finish = start + index_temp - - ! Allocate space for bank_position if this hasn't been done yet - if (.not. allocated(bank_position)) allocate(bank_position(n_procs)) - call MPI_ALLGATHER(start, 1, MPI_INTEGER8, bank_position, 1, & - MPI_INTEGER8, mpi_intracomm, mpi_err) -#else - start = 0_8 - finish = index_temp -#endif - - ! Now that the sampling is complete, we need to ensure that we have exactly - ! n_particles source sites. The way this is done in a reproducible manner is - ! to adjust only the source sites on the last processor. - - if (rank == n_procs - 1) then - if (finish > n_particles) then - ! If we have extra sites sampled, we will simply discard the extra - ! ones on the last processor - index_temp = n_particles - start - - elseif (finish < n_particles) then - ! If we have too few sites, repeat sites from the very end of the - ! fission bank - sites_needed = n_particles - finish - do i = 1, int(sites_needed,4) - index_temp = index_temp + 1 - temp_sites(index_temp) = fission_bank(n_bank - sites_needed + i) - end do - end if - - ! the last processor should not be sending sites to right - finish = work_index(rank + 1) - end if - - call time_bank_sample % stop() - call time_bank_sendrecv % start() - -#ifdef OPENMC_MPI - ! ========================================================================== - ! SEND BANK SITES TO NEIGHBORS - - index_local = 1 - n_request = 0 - - if (start < n_particles) then - ! Determine the index of the processor which has the first part of the - ! source_bank for the local processor - neighbor = binary_search(work_index, n_procs + 1, start) - 1 - - SEND_SITES: do while (start < finish) - ! Determine the number of sites to send - n = min(work_index(neighbor + 1), finish) - start - - ! Initiate an asynchronous send of source sites to the neighboring - ! process - if (neighbor /= rank) then - n_request = n_request + 1 - call MPI_ISEND(temp_sites(index_local), int(n), MPI_BANK, neighbor, & - rank, mpi_intracomm, request(n_request), mpi_err) - end if - - ! Increment all indices - start = start + n - index_local = index_local + n - neighbor = neighbor + 1 - - ! Check for sites out of bounds -- this only happens in the rare - ! circumstance that a processor close to the end has so many sites that - ! it would exceed the bank on the last processor - if (neighbor > n_procs - 1) exit - end do SEND_SITES - end if - - ! ========================================================================== - ! RECEIVE BANK SITES FROM NEIGHBORS OR TEMPORARY BANK - - start = work_index(rank) - index_local = 1 - - ! Determine what process has the source sites that will need to be stored at - ! the beginning of this processor's source bank. - - if (start >= bank_position(n_procs)) then - neighbor = n_procs - 1 - else - neighbor = binary_search(bank_position, n_procs, start) - 1 - end if - - RECV_SITES: do while (start < work_index(rank + 1)) - ! Determine how many sites need to be received - if (neighbor == n_procs - 1) then - n = work_index(rank + 1) - start - else - n = min(bank_position(neighbor + 2), work_index(rank + 1)) - start - end if - - if (neighbor /= rank) then - ! If the source sites are not on this processor, initiate an - ! asynchronous receive for the source sites - - n_request = n_request + 1 - call MPI_IRECV(source_bank(index_local), int(n), MPI_BANK, & - neighbor, neighbor, mpi_intracomm, request(n_request), mpi_err) - - else - ! If the source sites are on this procesor, we can simply copy them - ! from the temp_sites bank - - index_temp = start - bank_position(rank+1) + 1 - source_bank(index_local:index_local+n-1) = & - temp_sites(index_temp:index_temp+n-1) - end if - - ! Increment all indices - start = start + n - index_local = index_local + n - neighbor = neighbor + 1 - end do RECV_SITES - - ! Since we initiated a series of asynchronous ISENDs and IRECVs, now we have - ! to ensure that the data has actually been communicated before moving on to - ! the next generation - - call MPI_WAITALL(n_request, request, MPI_STATUSES_IGNORE, mpi_err) - - ! Deallocate space for bank_position on the very last generation - if (current_batch == n_max_batches .and. current_gen == gen_per_batch) & - deallocate(bank_position) -#else - source_bank = temp_sites(1:n_particles) -#endif - - call time_bank_sendrecv % stop() - - ! Deallocate space for the temporary source bank on the last generation - if (current_batch == n_max_batches .and. current_gen == gen_per_batch) & - deallocate(temp_sites) - - end subroutine synchronize_bank - !=============================================================================== ! CALCULATE_GENERATION_KEFF collects the single-processor tracklength k's onto ! the master processor and normalizes them. This should work whether or not the diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index f03820492..459204ea3 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -9,9 +9,14 @@ #include "openmc/hdf5_interface.h" #include "openmc/mesh.h" #include "openmc/message_passing.h" +#include "openmc/random_lcg.h" +#include "openmc/search.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include // for min +#include + namespace openmc { //============================================================================== @@ -25,6 +30,245 @@ xt::xtensor source_frac; // Non-member functions //============================================================================== +void synchronize_bank() +{ + // Get pointers to source/fission bank + Bank* source_bank; + Bank* fission_bank; + int64_t n; + openmc_source_bank(&source_bank, &n); + openmc_fission_bank(&fission_bank, &n); + + // In order to properly understand the fission bank algorithm, you need to + // think of the fission and source bank as being one global array divided + // over multiple processors. At the start, each processor has a random amount + // of fission bank sites -- each processor needs to know the total number of + // sites in order to figure out the probability for selecting + // sites. Furthermore, each proc also needs to know where in the 'global' + // fission bank its own sites starts in order to ensure reproducibility by + // skipping ahead to the proper seed. + +#ifdef OPENMC_MPI + int64_t start = 0; + MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm); + + // While we would expect the value of start on rank 0 to be 0, the MPI + // standard says that the receive buffer on rank 0 is undefined and not + // significant + if (mpi::rank == 0) start = 0; + + int64_t finish = start + n_bank; + int64_t total = finish; + MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm); + +#else + int64_t start = 0; + int64_t finish = n_bank; + int64_t total = n_bank; +#endif + + // If there are not that many particles per generation, it's possible that no + // fission sites were created at all on a single processor. Rather than add + // extra logic to treat this circumstance, we really want to ensure the user + // runs enough particles to avoid this in the first place. + + if (n_bank == 0) { + fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank)); + } + + // Make sure all processors start at the same point for random sampling. Then + // skip ahead in the sequence using the starting index in the 'global' + // fission bank for each processor. + + set_particle_seed(simulation::total_gen + overall_generation()); + advance_prn_seed(start); + + // Determine how many fission sites we need to sample from the source bank + // and the probability for selecting a site. + + int64_t sites_needed; + if (total < settings::n_particles) { + sites_needed = settings::n_particles % total; + } else { + sites_needed = settings::n_particles; + } + double p_sample = static_cast(sites_needed) / total; + + //time_bank_sample % start() + + // ========================================================================== + // SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES + + // Allocate temporary source bank + int64_t index_temp = 0; + Bank temp_sites[3*simulation::work]; + + for (int64_t i = 0; i < n_bank; ++i) { + // If there are less than n_particles particles banked, automatically add + // int(n_particles/total) sites to temp_sites. For example, if you need + // 1000 and 300 were banked, this would add 3 source sites per banked site + // and the remaining 100 would be randomly sampled. + if (total < settings::n_particles) { + for (int64_t j = 1; j <= settings::n_particles / total; ++j) { + temp_sites[index_temp] = fission_bank[i]; + ++index_temp; + } + } + + // Randomly sample sites needed + if (prn() < p_sample) { + temp_sites[index_temp] = fission_bank[i]; + ++index_temp; + } + } + + // At this point, the sampling of source sites is done and now we need to + // figure out where to send source sites. Since it is possible that one + // processor's share of the source bank spans more than just the immediate + // neighboring processors, we have to perform an ALLGATHER to determine the + // indices for all processors + +#ifdef OPENMC_MPI + // First do an exclusive scan to get the starting indices for + start = 0; + MPI_Exscan(&index_temp, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm); + finish = start + index_temp; + + // Allocate space for bank_position if this hasn't been done yet + int64_t bank_position[mpi::n_procs]; + MPI_Allgather(&start, 1, MPI_INT64_T, bank_position, 1, + MPI_INT64_T, mpi::intracomm); +#else + start = 0; + finish = index_temp; +#endif + + // Now that the sampling is complete, we need to ensure that we have exactly + // n_particles source sites. The way this is done in a reproducible manner is + // to adjust only the source sites on the last processor. + + if (mpi::rank == mpi::n_procs - 1) { + if (finish > settings::n_particles) { + // If we have extra sites sampled, we will simply discard the extra + // ones on the last processor + index_temp = settings::n_particles - start; + + } else if (finish < settings::n_particles) { + // If we have too few sites, repeat sites from the very end of the + // fission bank + sites_needed = settings::n_particles - finish; + for (int i = 0; i < sites_needed; ++i) { + temp_sites[index_temp] = fission_bank[n_bank - sites_needed + i]; + ++index_temp; + } + } + + // the last processor should not be sending sites to right + finish = simulation::work_index[mpi::rank + 1]; + } + + //time_bank_sample % stop() + //time_bank_sendrecv % start() + +#ifdef OPENMC_MPI + // ========================================================================== + // SEND BANK SITES TO NEIGHBORS + + int64_t index_local = 0; + std::vector requests; + + if (start < settings::n_particles) { + // Determine the index of the processor which has the first part of the + // source_bank for the local processor + int neighbor = upper_bound_index(simulation::work_index.begin(), + simulation::work_index.end(), start); + + while (start < finish) { + // Determine the number of sites to send + int64_t n = std::min(simulation::work_index[neighbor + 1], finish) - start; + + // Initiate an asynchronous send of source sites to the neighboring + // process + if (neighbor != mpi::rank) { + requests.emplace_back(); + MPI_Isend(&temp_sites[index_local], static_cast(n), mpi::bank, + neighbor, mpi::rank, mpi::intracomm, &requests.back()); + } + + // Increment all indices + start += n; + index_local += n; + ++neighbor; + + // Check for sites out of bounds -- this only happens in the rare + // circumstance that a processor close to the end has so many sites that + // it would exceed the bank on the last processor + if (neighbor > mpi::n_procs - 1) break; + } + } + + // ========================================================================== + // RECEIVE BANK SITES FROM NEIGHBORS OR TEMPORARY BANK + + start = simulation::work_index[mpi::rank]; + index_local = 0; + + // Determine what process has the source sites that will need to be stored at + // the beginning of this processor's source bank. + + int neighbor; + if (start >= bank_position[mpi::n_procs - 1]) { + neighbor = mpi::n_procs - 1; + } else { + neighbor = upper_bound_index(bank_position, bank_position + mpi::n_procs, start); + } + + while (start < simulation::work_index[mpi::rank + 1]) { + // Determine how many sites need to be received + int64_t n; + if (neighbor == mpi::n_procs - 1) { + n = simulation::work_index[mpi::rank + 1] - start; + } else { + n = std::min(bank_position[neighbor + 1], simulation::work_index[mpi::rank + 1]) - start; + } + + if (neighbor != mpi::rank) { + // If the source sites are not on this processor, initiate an + // asynchronous receive for the source sites + + requests.emplace_back(); + MPI_Irecv(&source_bank[index_local], static_cast(n), mpi::bank, + neighbor, neighbor, mpi::intracomm, &requests.back()); + + } else { + // If the source sites are on this procesor, we can simply copy them + // from the temp_sites bank + + index_temp = start - bank_position[mpi::rank]; + std::copy(&temp_sites[index_temp], &temp_sites[index_temp + n], + &source_bank[index_local]); + } + + // Increment all indices + start += n; + index_local += n; + ++neighbor; + } + + // Since we initiated a series of asynchronous ISENDs and IRECVs, now we have + // to ensure that the data has actually been communicated before moving on to + // the next generation + + int n_request = requests.size(); + MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE); + +#else + std::copy(temp_sites, temp_sites + settings::n_particles, source_bank); +#endif + + //time_bank_sendrecv % stop() +} + void shannon_entropy() { // Get pointer to entropy mesh diff --git a/src/initialize.F90 b/src/initialize.F90 index 0dbaa83bf..697fc1235 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -83,8 +83,9 @@ contains call time_initialize%start() #ifdef OPENMC_MPI - ! Setup MPI - call initialize_mpi(comm) + ! Indicate that MPI is turned on + mpi_enabled = .true. + mpi_intracomm = intracomm #endif #ifdef _OPENMP @@ -114,57 +115,6 @@ contains err = 0 end function openmc_init_f -#ifdef OPENMC_MPI -!=============================================================================== -! INITIALIZE_MPI starts up the Message Passing Interface (MPI) and determines -! the number of processors the problem is being run with as well as the rank of -! each processor. -!=============================================================================== - - subroutine initialize_mpi(intracomm) -#ifdef OPENMC_MPIF08 - type(MPI_Comm), intent(in) :: intracomm ! MPI intracommunicator -#else - integer, intent(in) :: intracomm ! MPI intracommunicator -#endif - - integer :: mpi_err ! MPI error code - integer :: bank_blocks(5) ! Count for each datatype -#ifdef OPENMC_MPIF08 - type(MPI_Datatype) :: bank_types(5) -#else - integer :: bank_types(5) ! Datatypes -#endif - integer(MPI_ADDRESS_KIND) :: bank_disp(5) ! Displacements - type(Bank) :: b - - ! Indicate that MPI is turned on - mpi_enabled = .true. - mpi_intracomm = intracomm - - ! ========================================================================== - ! CREATE MPI_BANK TYPE - - ! Determine displacements for MPI_BANK type - call MPI_GET_ADDRESS(b % wgt, bank_disp(1), mpi_err) - call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err) - call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err) - call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err) - call MPI_GET_ADDRESS(b % delayed_group, bank_disp(5), mpi_err) - - ! Adjust displacements - bank_disp = bank_disp - bank_disp(1) - - ! Define MPI_BANK for fission sites - bank_blocks = (/ 1, 3, 3, 1, 1 /) - bank_types = (/ MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_INTEGER /) - call MPI_TYPE_CREATE_STRUCT(5, bank_blocks, bank_disp, & - bank_types, MPI_BANK, mpi_err) - call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) - - end subroutine initialize_mpi -#endif - !=============================================================================== ! READ_COMMAND_LINE reads all parameters from the command line !=============================================================================== diff --git a/src/initialize.cpp b/src/initialize.cpp index ac7a5ef15..31fa8fb25 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -80,16 +80,18 @@ void initialize_mpi(MPI_Comm intracomm) // Create bank datatype Bank b; - MPI_Aint disp[] { - offsetof(Bank, wgt), - offsetof(Bank, xyz), - offsetof(Bank, uvw), - offsetof(Bank, E), - offsetof(Bank, delayed_group) - }; - int blocks[] {1, 3, 3, 1, 1}; - MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT}; - MPI_Type_create_struct(5, blocks, disp, types, &mpi::bank); + MPI_Aint disp[6]; + MPI_Get_address(&b.wgt, &disp[0]); + MPI_Get_address(&b.xyz, &disp[1]); + MPI_Get_address(&b.uvw, &disp[2]); + MPI_Get_address(&b.E, &disp[3]); + MPI_Get_address(&b.delayed_group, &disp[4]); + MPI_Get_address(&b.particle, &disp[5]); + for (int i = 5; i >= 0; --i) disp[i] -= disp[0]; + + int blocks[] {1, 3, 3, 1, 1, 1}; + MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT}; + MPI_Type_create_struct(6, blocks, disp, types, &mpi::bank); MPI_Type_commit(&mpi::bank); } #endif // OPENMC_MPI diff --git a/src/message_passing.F90 b/src/message_passing.F90 index 6ade0895a..05219f0f8 100644 --- a/src/message_passing.F90 +++ b/src/message_passing.F90 @@ -19,10 +19,8 @@ module message_passing logical(C_BOOL), bind(C, name='openmc_master') :: master = .true. ! master process? logical :: mpi_enabled = .false. ! is MPI in use and initialized? #ifdef OPENMC_MPIF08 - type(MPI_Datatype) :: MPI_BANK ! MPI datatype for fission bank type(MPI_Comm) :: mpi_intracomm ! MPI intra-communicator #else - integer :: MPI_BANK ! MPI datatype for fission bank integer :: mpi_intracomm ! MPI intra-communicator #endif diff --git a/src/simulation.F90 b/src/simulation.F90 index 767ec5244..2359306d2 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -11,7 +11,7 @@ module simulation use cmfd_header, only: cmfd_on use constants, only: ZERO use eigenvalue, only: calculate_average_keff, calculate_generation_keff, & - synchronize_bank, keff_generation, k_sum + keff_generation, k_sum #ifdef _OPENMP use eigenvalue, only: join_bank_from_threads #endif @@ -258,6 +258,9 @@ contains subroutine shannon_entropy() bind(C) end subroutine + + subroutine synchronize_bank() bind(C) + end subroutine end interface ! Update global tallies with the omp private accumulation variables diff --git a/src/simulation.cpp b/src/simulation.cpp index 386e00510..2935d71ac 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -86,8 +86,8 @@ void calculate_work() int64_t remainder = settings::n_particles % mpi::n_procs; int64_t i_bank = 0; - simulation::work_index.reserve(mpi::n_procs); - simulation::work_index.push_back(0); + simulation::work_index.resize(mpi::n_procs + 1); + simulation::work_index[0] = 0; for (int i = 0; i < mpi::n_procs; ++i) { // Number of particles for rank i int64_t work_i = i < remainder ? min_work + 1 : min_work; @@ -97,7 +97,7 @@ void calculate_work() // Set index into source bank for rank i i_bank += work_i; - simulation::work_index.push_back(i_bank); + simulation::work_index[i + 1] = i_bank; } } From d7da87ca1d98e6dc867e4af7c38dc9f56b031349 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 12 Oct 2018 07:03:20 -0500 Subject: [PATCH 15/23] Write C++ version of Timer class and use in synchronize_bank --- CMakeLists.txt | 1 + include/openmc/timer.h | 47 +++++++++++++++++++++++++++++++++ src/api.F90 | 4 +-- src/eigenvalue.cpp | 12 ++++++--- src/output.F90 | 6 ++--- src/simulation.F90 | 2 -- src/state_point.F90 | 6 ++--- src/timer.cpp | 60 ++++++++++++++++++++++++++++++++++++++++++ src/timer_header.F90 | 22 +++++++++++++--- 9 files changed, 142 insertions(+), 18 deletions(-) create mode 100644 include/openmc/timer.h create mode 100644 src/timer.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 99d4e4a5e..84f21427a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -443,6 +443,7 @@ add_library(libopenmc SHARED src/string_functions.cpp src/summary.cpp src/surface.cpp + src/timer.cpp src/thermal.cpp src/xml_interface.cpp src/xsdata.cpp diff --git a/include/openmc/timer.h b/include/openmc/timer.h new file mode 100644 index 000000000..11dfd1adf --- /dev/null +++ b/include/openmc/timer.h @@ -0,0 +1,47 @@ +#ifndef OPENMC_TIMER_H +#define OPENMC_TIMER_H + +#include + +namespace openmc { + +//============================================================================== +//! Class for measuring time elapsed +//============================================================================== + +class Timer { +public: + using clock = std::chrono::high_resolution_clock; + + Timer() {}; + + //! Start running the timer + void start (); + + //! Get total elapsed time in seconds + //! \return Elapsed time in [s] + double elapsed(); + + //! Stop running the timer + void stop(); + + //! Stop the timer and reset its elapsed time + void reset(); + +private: + bool running_ {false}; //!< is timer running? + std::chrono::time_point start_; //!< starting point for clock + double elapsed_; //!< elasped time in [s] +}; + +//============================================================================== +// Global variables +//============================================================================== + +extern Timer time_bank; +extern Timer time_bank_sample; +extern Timer time_bank_sendrecv; + +} // namespace openmc + +#endif // OPENMC_TIMER_H diff --git a/src/api.F90 b/src/api.F90 index d79aa2450..18162e600 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -276,14 +276,12 @@ contains call time_initialize % reset() call time_read_xs % reset() call time_unionize % reset() - call time_bank % reset() - call time_bank_sample % reset() - call time_bank_sendrecv % reset() call time_tallies % reset() call time_inactive % reset() call time_active % reset() call time_transport % reset() call time_finalize % reset() + call reset_timers() err = 0 end function openmc_reset diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 459204ea3..ba3ebfae5 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -13,6 +13,7 @@ #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include "openmc/timer.h" #include // for min #include @@ -32,6 +33,8 @@ xt::xtensor source_frac; void synchronize_bank() { + time_bank.start(); + // Get pointers to source/fission bank Bank* source_bank; Bank* fission_bank; @@ -94,7 +97,7 @@ void synchronize_bank() } double p_sample = static_cast(sites_needed) / total; - //time_bank_sample % start() + time_bank_sample.start(); // ========================================================================== // SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES @@ -167,8 +170,8 @@ void synchronize_bank() finish = simulation::work_index[mpi::rank + 1]; } - //time_bank_sample % stop() - //time_bank_sendrecv % start() + time_bank_sample.stop(); + time_bank_sendrecv.start(); #ifdef OPENMC_MPI // ========================================================================== @@ -266,7 +269,8 @@ void synchronize_bank() std::copy(temp_sites, temp_sites + settings::n_particles, source_bank); #endif - //time_bank_sendrecv % stop() + time_bank_sendrecv.stop(); + time_bank.stop(); } void shannon_entropy() diff --git a/src/output.F90 b/src/output.F90 index 4363731f0..051167511 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -511,9 +511,9 @@ contains end if write(ou,100) " Time in active batches", time_active % elapsed if (run_mode == MODE_EIGENVALUE) then - write(ou,100) " Time synchronizing fission bank", time_bank % elapsed - write(ou,100) " Sampling source sites", time_bank_sample % elapsed - write(ou,100) " SEND/RECV source sites", time_bank_sendrecv % elapsed + write(ou,100) " Time synchronizing fission bank", time_bank_elapsed() + write(ou,100) " Sampling source sites", time_bank_sample_elapsed() + write(ou,100) " SEND/RECV source sites", time_bank_sendrecv_elapsed() end if write(ou,100) " Time accumulating tallies", time_tallies % elapsed if (cmfd_run) write(ou,100) " Time in CMFD", time_cmfd % elapsed diff --git a/src/simulation.F90 b/src/simulation.F90 index 2359306d2..4bbf94ed1 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -294,9 +294,7 @@ contains #endif ! Distribute fission bank across processors evenly - call time_bank % start() call synchronize_bank() - call time_bank % stop() ! Calculate shannon entropy if (entropy_on) call shannon_entropy() diff --git a/src/state_point.F90 b/src/state_point.F90 index 11e4ad621..548db766c 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -419,11 +419,11 @@ contains time_active % get_value()) if (run_mode == MODE_EIGENVALUE) then call write_dataset(runtime_group, "synchronizing fission bank", & - time_bank % get_value()) + time_bank_elapsed()) call write_dataset(runtime_group, "sampling source sites", & - time_bank_sample % get_value()) + time_bank_sample_elapsed()) call write_dataset(runtime_group, "SEND-RECV source sites", & - time_bank_sendrecv % get_value()) + time_bank_sendrecv_elapsed()) end if call write_dataset(runtime_group, "accumulating tallies", & time_tallies % get_value()) diff --git a/src/timer.cpp b/src/timer.cpp new file mode 100644 index 000000000..85b1eb771 --- /dev/null +++ b/src/timer.cpp @@ -0,0 +1,60 @@ +#include "openmc/timer.h" + +namespace openmc { + +//============================================================================== +// Timer implementation +//============================================================================== + +void Timer::start () +{ + running_ = true; + start_ = clock::now(); +} + +void Timer::stop() +{ + elapsed_ = elapsed(); + running_ = false; +} + +void Timer::reset() +{ + running_ = false; + elapsed_ = 0.0; +} + +double Timer::elapsed() +{ + if (running_) { + std::chrono::duration diff = clock::now() - start_; + return elapsed_ += diff.count(); + } else { + return elapsed_; + } +} + +//============================================================================== +// Global variables +//============================================================================== + +Timer time_bank; +Timer time_bank_sample; +Timer time_bank_sendrecv; + +//============================================================================== +// Fortran compatibility +//============================================================================== + +extern "C" double time_bank_elapsed() { return time_bank.elapsed(); } +extern "C" double time_bank_sample_elapsed() { return time_bank_sample.elapsed(); } +extern "C" double time_bank_sendrecv_elapsed() { return time_bank_sendrecv.elapsed(); } + +extern "C" void reset_timers() +{ + time_bank.reset(); + time_bank_sample.reset(); + time_bank_sendrecv.reset(); +} + +} // namespace openmc diff --git a/src/timer_header.F90 b/src/timer_header.F90 index a147a551d..169154696 100644 --- a/src/timer_header.F90 +++ b/src/timer_header.F90 @@ -1,9 +1,28 @@ module timer_header + use, intrinsic :: ISO_C_BINDING + use constants, only: ZERO implicit none + interface + function time_bank_elapsed() result(t) bind(C) + import C_DOUBLE + real(C_DOUBLE) :: t + end function + function time_bank_sample_elapsed() result(t) bind(C) + import C_DOUBLE + real(C_DOUBLE) :: t + end function + function time_bank_sendrecv_elapsed() result(t) bind(C) + import C_DOUBLE + real(C_DOUBLE) :: t + end function + subroutine reset_timers() bind(C) + end subroutine + end interface + !=============================================================================== ! TIMER represents a timer that can be started and stopped to measure how long ! different routines run. The intrinsic routine system_clock is used to measure @@ -29,9 +48,6 @@ module timer_header type(Timer) :: time_initialize ! timer for initialization type(Timer) :: time_read_xs ! timer for reading cross sections type(Timer) :: time_unionize ! timer for material xs-energy grid union - type(Timer) :: time_bank ! timer for fission bank synchronization - type(Timer) :: time_bank_sample ! timer for fission bank sampling - type(Timer) :: time_bank_sendrecv ! timer for fission bank SEND/RECV type(Timer) :: time_tallies ! timer for accumulate tallies type(Timer) :: time_inactive ! timer for inactive batches type(Timer) :: time_active ! timer for active batches From 8ebd26b160cca2b9f0382b41f0a4ef188aa4aae6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 12 Oct 2018 08:01:23 -0500 Subject: [PATCH 16/23] Move code related to k_generation to C++ --- include/openmc/constants.h | 8 +++---- include/openmc/eigenvalue.h | 4 ++++ include/openmc/simulation.h | 3 +++ src/eigenvalue.F90 | 43 +++++++---------------------------- src/eigenvalue.cpp | 45 +++++++++++++++++++++++++++++++++++-- src/input_xml.F90 | 2 +- src/output.F90 | 4 ++-- src/simulation.F90 | 31 +++++-------------------- src/simulation.cpp | 26 ++++++++++++++++++++- src/simulation_header.F90 | 23 ++++++++++++++++--- src/state_point.F90 | 39 +++++++++----------------------- 11 files changed, 126 insertions(+), 102 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 6ead514d7..361cfc158 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -408,10 +408,10 @@ constexpr int STANDARD_DEVIATION {3}; // Global tally parameters constexpr int N_GLOBAL_TALLIES {4}; -constexpr int K_COLLISION {1}; -constexpr int K_ABSORPTION {2}; -constexpr int K_TRACKLENGTH {3}; -constexpr int LEAKAGE {4}; +constexpr int K_COLLISION {0}; +constexpr int K_ABSORPTION {1}; +constexpr int K_TRACKLENGTH {2}; +constexpr int LEAKAGE {3}; // Differential tally independent variables constexpr int DIFF_DENSITY {1}; diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index 59b448e29..cbf0baa2d 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -14,6 +14,7 @@ namespace openmc { // Global variables //============================================================================== +extern double keff_generation; //!< Single-generation k on each processor extern std::vector entropy; //!< Shannon entropy at each generation extern xt::xtensor source_frac; //!< Source fraction for UFS @@ -24,6 +25,9 @@ extern "C" int64_t n_bank; // Non-member functions //============================================================================== +//! Collect/normalize the tracklength keff from each process +extern "C" void calculate_generation_keff(); + //! Sample/redistribute source sites from accumulated fission sites extern "C" void synchronize_bank(); diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index ceb5eb92a..e21ea05d5 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -56,6 +56,9 @@ void calculate_work(); //! Initialize simulation extern "C" void openmc_simulation_init_c(); +//! Initialize a fission generation +extern "C" void initialize_generation(); + //! Determine overall generation number extern "C" int overall_generation(); diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index df108f6d7..92c731eae 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -16,42 +16,15 @@ module eigenvalue implicit none - real(8) :: keff_generation ! Single-generation k on each processor real(8) :: k_sum(2) = ZERO ! Used to reduce sum and sum_sq + interface + subroutine calculate_generation_keff() bind(C) + end subroutine + end interface + contains -!=============================================================================== -! CALCULATE_GENERATION_KEFF collects the single-processor tracklength k's onto -! the master processor and normalizes them. This should work whether or not the -! no-reduce method is being used. -!=============================================================================== - - subroutine calculate_generation_keff() - - real(8) :: keff_reduced -#ifdef OPENMC_MPI - integer :: mpi_err ! MPI error code -#endif - - ! Get keff for this generation by subtracting off the starting value - keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - keff_generation - -#ifdef OPENMC_MPI - ! Combine values across all processors - call MPI_ALLREDUCE(keff_generation, keff_reduced, 1, MPI_REAL8, & - MPI_SUM, mpi_intracomm, mpi_err) -#else - keff_reduced = keff_generation -#endif - - ! Normalize single batch estimate of k - ! TODO: This should be normalized by total_weight, not by n_particles - keff_reduced = keff_reduced / n_particles - call k_generation % push_back(keff_reduced) - - end subroutine calculate_generation_keff - !=============================================================================== ! CALCULATE_AVERAGE_KEFF calculates the mean and standard deviation of the mean ! of k-effective during active generations and broadcasts the mean to all @@ -76,12 +49,12 @@ contains if (n <= 0) then ! For inactive generations, use current generation k as estimate for next ! generation - keff = k_generation % data(i) + keff = k_generation(i) else ! Sample mean of keff - k_sum(1) = k_sum(1) + k_generation % data(i) - k_sum(2) = k_sum(2) + k_generation % data(i)**2 + k_sum(1) = k_sum(1) + k_generation(i) + k_sum(2) = k_sum(2) + k_generation(i)**2 ! Determine mean keff = k_sum(1) / n diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index ba3ebfae5..a9f99a946 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -5,6 +5,7 @@ #include "xtensor/xview.hpp" #include "openmc/capi.h" +#include "openmc/constants.h" #include "openmc/error.h" #include "openmc/hdf5_interface.h" #include "openmc/mesh.h" @@ -14,6 +15,7 @@ #include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/timer.h" +#include "openmc/tallies/tally.h" #include // for min #include @@ -24,6 +26,7 @@ namespace openmc { // Global variables //============================================================================== +double keff_generation; std::vector entropy; xt::xtensor source_frac; @@ -31,6 +34,28 @@ xt::xtensor source_frac; // Non-member functions //============================================================================== +void calculate_generation_keff() +{ + auto gt = global_tallies(); + + // Get keff for this generation by subtracting off the starting value + keff_generation = gt(K_TRACKLENGTH, RESULT_VALUE) - keff_generation; + + double keff_reduced; +#ifdef OPENMC_MPI + // Combine values across all processors + MPI_Allreduce(&keff_generation, &keff_reduced, 1, MPI_DOUBLE, + MPI_SUM, mpi::intracomm); +#else + keff_reduced = keff_generation; +#endif + + // Normalize single batch estimate of k + // TODO: This should be normalized by total_weight, not by n_particles + keff_reduced /= settings::n_particles; + simulation::k_generation.push_back(keff_reduced); +} + void synchronize_bank() { time_bank.start(); @@ -375,18 +400,34 @@ double ufs_get_weight(const Particle* p) } } -extern "C" void entropy_to_hdf5(hid_t group) +extern "C" void write_eigenvalue_hdf5(hid_t group) { + write_dataset(group, "n_inactive", settings::n_inactive); + write_dataset(group, "generations_per_batch", settings::gen_per_batch); + write_dataset(group, "k_generation", simulation::k_generation); if (settings::entropy_on) { write_dataset(group, "entropy", entropy); } + write_dataset(group, "k_col_abs", simulation::k_col_abs); + write_dataset(group, "k_col_tra", simulation::k_col_tra); + write_dataset(group, "k_abs_tra", simulation::k_abs_tra); + std::array k_combined; + openmc_get_keff(k_combined.data()); + write_dataset(group, "k_combined", k_combined); } -extern "C" void entropy_from_hdf5(hid_t group) +extern "C" void read_eigenvalue_hdf5(hid_t group) { + read_dataset(group, "generations_per_batch", settings::gen_per_batch); + int n = simulation::restart_batch*settings::gen_per_batch; + simulation::k_generation.resize(n); + read_dataset(group, "k_generation", simulation::k_generation); if (settings::entropy_on) { read_dataset(group, "entropy", entropy); } + read_dataset(group, "k_col_abs", simulation::k_col_abs); + read_dataset(group, "k_col_tra", simulation::k_col_tra); + read_dataset(group, "k_abs_tra", simulation::k_abs_tra); } extern "C" double entropy_c(int i) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 989f7b357..81f323351 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -223,7 +223,7 @@ contains if (run_mode == MODE_EIGENVALUE) then ! Preallocate space for keff and entropy by generation - call k_generation % reserve(n_max_batches*gen_per_batch) + call k_generation_reserve(n_max_batches*gen_per_batch) end if ! Particle tracks diff --git a/src/output.F90 b/src/output.F90 index 051167511..d775a2da7 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -339,7 +339,7 @@ contains ! write out information about batch and generation write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') & trim(to_str(current_batch)) // "/" // trim(to_str(current_gen)) - write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation % data(i) + write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation(i) ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & @@ -373,7 +373,7 @@ contains write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') & trim(to_str(current_batch)) // "/" // trim(to_str(gen_per_batch)) write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') & - k_generation % data(i) + k_generation(i) ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & diff --git a/src/simulation.F90 b/src/simulation.F90 index 4bbf94ed1..b94166f87 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -11,7 +11,7 @@ module simulation use cmfd_header, only: cmfd_on use constants, only: ZERO use eigenvalue, only: calculate_average_keff, calculate_generation_keff, & - keff_generation, k_sum + k_sum #ifdef _OPENMP use eigenvalue, only: join_bank_from_threads #endif @@ -57,6 +57,9 @@ module simulation subroutine initialize_source() bind(C) end subroutine + subroutine initialize_generation() bind(C) + end subroutine + function sample_external_source() result(site) bind(C) import Bank type(Bank) :: site @@ -222,30 +225,6 @@ contains end subroutine initialize_batch -!=============================================================================== -! INITIALIZE_GENERATION -!=============================================================================== - - subroutine initialize_generation() - - interface - subroutine ufs_count_sites() bind(C) - end subroutine - end interface - - if (run_mode == MODE_EIGENVALUE) then - ! Reset number of fission bank sites - n_bank = 0 - - ! Count source sites if using uniform fission source weighting - if (ufs) call ufs_count_sites() - - ! Store current value of tracklength k - keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - end if - - end subroutine initialize_generation - !=============================================================================== ! FINALIZE_GENERATION !=============================================================================== @@ -438,7 +417,7 @@ contains ! Reset global variables -- this is done before loading state point (as that ! will potentially populate k_generation and entropy) current_batch = 0 - call k_generation % clear() + call k_generation_clear() call entropy_clear() need_depletion_rx = .false. diff --git a/src/simulation.cpp b/src/simulation.cpp index 2935d71ac..0d94966ed 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1,6 +1,7 @@ #include "openmc/simulation.h" #include "openmc/capi.h" +#include "openmc/eigenvalue.h" #include "openmc/message_passing.h" #include "openmc/settings.h" #include "openmc/tallies/tally.h" @@ -62,7 +63,7 @@ int thread_id; //!< ID of a given thread } // namespace simulation //============================================================================== -// Functions +// Non-member functions //============================================================================== void openmc_simulation_init_c() @@ -71,6 +72,20 @@ void openmc_simulation_init_c() calculate_work(); } +void initialize_generation() +{ + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + // Reset number of fission bank sites + n_bank = 0; + + // Count source sites if using uniform fission source weighting + if (settings::ufs_on) ufs_count_sites(); + + // Store current value of tracklength k + keff_generation = global_tallies()(K_TRACKLENGTH, RESULT_VALUE); + } +} + int overall_generation() { using namespace simulation; @@ -139,4 +154,13 @@ void broadcast_triggers() } #endif +//============================================================================== +// Fortran compatibility +//============================================================================== + +extern "C" double k_generation(int i) { return simulation::k_generation.at(i - 1); } +extern "C" int k_generation_size() { return simulation::k_generation.size(); } +extern "C" void k_generation_clear() { simulation::k_generation.clear(); } +extern "C" void k_generation_reserve(int i) { simulation::k_generation.reserve(i); } + } // namespace openmc diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 29ee34ccf..65a36654f 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -39,7 +39,6 @@ module simulation_header ! K-EIGENVALUE SIMULATION VARIABLES ! Temporary k-effective values - type(VectorReal) :: k_generation ! single-generation estimates of k real(C_DOUBLE), bind(C) :: keff ! average k over active batches real(C_DOUBLE), bind(C) :: keff_std ! standard deviation of average k real(C_DOUBLE), bind(C) :: k_col_abs ! sum over batches of k_collision * k_absorption @@ -71,6 +70,25 @@ module simulation_header import C_INT integer(C_INT) :: gen end function overall_generation + + function k_generation(i) result(k) bind(C) + import C_DOUBLE, C_INT + integer(C_INT), value :: i + real(C_DOUBLE) :: k + end function + + function k_generation_size() result(sz) bind(C) + import C_INT + integer(C_INT) :: sz + end function + + subroutine k_generation_clear() bind(C) + end subroutine + + subroutine k_generation_reserve(i) bind(C) + import C_INT + integer(C_INT), value :: i + end subroutine end interface contains @@ -83,8 +101,7 @@ contains if (allocated(work_index)) deallocate(work_index) - call k_generation % clear() - call k_generation % shrink_to_fit() + call k_generation_clear() call entropy_clear() end subroutine free_memory_simulation diff --git a/src/state_point.F90 b/src/state_point.F90 index 548db766c..f861f871a 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -63,15 +63,13 @@ contains integer(C_INT) :: err logical :: write_source_ - integer :: i, j, k + integer :: i, j integer :: i_xs integer, allocatable :: id_array(:) integer(HID_T) :: file_id integer(HID_T) :: cmfd_group, tallies_group, tally_group, & filters_group, filter_group, derivs_group, & deriv_group, runtime_group - integer(C_INT) :: ignored_err - real(C_DOUBLE) :: k_combined(2) character(MAX_WORD_LEN), allocatable :: str_array(:) character(C_CHAR), pointer :: string(:) character(len=:, kind=C_CHAR), allocatable :: filename_ @@ -83,7 +81,7 @@ contains import HID_T integer(HID_T), value :: group end subroutine - subroutine entropy_to_hdf5(group) bind(C) + subroutine write_eigenvalue_hdf5(group) bind(C) import HID_T integer(HID_T), value :: group end subroutine @@ -173,16 +171,7 @@ contains ! Write out information for eigenvalue run if (run_mode == MODE_EIGENVALUE) then - call write_dataset(file_id, "n_inactive", n_inactive) - call write_dataset(file_id, "generations_per_batch", gen_per_batch) - k = k_generation % size() - call write_dataset(file_id, "k_generation", k_generation % data(1:k)) - call entropy_to_hdf5(file_id) - call write_dataset(file_id, "k_col_abs", k_col_abs) - call write_dataset(file_id, "k_col_tra", k_col_tra) - call write_dataset(file_id, "k_abs_tra", k_abs_tra) - ignored_err = openmc_get_keff(k_combined) - call write_dataset(file_id, "k_combined", k_combined) + call write_eigenvalue_hdf5(file_id) ! Write out CMFD info if (cmfd_on) then @@ -513,7 +502,9 @@ contains character(MAX_WORD_LEN) :: word interface - subroutine entropy_from_hdf5() bind(C) + subroutine read_eigenvalue_hdf5(group) bind(C) + import HID_T + integer(HID_T), value :: group end subroutine end interface @@ -592,16 +583,7 @@ contains ! Read information specific to eigenvalue run if (run_mode == MODE_EIGENVALUE) then call read_dataset(int_array(1), file_id, "n_inactive") - call read_dataset(gen_per_batch, file_id, "generations_per_batch") - - n = restart_batch*gen_per_batch - call k_generation % resize(n) - call read_dataset(k_generation % data(1:n), file_id, "k_generation") - - call entropy_from_hdf5() - call read_dataset(k_col_abs, file_id, "k_col_abs") - call read_dataset(k_col_tra, file_id, "k_col_tra") - call read_dataset(k_abs_tra, file_id, "k_abs_tra") + call read_eigenvalue_hdf5(file_id) ! Take maximum of statepoint n_inactive and input n_inactive n_inactive = max(n_inactive, int_array(1)) @@ -634,13 +616,14 @@ contains ! of active cycle or inactive cycle if (restart_batch > n_inactive) then do i = n_inactive + 1, restart_batch - k_sum(1) = k_sum(1) + k_generation % data(i) - k_sum(2) = k_sum(2) + k_generation % data(i)**2 + k_sum(1) = k_sum(1) + k_generation(i) + k_sum(2) = k_sum(2) + k_generation(i)**2 end do n = gen_per_batch*n_realizations keff = k_sum(1) / n else - keff = k_generation % data(n) + n = k_generation_size() + keff = k_generation(n) end if current_batch = restart_batch From bbccbf436d2523feec7e6454290dbebd380f1e69 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 14 Oct 2018 16:05:12 -0500 Subject: [PATCH 17/23] Use namespace in finalize.h/cpp --- include/openmc/finalize.h | 4 ++++ src/finalize.cpp | 6 +++++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/include/openmc/finalize.h b/include/openmc/finalize.h index 58cf75606..5a8238841 100644 --- a/include/openmc/finalize.h +++ b/include/openmc/finalize.h @@ -1,6 +1,10 @@ #ifndef OPENMC_FINALIZE_H #define OPENMC_FINALIZE_H +namespace openmc { + extern "C" void openmc_free_bank(); +} // namespace openmc + #endif // OPENMC_FINALIZE_H diff --git a/src/finalize.cpp b/src/finalize.cpp index 3b4f26764..7f07ddac0 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -2,9 +2,13 @@ #include "openmc/message_passing.h" +namespace openmc { + void openmc_free_bank() { #ifdef OPENMC_MPI - MPI_Type_free(&openmc::mpi::bank); + MPI_Type_free(&mpi::bank); #endif } + +} // namespace openmc From 242895e99bbe9f76e050146d8a38f00230c05e40 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 15 Oct 2018 07:45:26 -0500 Subject: [PATCH 18/23] Remove work_index on Fortran side --- include/openmc/state_point.h | 6 ++--- src/simulation.F90 | 43 ------------------------------------ src/simulation.cpp | 1 + src/simulation_header.F90 | 9 +++++--- src/source.cpp | 4 ++-- src/state_point.F90 | 12 +++++----- src/state_point.cpp | 25 +++++++++++---------- 7 files changed, 29 insertions(+), 71 deletions(-) diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index 6541a67fc..058edc09c 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -9,10 +9,8 @@ namespace openmc { -extern "C" void write_source_bank(hid_t group_id, int64_t* work_index, - Bank* source_bank); -extern "C" void read_source_bank(hid_t group_id, int64_t* work_index, - Bank* source_bank); +extern "C" void write_source_bank(hid_t group_id, Bank* source_bank); +extern "C" void read_source_bank(hid_t group_id, Bank* source_bank); extern "C" void write_tally_results_nr(hid_t file_id); diff --git a/src/simulation.F90 b/src/simulation.F90 index b94166f87..0d6c534f1 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -383,9 +383,6 @@ contains ! Set up tally procedure pointers call init_tally_routines() - ! Determine how much work each processor should do - call calculate_work() - ! Allocate source bank, and for eigenvalue simulations also allocate the ! fission bank call allocate_banks() @@ -516,46 +513,6 @@ contains end function openmc_simulation_finalize -!=============================================================================== -! CALCULATE_WORK determines how many particles each processor should simulate -!=============================================================================== - - subroutine calculate_work() - - integer :: i ! loop index - integer :: remainder ! Number of processors with one extra particle - integer(8) :: i_bank ! Running count of number of particles - integer(8) :: min_work ! Minimum number of particles on each proc - integer(8) :: work_i ! Number of particles on rank i - - if (.not. allocated(work_index)) allocate(work_index(0:n_procs)) - - ! Determine minimum amount of particles to simulate on each processor - min_work = n_particles/n_procs - - ! Determine number of processors that have one extra particle - remainder = int(mod(n_particles, int(n_procs,8)), 4) - - i_bank = 0 - work_index(0) = 0 - do i = 0, n_procs - 1 - ! Number of particles for rank i - if (i < remainder) then - work_i = min_work + 1 - else - work_i = min_work - end if - - ! Set number of particles - if (rank == i) work = work_i - - ! Set index into source bank for rank i - i_bank = i_bank + work_i - work_index(i+1) = i_bank - end do - - end subroutine calculate_work - !=============================================================================== ! ALLOCATE_BANKS allocates memory for the fission and source banks !=============================================================================== diff --git a/src/simulation.cpp b/src/simulation.cpp index 0d94966ed..13fe6fa50 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -162,5 +162,6 @@ extern "C" double k_generation(int i) { return simulation::k_generation.at(i - 1 extern "C" int k_generation_size() { return simulation::k_generation.size(); } extern "C" void k_generation_clear() { simulation::k_generation.clear(); } extern "C" void k_generation_reserve(int i) { simulation::k_generation.reserve(i); } +extern "C" int64_t work_index(int rank) { return simulation::work_index[rank]; } } // namespace openmc diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 65a36654f..1844c8da6 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -32,7 +32,6 @@ module simulation_header logical(C_BOOL), bind(C) :: satisfy_triggers ! whether triggers are satisfied integer(C_INT64_T), bind(C) :: work ! number of particles per processor - integer(C_INT64_T), allocatable :: work_index(:) ! starting index in source bank for each process integer(C_INT64_T), bind(C) :: current_work ! index in source bank of current history simulated ! ============================================================================ @@ -89,6 +88,12 @@ module simulation_header import C_INT integer(C_INT), value :: i end subroutine + + function work_index(rank) result(i) bind(C) + import C_INT, C_INT64_T + integer(C_INT), value :: rank + integer(C_INT64_T) :: i + end function end interface contains @@ -99,8 +104,6 @@ contains subroutine free_memory_simulation() - if (allocated(work_index)) deallocate(work_index) - call k_generation_clear() call entropy_clear() end subroutine free_memory_simulation diff --git a/src/source.cpp b/src/source.cpp index 1290c5f0d..3427579ee 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -264,7 +264,7 @@ void initialize_source() } // Read in the source bank - read_source_bank(file_id, simulation::work_index.data(), source_bank); + read_source_bank(file_id, source_bank); // Close file file_close(file_id); @@ -287,7 +287,7 @@ void initialize_source() write_message("Writing out initial source...", 5); std::string filename = settings::path_output + "initial_source.h5"; hid_t file_id = file_open(filename, 'w', true); - write_source_bank(file_id, simulation::work_index.data(), source_bank); + write_source_bank(file_id, source_bank); file_close(file_id); } } diff --git a/src/state_point.F90 b/src/state_point.F90 index f861f871a..52b01557d 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -36,17 +36,15 @@ module state_point implicit none interface - subroutine write_source_bank(group_id, work_index, bank_) bind(C) + subroutine write_source_bank(group_id, bank_) bind(C) import HID_T, C_INT64_T, Bank integer(HID_T), value :: group_id - integer(C_INT64_T), intent(in) :: work_index(*) type(Bank), intent(in) :: bank_(*) end subroutine write_source_bank - subroutine read_source_bank(group_id, work_index, bank_) bind(C) + subroutine read_source_bank(group_id, bank_) bind(C) import HID_T, C_INT64_T, Bank integer(HID_T), value :: group_id - integer(C_INT64_T), intent(in) :: work_index(*) type(Bank), intent(out) :: bank_(*) end subroutine read_source_bank end interface @@ -440,7 +438,7 @@ contains if (master .or. parallel) then file_id = file_open(filename_, 'a', parallel=.true.) end if - call write_source_bank(file_id, work_index, source_bank) + call write_source_bank(file_id, source_bank) if (master .or. parallel) call file_close(file_id) end if end function openmc_statepoint_write @@ -478,7 +476,7 @@ contains file_id = file_open(filename_, 'w', parallel=.true.) call write_attribute(file_id, "filetype", 'source') end if - call write_source_bank(file_id, work_index, source_bank) + call write_source_bank(file_id, source_bank) if (master .or. parallel) call file_close(file_id) end subroutine write_source_point @@ -684,7 +682,7 @@ contains end if ! Write out source - call read_source_bank(file_id, work_index, source_bank) + call read_source_bank(file_id, source_bank) end if diff --git a/src/state_point.cpp b/src/state_point.cpp index 4bd91e121..2383a5ed1 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -38,7 +38,7 @@ hid_t h5banktype() { void -write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) +write_source_bank(hid_t group_id, Bank* source_bank) { hid_t banktype = h5banktype(); @@ -54,7 +54,7 @@ write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) hid_t memspace = H5Screate_simple(1, count, nullptr); // Select hyperslab for this dataspace - hsize_t start[] {static_cast(work_index[openmc::mpi::rank])}; + hsize_t start[] {static_cast(simulation::work_index[mpi::rank])}; H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr); // Set up the property list for parallel writing @@ -84,21 +84,22 @@ write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) std::vector temp_source {source_bank, source_bank + simulation::work}; #endif - for (int i = 0; i < openmc::mpi::n_procs; ++i) { + for (int i = 0; i < mpi::n_procs; ++i) { // Create memory space - hsize_t count[] {static_cast(work_index[i+1] - work_index[i])}; + hsize_t count[] {static_cast(simulation::work_index[i+1] - + simulation::work_index[i])}; hid_t memspace = H5Screate_simple(1, count, nullptr); #ifdef OPENMC_MPI // Receive source sites from other processes if (i > 0) - MPI_Recv(source_bank, count[0], openmc::mpi::bank, i, i, - openmc::mpi::intracomm, MPI_STATUS_IGNORE); + MPI_Recv(source_bank, count[0], mpi::bank, i, i, + mpi::intracomm, MPI_STATUS_IGNORE); #endif // Select hyperslab for this dataspace dspace = H5Dget_space(dset); - hsize_t start[] {static_cast(work_index[i])}; + hsize_t start[] {static_cast(simulation::work_index[i])}; H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr); // Write data to hyperslab @@ -117,8 +118,8 @@ write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) #endif } else { #ifdef OPENMC_MPI - MPI_Send(source_bank, simulation::work, openmc::mpi::bank, 0, openmc::mpi::rank, - openmc::mpi::intracomm); + MPI_Send(source_bank, simulation::work, mpi::bank, 0, mpi::rank, + mpi::intracomm); #endif } #endif @@ -127,7 +128,7 @@ write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) } -void read_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) +void read_source_bank(hid_t group_id, Bank* source_bank) { hid_t banktype = h5banktype(); @@ -142,13 +143,13 @@ void read_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) hid_t dspace = H5Dget_space(dset); hsize_t dims_all[1]; H5Sget_simple_extent_dims(dspace, dims_all, nullptr); - if (work_index[openmc::mpi::n_procs] > dims_all[0]) { + if (simulation::work_index[mpi::n_procs] > dims_all[0]) { fatal_error("Number of source sites in source file is less " "than number of source particles per generation."); } // Select hyperslab for each process - hsize_t start[] {static_cast(work_index[openmc::mpi::rank])}; + hsize_t start[] {static_cast(simulation::work_index[mpi::rank])}; H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, dims, nullptr); #ifdef PHDF5 From 0e6e082bcbca28fde25776d8f7f4df3b18383b8c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 15 Oct 2018 10:24:30 -0500 Subject: [PATCH 19/23] Perform broadcasts for CMFD from C++ --- src/cmfd_execute.F90 | 21 ++++++++++++--------- src/cmfd_execute.cpp | 10 ++++++++++ src/cmfd_header.F90 | 4 ++-- 3 files changed, 24 insertions(+), 11 deletions(-) diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 18f5ebec1..bcd7e9f5c 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -13,6 +13,16 @@ module cmfd_execute private public :: execute_cmfd, cmfd_init_batch, cmfd_tally_init +#ifdef OPENMC_MPI + interface + subroutine cmfd_broadcast(n, buffer) bind(C) + import C_DOUBLE, C_INT + integer(C_INT), value :: n + real(C_DOUBLE), intent(out) :: buffer + end subroutine + end interface +#endif + contains !============================================================================== @@ -105,9 +115,6 @@ contains real(8) :: hxyz(3) ! cell dimensions of current ijk cell real(8) :: vol ! volume of cell real(8),allocatable :: source(:,:,:,:) ! tmp source array for entropy -#ifdef OPENMC_MPI - integer :: mpi_err ! MPI error code -#endif ! Get maximum of spatial and group indices nx = cmfd % indices(1) @@ -199,7 +206,7 @@ contains #ifdef OPENMC_MPI ! Broadcast full source to all procs - call MPI_BCAST(cmfd % cmfd_src, n, MPI_REAL8, 0, mpi_intracomm, mpi_err) + call cmfd_broadcast(n, cmfd % cmfd_src(1,1,1,1)) #endif end subroutine calc_fission_source @@ -232,9 +239,6 @@ contains real(8) :: norm ! normalization factor logical(C_BOOL) :: outside ! any source sites outside mesh logical :: in_mesh ! source site is inside mesh -#ifdef OPENMC_MPI - integer :: mpi_err -#endif interface subroutine cmfd_populate_sourcecounts(ng, energies, source_counts, outside) bind(C) @@ -300,8 +304,7 @@ contains ! Broadcast weight factors to all procs #ifdef OPENMC_MPI - call MPI_BCAST(cmfd % weightfactors, ng*nx*ny*nz, MPI_REAL8, 0, & - mpi_intracomm, mpi_err) + call cmfd_broadcast(ng*nx*ny*nz, cmfd % weightfactors(1,1,1,1)) #endif end if diff --git a/src/cmfd_execute.cpp b/src/cmfd_execute.cpp index 00d162069..aa9907a5b 100644 --- a/src/cmfd_execute.cpp +++ b/src/cmfd_execute.cpp @@ -7,6 +7,7 @@ #include "openmc/capi.h" #include "openmc/mesh.h" +#include "openmc/message_passing.h" #include "openmc/simulation.h" namespace openmc { @@ -30,4 +31,13 @@ cmfd_populate_sourcecounts(int n_energy, const double* energies, std::copy(counts.begin(), counts.end(), source_counts); } +#ifdef OPENMC_MPI +extern "C" void +cmfd_broadcast(int n, double* buffer) +{ + MPI_Bcast(buffer, n, MPI_DOUBLE, 0, mpi::intracomm); +} +#endif + + } // namespace openmc diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index e5e5cc423..036e193b0 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -51,8 +51,8 @@ module cmfd_header real(8), allocatable :: hxyz(:,:,:,:) ! Source distributions - real(8), allocatable :: cmfd_src(:,:,:,:) - real(8), allocatable :: openmc_src(:,:,:,:) + real(C_DOUBLE), allocatable :: cmfd_src(:,:,:,:) + real(C_DOUBLE), allocatable :: openmc_src(:,:,:,:) ! Source sites in each mesh box real(C_DOUBLE), allocatable :: sourcecounts(:,:) From 18d6e25e20cc227e647905932067216a50a4d205 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 15 Oct 2018 10:41:03 -0500 Subject: [PATCH 20/23] Remove last MPI calls on the Fortran side --- src/message_passing.cpp | 19 +++++++++++++++++++ src/volume_calc.F90 | 35 ++++++++++++++++++++++++++--------- 2 files changed, 45 insertions(+), 9 deletions(-) diff --git a/src/message_passing.cpp b/src/message_passing.cpp index dc287e3b3..8d8779cca 100644 --- a/src/message_passing.cpp +++ b/src/message_passing.cpp @@ -13,4 +13,23 @@ MPI_Datatype bank; #endif } // namespace mpi + +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +#ifdef OPENMC_MPI +extern "C" void +send_int(void* buffer, int count, int dest, int tag) +{ + MPI_Send(buffer, count, MPI_INTEGER, dest, tag, mpi::intracomm); +} + +extern "C" void +recv_int(void* buffer, int count, int source, int tag) +{ + MPI_Recv(buffer, count, MPI_INTEGER, source, tag, mpi::intracomm, MPI_STATUS_IGNORE); +} +#endif + } // namespace openmc diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index f7ded7361..5379a3b2e 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -138,16 +138,35 @@ contains integer :: min_samples ! minimum number of samples per process integer :: remainder ! leftover samples from uneven divide #ifdef OPENMC_MPI - integer :: mpi_err ! MPI error code integer :: m ! index over materials - integer :: n ! number of materials - integer, allocatable :: data(:) ! array used to send number of hits + integer(C_INT) :: n ! number of materials + integer(C_INT), allocatable :: data(:) ! array used to send number of hits #endif real(8) :: f ! fraction of hits real(8) :: var_f ! variance of fraction of hits real(8) :: volume_sample ! total volume of sampled region real(8) :: atoms(2, size(nuclides)) +#ifdef OPENMC_MPI + interface + subroutine send_int(buffer, count, dest, tag) bind(C) + import C_INT + integer(C_INT), intent(in) :: buffer + integer(C_INT), value :: count + integer(C_INT), value :: dest + integer(C_INT), value :: tag + end subroutine + + subroutine recv_int(buffer, count, source, tag) bind(C) + import C_INT + integer(C_INT), intent(out) :: buffer + integer(C_INT), value :: count + integer(C_INT), value :: source + integer(C_INT), value :: tag + end subroutine + end interface +#endif + ! Divide work over MPI processes min_samples = this % samples / n_procs remainder = mod(this % samples, n_procs) @@ -285,12 +304,10 @@ contains if (master) then #ifdef OPENMC_MPI do j = 1, n_procs - 1 - call MPI_RECV(n, 1, MPI_INTEGER, j, 0, mpi_intracomm, & - MPI_STATUS_IGNORE, mpi_err) + call recv_int(n, 1, j, 0) allocate(data(2*n)) - call MPI_RECV(data, 2*n, MPI_INTEGER, j, 1, mpi_intracomm, & - MPI_STATUS_IGNORE, mpi_err) + call recv_int(data(1), 2*n, j, 1) do k = 0, n - 1 do m = 1, master_indices(i_domain) % size() if (data(2*k + 1) == master_indices(i_domain) % data(m)) then @@ -353,8 +370,8 @@ contains data(2*k + 2) = master_hits(i_domain) % data(k + 1) end do - call MPI_SEND(n, 1, MPI_INTEGER, 0, 0, mpi_intracomm, mpi_err) - call MPI_SEND(data, 2*n, MPI_INTEGER, 0, 1, mpi_intracomm, mpi_err) + call send_int(n, 1, 0, 0) + call send_int(data(1), 2*n, 0, 1) deallocate(data) #endif end if From bb056d597e259a2c2d61a3470619ab5ff283a6b7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 15 Oct 2018 10:52:27 -0500 Subject: [PATCH 21/23] Remove MPI intracomm from Fortran side --- CMakeLists.txt | 4 ---- src/hdf5_interface.F90 | 3 --- src/initialize.F90 | 30 +----------------------------- src/message_passing.F90 | 20 +++----------------- 4 files changed, 4 insertions(+), 53 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 84f21427a..c16d72662 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,7 +18,6 @@ option(profile "Compile with profiling flags" OFF) 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(mpif08 "Use Fortran 2008 MPI interface" OFF) option(dagmc "Enable support for DAGMC (CAD) geometry" OFF) # Maximum number of nested coordinates levels @@ -474,9 +473,6 @@ if (HDF5_IS_PARALLEL) endif() if (MPI_ENABLED) target_compile_definitions(libopenmc PUBLIC -DOPENMC_MPI) - if (mpif08) - target_compile_definitions(libopenmc PRIVATE -DOPENMC_MPIF08) - endif() endif() # Set git SHA1 hash as a compile definition diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 8872a84f6..7856abb50 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -13,9 +13,6 @@ module hdf5_interface use, intrinsic :: ISO_C_BINDING use error, only: fatal_error -#ifdef PHDF5 - use message_passing, only: mpi_intracomm, MPI_INFO_NULL -#endif use string, only: to_c_string, to_f_string implicit none diff --git a/src/initialize.F90 b/src/initialize.F90 index 697fc1235..3d805c97d 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -49,45 +49,17 @@ contains ! setting up timers, etc. !=============================================================================== - function openmc_init_f(intracomm) result(err) bind(C) - integer, intent(in), optional :: intracomm ! MPI intracommunicator + function openmc_init_f() result(err) bind(C) integer(C_INT) :: err #ifdef _OPENMP character(MAX_WORD_LEN) :: envvar #endif - ! Copy the communicator to a new variable. This is done to avoid changing - ! the signature of this subroutine. If MPI is being used but no communicator - ! was passed, assume MPI_COMM_WORLD. -#ifdef OPENMC_MPI -#ifdef OPENMC_MPIF08 - type(MPI_Comm), intent(in) :: comm ! MPI intracommunicator - if (present(intracomm)) then - comm % MPI_VAL = intracomm - else - comm = MPI_COMM_WORLD - end if -#else - integer :: comm - if (present(intracomm)) then - comm = intracomm - else - comm = MPI_COMM_WORLD - end if -#endif -#endif - ! Start total and initialization timer call time_total%start() call time_initialize%start() -#ifdef OPENMC_MPI - ! Indicate that MPI is turned on - mpi_enabled = .true. - mpi_intracomm = intracomm -#endif - #ifdef _OPENMP ! Change schedule of main parallel-do loop if OMP_SCHEDULE is set call get_environment_variable("OMP_SCHEDULE", envvar) diff --git a/src/message_passing.F90 b/src/message_passing.F90 index 05219f0f8..37d203447 100644 --- a/src/message_passing.F90 +++ b/src/message_passing.F90 @@ -2,26 +2,12 @@ module message_passing use, intrinsic :: ISO_C_BINDING -#ifdef OPENMC_MPI -#ifdef OPENMC_MPIF08 - use mpi_f08 -#else - use mpi -#endif -#endif - - ! The defaults set here for the number of processors, rank, and master and - ! mpi_enabled flag are for when MPI is not being used at all, i.e. a serial - ! run. In this case, these variables are still used at times. + ! The defaults set here for the number of processors, rank, and master and are + ! for when MPI is not being used at all, i.e. a serial run. In this case, these + ! variables are still used at times. integer(C_INT), bind(C, name='openmc_n_procs') :: n_procs = 1 ! number of processes integer(C_INT), bind(C, name='openmc_rank') :: rank = 0 ! rank of process logical(C_BOOL), bind(C, name='openmc_master') :: master = .true. ! master process? - logical :: mpi_enabled = .false. ! is MPI in use and initialized? -#ifdef OPENMC_MPIF08 - type(MPI_Comm) :: mpi_intracomm ! MPI intra-communicator -#else - integer :: mpi_intracomm ! MPI intra-communicator -#endif end module message_passing From 77951e1b7012a6521bc531f5a09536c0ed058283 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 16 Oct 2018 14:43:19 -0500 Subject: [PATCH 22/23] Fix signature on openmc_init_f --- include/openmc/capi.h | 1 - src/initialize.cpp | 8 ++------ 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index f7d69290c..05512b379 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -52,7 +52,6 @@ extern "C" { int openmc_global_tallies(double** ptr); int openmc_hard_reset(); int openmc_init(int argc, char* argv[], const void* intracomm); - int openmc_init_f(const int* intracomm); int openmc_legendre_filter_get_order(int32_t index, int* order); int openmc_legendre_filter_set_order(int32_t index, int order); int openmc_load_nuclide(const char name[]); diff --git a/src/initialize.cpp b/src/initialize.cpp index 31fa8fb25..26d21bf09 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -19,6 +19,7 @@ #include "openmc/string_utils.h" // data/functions from Fortran side +extern "C" void openmc_init_f(); extern "C" void print_usage(); extern "C" void print_version(); @@ -47,12 +48,7 @@ int openmc_init(int argc, char* argv[], const void* intracomm) if (err) return err; // Continue with rest of initialization -#ifdef OPENMC_MPI - MPI_Fint fcomm = MPI_Comm_c2f(comm); - openmc_init_f(&fcomm); -#else - openmc_init_f(nullptr); -#endif + openmc_init_f(); return 0; } From 81858c119c8b6d3b6a77205369cea5cb9157c6bd Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 17 Oct 2018 13:48:30 -0500 Subject: [PATCH 23/23] Revert include ordering in state_point.cpp. Add comment about adaptor_type --- include/openmc/tallies/tally.h | 2 ++ src/state_point.cpp | 14 +++++++------- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 96b081615..c323d3bac 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -22,6 +22,8 @@ extern "C" double global_tally_absorption; // Non-member functions //============================================================================== +// Alias for the type returned by xt::adapt(...). N is the dimension of the +// multidimensional array template using adaptor_type = xt::xtensor_adaptor, N>; diff --git a/src/state_point.cpp b/src/state_point.cpp index 2383a5ed1..63c8264a3 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -1,5 +1,12 @@ #include "openmc/state_point.h" +#include +#include +#include + +#include "xtensor/xbuilder.hpp" // for empty_like +#include "xtensor/xview.hpp" + #include "openmc/capi.h" #include "openmc/constants.h" #include "openmc/error.h" @@ -9,13 +16,6 @@ #include "openmc/simulation.h" #include "openmc/tallies/tally.h" -#include "xtensor/xbuilder.hpp" // for empty_like -#include "xtensor/xview.hpp" - -#include -#include -#include - namespace openmc {