From a3065b479c0a47caae6229fa8657da59668be1e9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 22 Dec 2022 11:33:29 -0600 Subject: [PATCH] Move FileSource MCPL reading to mcpl_interface.h/cpp --- CMakeLists.txt | 1 + include/openmc/mcpl_interface.h | 25 +++++++++ include/openmc/source.h | 10 +--- src/mcpl_interface.cpp | 93 +++++++++++++++++++++++++++++++++ src/settings.cpp | 11 ++-- src/source.cpp | 69 ------------------------ 6 files changed, 126 insertions(+), 83 deletions(-) create mode 100644 include/openmc/mcpl_interface.h create mode 100644 src/mcpl_interface.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 24018eb0f..af3a5b4af 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -331,6 +331,7 @@ list(APPEND libopenmc_SOURCES src/lattice.cpp src/material.cpp src/math_functions.cpp + src/mcpl_interface.cpp src/mesh.cpp src/message_passing.cpp src/mgxs.cpp diff --git a/include/openmc/mcpl_interface.h b/include/openmc/mcpl_interface.h new file mode 100644 index 000000000..89dab16bc --- /dev/null +++ b/include/openmc/mcpl_interface.h @@ -0,0 +1,25 @@ +#ifndef OPENMC_MCPL_INTERFACE_H +#define OPENMC_MCPL_INTERFACE_H + +#include "openmc/particle_data.h" + +#include +#include + +#ifdef OPENMC_MCPL +#include +#endif + +namespace openmc { + +extern "C" const bool MCPL_ENABLED; + +//! 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); + +} // namespace openmc + +#endif // OPENMC_MCPL_INTERFACE_H diff --git a/include/openmc/source.h b/include/openmc/source.h index 051167b3f..96d659dcb 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -14,10 +14,6 @@ #include "openmc/particle.h" #include "openmc/vector.h" -#ifdef OPENMC_MCPL -#include -#endif - namespace openmc { //============================================================================== @@ -32,7 +28,6 @@ constexpr double EXTSRC_REJECT_FRACTION {0.05}; //============================================================================== // Global variables //============================================================================== -extern "C" const bool MCPL_ENABLED; class Source; @@ -107,9 +102,8 @@ class FileSource : public Source { public: // Constructors explicit FileSource(std::string path); -#ifdef OPENMC_MCPL - explicit FileSource(mcpl_file_t mcpl_file); -#endif + explicit FileSource(const vector& sites) : sites_ {sites} {} + // Methods SourceSite sample(uint64_t* seed) const override; diff --git a/src/mcpl_interface.cpp b/src/mcpl_interface.cpp new file mode 100644 index 000000000..d55cc5af7 --- /dev/null +++ b/src/mcpl_interface.cpp @@ -0,0 +1,93 @@ +#include "openmc/mcpl_interface.h" + +#include "openmc/error.h" + +#ifdef OPENMC_MCPL +#include +#endif + +namespace openmc { + +#ifdef OPENMC_MCPL +const bool MCPL_ENABLED = true; +#else +const bool MCPL_ENABLED = false; +#endif + +#ifdef OPENMC_MCPL +SourceSite mcpl_particle_to_site(const mcpl_particle* particle) +{ + SourceSite site; + + switch (particle->pdgcode) { + 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; + } + + // 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]; + + // mcpl stores kinetic energy in MeV + site.E = mcpl_particle->ekin * 1e6; + // mcpl stores time in ms + site.time = mcpl_particle->time * 1e-3; + site.wgt = mcpl_particle->weight; + + return site; +} +#endif + +vector mcpl_source_sites(std::string path) +{ + vector sites; + +#ifdef OPENMC_MCPL + 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; + } + + // Convert to source site and add to vector + sites.push_back(mcpl_particle_to_site(particle)); + } + + // Check that some sites were read + if (sites_.empty()) { + fatal_error("MCPL file contained no neutron, photon, electron, or positron " + "source particles."); + } + + mcpl_close_file(mcpl_file); +#else + fatal_error( + "Your build of OpenMC does not support reading MCPL source files."); +#endif + + return sites; +} + +} // namespace openmc diff --git a/src/settings.cpp b/src/settings.cpp index de298c58d..bc31e2c4d 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -18,6 +18,7 @@ #include "openmc/eigenvalue.h" #include "openmc/error.h" #include "openmc/file_utils.h" +#include "openmc/mcpl_interface.h" #include "openmc/mesh.h" #include "openmc/message_passing.h" #include "openmc/output.h" @@ -431,16 +432,14 @@ void read_settings_xml() for (pugi::xml_node node : root.children("source")) { if (check_for_node(node, "file")) { auto path = get_node_value(node, "file", false, true); -#ifdef OPENMC_MCPL if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) { - model::external_sources.push_back( - make_unique(mcpl_open_file(path.c_str()))); +#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)); } -#else - model::external_sources.push_back(make_unique(path)); -#endif } else if (check_for_node(node, "library")) { // Get shared library path and parameters auto path = get_node_value(node, "library", false, true); diff --git a/src/source.cpp b/src/source.cpp index 4ab087e12..d201e3c03 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -33,22 +33,12 @@ #include "openmc/state_point.h" #include "openmc/xml_interface.h" -#ifdef OPENMC_MCPL -#include -#endif - namespace openmc { //============================================================================== // Global variables //============================================================================== -#ifdef OPENMC_MCPL -const bool MCPL_ENABLED = true; -#else -const bool MCPL_ENABLED = false; -#endif - namespace model { vector> external_sources; @@ -330,65 +320,6 @@ FileSource::FileSource(std::string path) file_close(file_id); } -#ifdef OPENMC_MCPL -FileSource::FileSource(mcpl_file_t mcpl_file) -{ - size_t n_sites = mcpl_hdr_nparticles(mcpl_file); - - for (int i = 0; i < n_sites; i++) { - SourceSite site; - - const mcpl_particle_t* mcpl_particle; - // extract particle from mcpl-file - mcpl_particle = mcpl_read(mcpl_file); - // 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); - pdg = mcpl_particle->pdgcode; - } - - 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; - } - - // 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]; - - // mcpl stores kinetic energy in MeV - site.E = mcpl_particle->ekin * 1e6; - // mcpl stores time in ms - site.time = mcpl_particle->time * 1e-3; - site.wgt = mcpl_particle->weight; - sites_.push_back(site); - } - - // Check that some sites were read - if (sites_.empty()) { - fatal_error("MCPL file contained no neutron, photon, electron, or positron " - "source particles."); - } - - mcpl_close_file(mcpl_file); -} -#endif // OPENMC_MCPL - SourceSite FileSource::sample(uint64_t* seed) const { size_t i_site = sites_.size() * prn(seed);