diff --git a/CMakeLists.txt b/CMakeLists.txt index af3a5b4af..466816b30 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -164,6 +164,7 @@ endif() #=============================================================================== if (OPENMC_USE_MCPL) find_package(MCPL REQUIRED) + message(STATUS "Found MCPL: ${MCPL_DIR} (found version \"${MCPL_VERSION}\")") endif() #=============================================================================== diff --git a/include/openmc/mcpl_interface.h b/include/openmc/mcpl_interface.h index 89dab16bc..c931f1774 100644 --- a/include/openmc/mcpl_interface.h +++ b/include/openmc/mcpl_interface.h @@ -6,20 +6,31 @@ #include #include -#ifdef OPENMC_MCPL -#include -#endif - namespace openmc { +//============================================================================== +// Constants +//============================================================================== + extern "C" const bool MCPL_ENABLED; +//============================================================================== +// Functions +//============================================================================== + //! Get a vector of source sites from an MCPL file // //! \param[in] path Path to MCPL file //! \return Vector of source sites vector mcpl_source_sites(std::string path); +//! Write an MCPL source file +// +//! \param[in] filename Path to MCPL file +//! \param[in] surf_source_bank Whether to use the surface source bank +void write_mcpl_source_point( + const char* filename, bool surf_source_bank = false); + } // namespace openmc #endif // OPENMC_MCPL_INTERFACE_H diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index 730127554..5ada2bf88 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -9,10 +9,6 @@ #include "openmc/particle.h" #include "openmc/vector.h" -#ifdef OPENMC_MCPL -#include -#endif - namespace openmc { void load_state_point(); @@ -25,11 +21,5 @@ void write_tally_results_nr(hid_t file_id); void restart_set_keff(); void write_unstructured_mesh_results(); -#ifdef OPENMC_MCPL -void write_mcpl_source_point( - const char* filename, bool surf_source_bank = false); -void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank); -#endif - } // namespace openmc #endif // OPENMC_STATE_POINT_H diff --git a/src/mcpl_interface.cpp b/src/mcpl_interface.cpp index d55cc5af7..e9bc64830 100644 --- a/src/mcpl_interface.cpp +++ b/src/mcpl_interface.cpp @@ -1,6 +1,13 @@ #include "openmc/mcpl_interface.h" +#include "openmc/bank.h" #include "openmc/error.h" +#include "openmc/message_passing.h" +#include "openmc/settings.h" +#include "openmc/simulation.h" +#include "openmc/state_point.h" + +#include #ifdef OPENMC_MCPL #include @@ -8,14 +15,22 @@ namespace openmc { +//============================================================================== +// Constants +//============================================================================== + #ifdef OPENMC_MCPL const bool MCPL_ENABLED = true; #else const bool MCPL_ENABLED = false; #endif +//============================================================================== +// Functions +//============================================================================== + #ifdef OPENMC_MCPL -SourceSite mcpl_particle_to_site(const mcpl_particle* particle) +SourceSite mcpl_particle_to_site(const mcpl_particle_t* particle) { SourceSite site; @@ -35,40 +50,42 @@ SourceSite mcpl_particle_to_site(const mcpl_particle* particle) } // Copy position and direction - site.r.x = mcpl_particle->position[0]; - site.r.y = mcpl_particle->position[1]; - site.r.z = mcpl_particle->position[2]; - site.u.x = mcpl_particle->direction[0]; - site.u.y = mcpl_particle->direction[1]; - site.u.z = mcpl_particle->direction[2]; + site.r.x = particle->position[0]; + site.r.y = particle->position[1]; + site.r.z = particle->position[2]; + site.u.x = particle->direction[0]; + site.u.y = particle->direction[1]; + site.u.z = particle->direction[2]; // mcpl stores kinetic energy in MeV - site.E = mcpl_particle->ekin * 1e6; + site.E = particle->ekin * 1e6; // mcpl stores time in ms - site.time = mcpl_particle->time * 1e-3; - site.wgt = mcpl_particle->weight; + site.time = particle->time * 1e-3; + site.wgt = particle->weight; return site; } #endif +//============================================================================== + vector mcpl_source_sites(std::string path) { vector sites; #ifdef OPENMC_MCPL + // Open MCPL file and determine number of particles + auto mcpl_file = mcpl_open_file(path.c_str()); size_t n_sites = mcpl_hdr_nparticles(mcpl_file); for (int i = 0; i < n_sites; i++) { - SourceSite site; - // Extract particle from mcpl-file, checking if it is a neutron, photon, // electron, or positron. Otherwise skip. const mcpl_particle_t* particle; int pdg = 0; while (pdg != 2112 && pdg != 22 && pdg != 11 && pdg != -11) { particle = mcpl_read(mcpl_file); - pdg = mcpl_particle->pdgcode; + pdg = particle->pdgcode; } // Convert to source site and add to vector @@ -76,7 +93,7 @@ vector mcpl_source_sites(std::string path) } // Check that some sites were read - if (sites_.empty()) { + if (sites.empty()) { fatal_error("MCPL file contained no neutron, photon, electron, or positron " "source particles."); } @@ -90,4 +107,144 @@ vector mcpl_source_sites(std::string path) return sites; } +//============================================================================== + +#ifdef OPENMC_MCPL +void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank) +{ + int64_t dims_size = settings::n_particles; + int64_t count_size = simulation::work_per_rank; + + // Set vectors for source bank and starting bank index of each process + vector* bank_index = &simulation::work_index; + vector* source_bank = &simulation::source_bank; + vector surf_source_index_vector; + vector surf_source_bank_vector; + + if (surf_source_bank) { + surf_source_index_vector = calculate_surf_source_size(); + dims_size = surf_source_index_vector[mpi::n_procs]; + count_size = simulation::surf_source_bank.size(); + + bank_index = &surf_source_index_vector; + + // Copy data in a SharedArray into a vector. + surf_source_bank_vector.resize(count_size); + surf_source_bank_vector.assign(simulation::surf_source_bank.data(), + simulation::surf_source_bank.data() + count_size); + source_bank = &surf_source_bank_vector; + } + + if (mpi::master) { + // Particles are writeen to disk from the master node only + + // Save source bank sites since the array is overwritten below +#ifdef OPENMC_MPI + vector temp_source {source_bank->begin(), source_bank->end()}; +#endif + + // loop over the other nodes and receive data - then write those. + for (int i = 0; i < mpi::n_procs; ++i) { + // number of particles for node node i + size_t count[] { + static_cast((*bank_index)[i + 1] - (*bank_index)[i])}; + +#ifdef OPENMC_MPI + if (i > 0) + MPI_Recv(source_bank->data(), count[0], mpi::source_site, i, i, + mpi::intracomm, MPI_STATUS_IGNORE); +#endif + // now write the source_bank data again. + for (vector::iterator _site = source_bank->begin(); + _site != source_bank->end(); _site++) { + // particle is now at the iterator + // write it to the mcpl-file + mcpl_particle_t p; + p.position[0] = _site->r.x; + p.position[1] = _site->r.y; + p.position[2] = _site->r.z; + + // mcpl requires that the direction vector is unit length + // which is also the case in openmc + p.direction[0] = _site->u.x; + p.direction[1] = _site->u.y; + p.direction[2] = _site->u.z; + + // mcpl stores kinetic energy in MeV + p.ekin = _site->E * 1e-6; + + p.time = _site->time * 1e3; + + p.weight = _site->wgt; + + switch (_site->particle) { + case ParticleType::neutron: + p.pdgcode = 2112; + break; + case ParticleType::photon: + p.pdgcode = 22; + break; + case ParticleType::electron: + p.pdgcode = 11; + break; + case ParticleType::positron: + p.pdgcode = -11; + break; + } + + mcpl_add_particle(file_id, &p); + } + } +#ifdef OPENMC_MPI + // Restore state of source bank + std::copy(temp_source.begin(), temp_source.end(), source_bank->begin()); +#endif + } else { +#ifdef OPENMC_MPI + MPI_Send(source_bank->data(), count_size, mpi::source_site, 0, mpi::rank, + mpi::intracomm); +#endif + } +} +#endif + +//============================================================================== + +void write_mcpl_source_point(const char* filename, bool surf_source_bank) +{ + std::string filename_; + if (filename) { + filename_ = filename; + } else { + // Determine width for zero padding + int w = std::to_string(settings::n_max_batches).size(); + + filename_ = fmt::format("{0}source.{1:0{2}}.mcpl", settings::path_output, + simulation::current_batch, w); + } + +#ifdef OPENMC_MCPL + mcpl_outfile_t file_id; + + std::string line; + if (mpi::master) { + file_id = mcpl_create_outfile(filename_.c_str()); + if (VERSION_DEV) { + line = fmt::format("OpenMC {0}.{1}.{2}-development", VERSION_MAJOR, + VERSION_MINOR, VERSION_RELEASE); + } else { + line = fmt::format( + "OpenMC {0}.{1}.{2}", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE); + } + mcpl_hdr_set_srcname(file_id, line.c_str()); + } + + write_mcpl_source_bank(file_id, surf_source_bank); + + if (mpi::master) { + mcpl_close_outfile(file_id); + } +#endif +} + } // namespace openmc diff --git a/src/settings.cpp b/src/settings.cpp index bc31e2c4d..cb39d0c8b 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -433,10 +433,8 @@ void read_settings_xml() if (check_for_node(node, "file")) { auto path = get_node_value(node, "file", false, true); if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) { -#ifdef OPENMC_MCPL auto sites = mcpl_source_sites(path); model::external_sources.push_back(make_unique(sites)); -#endif } else { model::external_sources.push_back(make_unique(path)); } @@ -659,6 +657,12 @@ void read_settings_xml() } if (check_for_node(node_sp, "mcpl")) { source_mcpl_write = get_node_value_bool(node_sp, "mcpl"); + + // Make sure MCPL support is enabled + if (source_mcpl_write && !MCPL_ENABLED) { + fatal_error( + "Your build of OpenMC does not support writing MCPL source files."); + } } if (check_for_node(node_sp, "overwrite_latest")) { source_latest = get_node_value_bool(node_sp, "overwrite_latest"); @@ -690,11 +694,15 @@ void read_settings_xml() max_surface_particles = std::stoll(get_node_value(node_ssw, "max_particles")); } -#ifdef OPENMC_MCPL if (check_for_node(node_ssw, "mcpl")) { - surf_mcpl_write = true; + surf_mcpl_write = get_node_value_bool(node_ssw, "mcpl"); + + // Make sure MCPL support is enabled + if (surf_mcpl_write && !MCPL_ENABLED) { + fatal_error("Your build of OpenMC does not support writing MCPL " + "surface source files."); + } } -#endif } // If source is not separate and is to be written out in the statepoint file, diff --git a/src/simulation.cpp b/src/simulation.cpp index 2c1d1c988..497b26124 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -8,6 +8,7 @@ #include "openmc/event.h" #include "openmc/geometry_aux.h" #include "openmc/material.h" +#include "openmc/mcpl_interface.h" #include "openmc/message_passing.h" #include "openmc/nuclide.h" #include "openmc/output.h" @@ -389,9 +390,7 @@ void finalize_batch() if (settings::run_mode == RunMode::EIGENVALUE) { // Write out a separate source point if it's been specified for this batch -#ifdef OPENMC_MCPL if (!settings::source_mcpl_write) { -#endif if (contains(settings::sourcepoint_batch, simulation::current_batch) && settings::source_write && settings::source_separate) { write_source_point(nullptr); @@ -402,7 +401,6 @@ void finalize_batch() auto filename = settings::path_output + "source.h5"; write_source_point(filename.c_str()); } -#ifdef OPENMC_MCPL } else { if (contains(settings::sourcepoint_batch, simulation::current_batch) && settings::source_mcpl_write && settings::source_separate) { @@ -415,7 +413,6 @@ void finalize_batch() write_mcpl_source_point(filename.c_str()); } } -#endif } // Write out surface source if requested. @@ -424,13 +421,11 @@ void finalize_batch() auto filename = settings::path_output + "surface_source.h5"; write_source_point(filename.c_str(), true); } -#ifdef OPENMC_MCPL if (settings::surf_mcpl_write && simulation::current_batch == settings::n_batches) { auto filename = settings::path_output + "surface_source.mcpl"; write_mcpl_source_point(filename.c_str(), true); } -#endif } void initialize_generation() diff --git a/src/state_point.cpp b/src/state_point.cpp index 9dd0f1374..470dd7718 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -28,10 +28,6 @@ #include "openmc/timer.h" #include "openmc/vector.h" -#ifdef OPENMC_MCPL -#include -#endif - namespace openmc { extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) @@ -603,43 +599,6 @@ void write_source_point(const char* filename, bool surf_source_bank) file_close(file_id); } -#ifdef OPENMC_MCPL -void write_mcpl_source_point(const char* filename, bool surf_source_bank) -{ - std::string filename_; - if (filename) { - filename_ = filename; - } else { - // Determine width for zero padding - int w = std::to_string(settings::n_max_batches).size(); - - filename_ = fmt::format("{0}source.{1:0{2}}.mcpl", settings::path_output, - simulation::current_batch, w); - } - - mcpl_outfile_t file_id; - - std::string line; - if (mpi::master) { - file_id = mcpl_create_outfile(filename_.c_str()); - if (VERSION_DEV) { - line = fmt::format("OpenMC {0}.{1}.{2}-development", VERSION_MAJOR, - VERSION_MINOR, VERSION_RELEASE); - } else { - line = fmt::format( - "OpenMC {0}.{1}.{2}", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE); - } - mcpl_hdr_set_srcname(file_id, line.c_str()); - } - - write_mcpl_source_bank(file_id, surf_source_bank); - - if (mpi::master) { - mcpl_close_outfile(file_id); - } -} -#endif - void write_source_bank(hid_t group_id, bool surf_source_bank) { hid_t banktype = h5banktype(); @@ -756,105 +715,6 @@ void write_source_bank(hid_t group_id, bool surf_source_bank) H5Tclose(banktype); } -#ifdef OPENMC_MCPL -void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank) -{ - int64_t dims_size = settings::n_particles; - int64_t count_size = simulation::work_per_rank; - - // Set vectors for source bank and starting bank index of each process - vector* bank_index = &simulation::work_index; - vector* source_bank = &simulation::source_bank; - vector surf_source_index_vector; - vector surf_source_bank_vector; - - if (surf_source_bank) { - surf_source_index_vector = calculate_surf_source_size(); - dims_size = surf_source_index_vector[mpi::n_procs]; - count_size = simulation::surf_source_bank.size(); - - bank_index = &surf_source_index_vector; - - // Copy data in a SharedArray into a vector. - surf_source_bank_vector.resize(count_size); - surf_source_bank_vector.assign(simulation::surf_source_bank.data(), - simulation::surf_source_bank.data() + count_size); - source_bank = &surf_source_bank_vector; - } - - if (mpi::master) { - // Particles are writeen to disk from the master node only - - // Save source bank sites since the array is overwritten below -#ifdef OPENMC_MPI - vector temp_source {source_bank->begin(), source_bank->end()}; -#endif - - // loop over the other nodes and receive data - then write those. - for (int i = 0; i < mpi::n_procs; ++i) { - // number of particles for node node i - size_t count[] { - static_cast((*bank_index)[i + 1] - (*bank_index)[i])}; - -#ifdef OPENMC_MPI - if (i > 0) - MPI_Recv(source_bank->data(), count[0], mpi::source_site, i, i, - mpi::intracomm, MPI_STATUS_IGNORE); -#endif - // now write the source_bank data again. - for (vector::iterator _site = source_bank->begin(); - _site != source_bank->end(); _site++) { - // particle is now at the iterator - // write it to the mcpl-file - mcpl_particle_t p; - p.position[0] = _site->r.x; - p.position[1] = _site->r.y; - p.position[2] = _site->r.z; - - // mcpl requires that the direction vector is unit length - // which is also the case in openmc - p.direction[0] = _site->u.x; - p.direction[1] = _site->u.y; - p.direction[2] = _site->u.z; - - // mcpl stores kinetic energy in MeV - p.ekin = _site->E * 1e-6; - - p.time = _site->time * 1e3; - - p.weight = _site->wgt; - - switch (_site->particle) { - case ParticleType::neutron: - p.pdgcode = 2112; - break; - case ParticleType::photon: - p.pdgcode = 22; - break; - case ParticleType::electron: - p.pdgcode = 11; - break; - case ParticleType::positron: - p.pdgcode = -11; - break; - } - - mcpl_add_particle(file_id, &p); - } - } -#ifdef OPENMC_MPI - // Restore state of source bank - std::copy(temp_source.begin(), temp_source.end(), source_bank->begin()); -#endif - } else { -#ifdef OPENMC_MPI - MPI_Send(source_bank->data(), count_size, mpi::source_site, 0, mpi::rank, - mpi::intracomm); -#endif - } -} -#endif - // Determine member names of a compound HDF5 datatype std::string dtype_member_names(hid_t dtype_id) {