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..5b755784e4 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -47,7 +47,9 @@ extern "C" bool run_CE; //!< run with continuous-energy data? 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 source_mcpl_write; //!< write source in mcpl 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..9098eb6205 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -60,7 +60,9 @@ bool run_CE {true}; bool source_latest {false}; bool source_separate {false}; bool source_write {true}; +bool source_mcpl_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}; @@ -652,6 +654,9 @@ void read_settings_xml() if (check_for_node(node_sp, "write")) { source_write = get_node_value_bool(node_sp, "write"); } + if (check_for_node(node_sp, "mcpl")) { + source_write = get_node_value_bool(node_sp, "mcpl"); + } if (check_for_node(node_sp, "overwrite_latest")) { source_latest = get_node_value_bool(node_sp, "overwrite_latest"); source_separate = source_latest; @@ -682,9 +687,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..d7cc2237d5 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -389,16 +389,33 @@ void finalize_batch() if (settings::run_mode == RunMode::EIGENVALUE) { // Write out a separate source point if it's been specified for this batch - if (contains(settings::sourcepoint_batch, simulation::current_batch) && - settings::source_write && settings::source_separate) { - write_source_point(nullptr); - } +#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); + } - // Write a continously-overwritten source point if requested. - if (settings::source_latest) { - auto filename = settings::path_output + "source.h5"; - write_source_point(filename.c_str()); + // Write a continously-overwritten source point if requested. + if (settings::source_latest) { + 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) { + write_mcpl_source_point(nullptr); + } + + // Write a continously-overwritten source point if requested. + if (settings::source_latest && settings::source_mcpl_write) { + auto filename = settings::path_output + "source.mcpl"; + write_mcpl_source_point(filename.c_str()); + } } +#endif } // Write out surface source if requested. @@ -407,6 +424,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..d0fe92bdf5 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_close_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) { diff --git a/tests/regression_tests/source_mcpl_file/results_true.dat b/tests/regression_tests/source_mcpl_file/results_true.dat index f2209006a1..19460623f5 100644 --- a/tests/regression_tests/source_mcpl_file/results_true.dat +++ b/tests/regression_tests/source_mcpl_file/results_true.dat @@ -1,2 +1,2 @@ k-combined: -3.017557E-01 3.398770E-03 +3.039964E-01 3.654869E-04 diff --git a/tests/regression_tests/source_mcpl_file/settings.xml b/tests/regression_tests/source_mcpl_file/settings.xml index 8b1510e0ce..47b010bffe 100644 --- a/tests/regression_tests/source_mcpl_file/settings.xml +++ b/tests/regression_tests/source_mcpl_file/settings.xml @@ -1,11 +1,15 @@ + + 10 5 1000 - source.10.mcpl + + -4 -4 -4 4 4 4 + diff --git a/tests/regression_tests/source_mcpl_file/test.py b/tests/regression_tests/source_mcpl_file/test.py index f55a769047..a005fc2d83 100644 --- a/tests/regression_tests/source_mcpl_file/test.py +++ b/tests/regression_tests/source_mcpl_file/test.py @@ -1,65 +1,101 @@ -import os -import subprocess +#!/usr/bin/env python + import glob +import os -from tests.testing_harness import TestHarness +from tests.testing_harness import * -class SourceMCPLFileTestHarness(TestHarness): - def execute_test(self): - """Run OpenMC with the appropriate arguments and check the outputs.""" - try: - self._create_input() - self._test_input_created() - self._run_openmc() - self._test_output_created() - results = self._get_results() - self._write_results(results) - self._compare_results() - finally: - self._cleanup() +settings1=""" + + + + + 10 + 5 + 1000 + + + + -4 -4 -4 4 4 4 + + + +""" - def _create_input(self): - compiled=subprocess.run(['gcc','-o','gen_dummy_mcpl.out','gen_dummy_mcpl.c','-lm','-lmcpl']) - assert compiled==0, 'Could not compile mcpl-file generator code' - subprocess.run(['./gen_dummy_mcpl.out']) +settings2 = """ + + + 10 + 5 + 1000 + + + source.10.{0} + + +""" - def update_results(self): - """Update the results_true using the current version of OpenMC.""" - try: - self._create_input() - self._test_input_created() - self._run_openmc() - self._test_output_created() - results = self._get_results() - self._write_results(results) - self._overwrite_results() - finally: - self._cleanup() - def _test_input_created(self): - """Check that the input mcpl.file was generated as it should""" - mcplfile=glob.glob(os.path.join(os.getcwd(),'source.10.mcpl')) - assert len(mcplfile) == 1, 'Either multiple or no mcpl files ' \ - 'exist.' - assert mcplfile[0].endswith('mcpl'), \ - 'output file does not end with mcpl.' +class SourceFileTestHarness(TestHarness): + def execute_test(self): + """Run OpenMC with the appropriate arguments and check the outputs.""" + try: + self._run_openmc() + self._test_output_created() + self._run_openmc_restart() + results = self._get_results() + self._write_results(results) + self._compare_results() + finally: + self._cleanup() - def _test_output_created(self): - """Check that the output files were created""" - statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name)) - assert len(statepoint) == 1, 'Either multiple or no statepoint files ' \ - 'exist.' - assert statepoint[0].endswith('h5'), \ - 'statepoint file does not end with h5.' + def update_results(self): + """Update the results_true using the current version of OpenMC.""" + try: + self._run_openmc() + self._test_output_created() + self._run_openmc_restart() + results = self._get_results() + self._write_results(results) + self._overwrite_results() + finally: + self._cleanup() - def _cleanup(self): - super()._cleanup() - source_mcpl=glob.glob(os.path.join(os.getcwd(),'source*.mcpl')) - for f in source_mcpl: - if (os.path.exists(f)): - os.remove(f) + def _test_output_created(self): + """Make sure statepoint and source files have been created.""" + statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name)) + assert len(statepoint) == 1, 'Either multiple or no statepoint files ' \ + 'exist.' + assert statepoint[0].endswith('h5'), \ + 'Statepoint file is not a HDF5 file.' -def test_mcpl_source_file(): - harness = SourceMCPLFileTestHarness('source.10.mcpl') - harness.main() + source = glob.glob(os.path.join(os.getcwd(), 'source.10.mcpl*')) + assert len(source) == 1, 'Either multiple or no source files exist.' + assert source[0].endswith('mcpl') or source[0].endswith('mcpl.gz'), \ + 'Source file is not a MCPL file.' + + def _run_openmc_restart(self): + # Get the name of the source file. + source = glob.glob(os.path.join(os.getcwd(), 'source.10.*')) + + # Write the new settings.xml file. + with open('settings.xml','w') as fh: + fh.write(settings2.format(source[0].split('.')[-1])) + + # Run OpenMC. + self._run_openmc() + + def _cleanup(self): + TestHarness._cleanup(self) + output = glob.glob(os.path.join(os.getcwd(), 'source.*')) + #for f in output: + # if os.path.exists(f): + # os.remove(f) + with open('settings.xml','w') as fh: + fh.write(settings1) + + +def test_source_file(): + harness = SourceFileTestHarness('statepoint.10.h5') + harness.main()