Move FileSource MCPL reading to mcpl_interface.h/cpp

This commit is contained in:
Paul Romano 2022-12-22 11:33:29 -06:00
parent 987293cfff
commit a3065b479c
6 changed files with 126 additions and 83 deletions

View file

@ -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

View file

@ -0,0 +1,25 @@
#ifndef OPENMC_MCPL_INTERFACE_H
#define OPENMC_MCPL_INTERFACE_H
#include "openmc/particle_data.h"
#include <string>
#include <vector>
#ifdef OPENMC_MCPL
#include <mcpl.h>
#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<SourceSite> mcpl_source_sites(std::string path);
} // namespace openmc
#endif // OPENMC_MCPL_INTERFACE_H

View file

@ -14,10 +14,6 @@
#include "openmc/particle.h"
#include "openmc/vector.h"
#ifdef OPENMC_MCPL
#include <mcpl.h>
#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<SourceSite>& sites) : sites_ {sites} {}
// Methods
SourceSite sample(uint64_t* seed) const override;

93
src/mcpl_interface.cpp Normal file
View file

@ -0,0 +1,93 @@
#include "openmc/mcpl_interface.h"
#include "openmc/error.h"
#ifdef OPENMC_MCPL
#include <mcpl.h>
#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<SourceSite> mcpl_source_sites(std::string path)
{
vector<SourceSite> 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

View file

@ -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<FileSource>(mcpl_open_file(path.c_str())));
#ifdef OPENMC_MCPL
auto sites = mcpl_source_sites(path);
model::external_sources.push_back(make_unique<FileSource>(sites));
#endif
} else {
model::external_sources.push_back(make_unique<FileSource>(path));
}
#else
model::external_sources.push_back(make_unique<FileSource>(path));
#endif
} else if (check_for_node(node, "library")) {
// Get shared library path and parameters
auto path = get_node_value(node, "library", false, true);

View file

@ -33,22 +33,12 @@
#include "openmc/state_point.h"
#include "openmc/xml_interface.h"
#ifdef OPENMC_MCPL
#include <mcpl.h>
#endif
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
#ifdef OPENMC_MCPL
const bool MCPL_ENABLED = true;
#else
const bool MCPL_ENABLED = false;
#endif
namespace model {
vector<unique_ptr<Source>> 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);