diff --git a/CMakeLists.txt b/CMakeLists.txt index 32beb41ecd..24018eb0f1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -159,6 +159,13 @@ if(${HDF5_VERSION} VERSION_GREATER_EQUAL 1.12.0) list(APPEND cxxflags -DH5Oget_info_by_idx_vers=1 -DH5O_info_t_vers=1) endif() +#=============================================================================== +# MCPL +#=============================================================================== +if (OPENMC_USE_MCPL) + find_package(MCPL REQUIRED) +endif() + #=============================================================================== # Set compile/link flags based on which compiler is being used #=============================================================================== @@ -186,10 +193,6 @@ if(OPENMC_ENABLE_COVERAGE) list(APPEND ldflags --coverage) endif() -if(OPENMC_USE_MCPL) - list(APPEND ldflags -lmcpl) -endif() - # Show flags being used message(STATUS "OpenMC C++ flags: ${cxxflags}") message(STATUS "OpenMC Linker flags: ${ldflags}") @@ -444,10 +447,6 @@ endif() if (OPENMC_USE_MPI) target_compile_definitions(libopenmc PUBLIC -DOPENMC_MPI) endif() -if(OPENMC_USE_MCPL) - target_compile_definitions(libopenmc PUBLIC -DOPENMC_MCPL) -endif() - # Set git SHA1 hash as a compile definition if(GIT_FOUND) @@ -491,6 +490,11 @@ if (OPENMC_USE_MPI) target_link_libraries(libopenmc MPI::MPI_CXX) endif() +if (OPENMC_USE_MCPL) + target_compile_definitions(libopenmc PUBLIC OPENMC_MCPL) + target_link_libraries(libopenmc MCPL::mcpl) +endif() + #=============================================================================== # openmc executable #=============================================================================== diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 1e061b235b..66825867d8 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -48,6 +48,7 @@ extern bool source_latest; //!< write latest source at each batch? extern bool source_separate; //!< write source to separate file? extern bool source_write; //!< write source in HDF5 files? extern bool surf_source_write; //!< write surface source file? +extern bool surf_mcpl_write; //!< write surface mcpl file? extern bool surf_source_read; //!< read surface source file? extern bool survival_biasing; //!< use survival biasing? extern bool temperature_multipole; //!< use multipole data? diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index 5ada2bf88d..2d92c935a6 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -9,6 +9,10 @@ #include "openmc/particle.h" #include "openmc/vector.h" +#ifdef OPENMC_MCPL +#include +#endif + namespace openmc { void load_state_point(); @@ -21,5 +25,10 @@ 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/openmc/settings.py b/openmc/settings.py index ad5a9b28ae..0937a3e480 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -171,6 +171,7 @@ class Settings: banked (int) :max_particles: Maximum number of particles to be banked on surfaces per process (int) + :mcpl: Output in the form of an MCPL-file (bool) survival_biasing : bool Indicate whether survival biasing is to be used tabular_legendre : dict @@ -652,7 +653,7 @@ class Settings: cv.check_type('surface source writing options', surf_source_write, Mapping) for key, value in surf_source_write.items(): cv.check_value('surface source writing key', key, - ('surface_ids', 'max_particles')) + ('surface_ids', 'max_particles', 'mcpl')) if key == 'surface_ids': cv.check_type('surface ids for source banking', value, Iterable, Integral) @@ -664,6 +665,9 @@ class Settings: value, Integral) cv.check_greater_than('maximum particle banks on surfaces per process', value, 0) + elif key == 'mcpl': + cv.check_type('write to an MCPL-format file', value, bool) + self._surf_source_write = surf_source_write @confidence_intervals.setter @@ -1032,6 +1036,9 @@ class Settings: if 'max_particles' in self._surf_source_write: subelement = ET.SubElement(element, "max_particles") subelement.text = str(self._surf_source_write['max_particles']) + if 'mcpl' in self._surf_source_write: + subelement = ET.SubElement(element, "mcpl") + subelement.text = str(self._surf_source_write['mcpl']).lower() def _create_confidence_intervals(self, root): if self._confidence_intervals is not None: diff --git a/openmc/source.py b/openmc/source.py index 48e14f6f71..fdb2d5bc25 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -223,7 +223,10 @@ class Source: if self.particle != 'neutron': element.set("particle", self.particle) if self.file is not None: - element.set("file", self.file) + if (self.file.endswith('.mcpl')): + element.set("mcpl", self.file) + else: + element.set("file", self.file) if self.library is not None: element.set("library", self.library) if self.parameters is not None: diff --git a/src/settings.cpp b/src/settings.cpp index 11dd96e025..d1762596c6 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -61,6 +61,7 @@ bool source_latest {false}; bool source_separate {false}; bool source_write {true}; bool surf_source_write {false}; +bool surf_mcpl_write {false}; bool surf_source_read {false}; bool survival_biasing {false}; bool temperature_multipole {false}; @@ -682,9 +683,14 @@ 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; + } +#endif } - // If source is not seperate and is to be written out in the statepoint file, + // If source is not separate and is to be written out in the statepoint file, // make sure that the sourcepoint batch numbers are contained in the // statepoint list if (!source_separate) { diff --git a/src/simulation.cpp b/src/simulation.cpp index d7292cb9d9..d90c33ffd6 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -407,6 +407,13 @@ 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/source.cpp b/src/source.cpp index e43f03929b..1c46e4939c 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -343,17 +343,28 @@ FileSource::FileSource(mcpl_file_t mcpl_file) const mcpl_particle_t *mcpl_particle; //extract particle from mcpl-file mcpl_particle=mcpl_read(mcpl_file); - // check if it is a neutron or a photon. otherwise skip - while ( mcpl_particle->pdgcode!=2112 && mcpl_particle->pdgcode!=22 ) { + // check if it is a neutron, photon, electron, or positron. Otherwise skip. + int pdg=mcpl_particle->pdgcode; + while ( pdg!=2112 && pdg!=22 && pdg!=11 && pdg!=-11) { mcpl_particle=mcpl_read(mcpl_file); - //should check for file exhaustion This could happen if particles are other than - //neutrons or photons + pdg=mcpl_particle->pdgcode; + //should check for file exhaustion. This could happen if particles are other than + //neutrons, photons, electrons, or positrons. } - if(mcpl_particle->pdgcode==2112) { - site_.particle=ParticleType::neutron; - } else if (mcpl_particle->pdgcode==22) { - site_.particle=ParticleType::photon; + switch(pdg){ + case 2112: + site_.particle=ParticleType::neutron; + break; + case 22: + site_.particle=ParticleType::photon; + break; + case 11: + site_.particle=ParticleType::electron; + break; + case -11: + site_.particle=ParticleType::positron; + break; } //particle is good, convert to openmc-formalism diff --git a/src/state_point.cpp b/src/state_point.cpp index 470dd77182..983147c62e 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -28,6 +28,10 @@ #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) @@ -599,6 +603,46 @@ 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) { + // this must be rewritten: file_open is h5-specific + file_id = mcpl_create_outfile(filename_.c_str()); + //write_attribute(file_id, "filetype", "source"); + //write header stuff (oopy in xml-files as binary blobs for instance)) + 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) { + //change this - this is h5 specific + mcpl_closeandgzip_outfile(file_id); + } + +} +#endif + void write_source_bank(hid_t group_id, bool surf_source_bank) { hid_t banktype = h5banktype(); @@ -715,6 +759,105 @@ 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) {