refactor MCPL file read part of FileSource

A new constructor is added to FileSource that reads from
a given MCPL-file.
This commit is contained in:
erkn 2022-08-09 08:12:19 +02:00
parent 1ba6967414
commit 380cfc150b
3 changed files with 52 additions and 108 deletions

View file

@ -96,7 +96,9 @@ class FileSource : public Source {
public:
// Constructors
explicit FileSource(std::string path);
#ifdef OPENMC_MCPL
explicit FileSource(mcpl_file_t mcpl_file);
#endif
// Methods
SourceSite sample(uint64_t* seed) const override;
@ -104,32 +106,6 @@ private:
vector<SourceSite> sites_; //!< Source sites from a file
};
#ifdef OPENMC_MCPL
//==============================================================================
// MCPL-file input source
//==============================================================================
class MCPLFileSource : public Source {
public:
// Constructors, destructors
MCPLFileSource(std::string path);
~MCPLFileSource();
// Methods
//! Sample from the external source distribution
//! \param[inout] seed Pseudorandom seed pointer
//! \return Site read from MCPL-file
SourceSite sample(uint64_t* seed) const;
private:
SourceSite read_single_particle() const;
void read_source_bank(vector<SourceSite> &sites_);
vector <SourceSite> sites_; //!<source sites from an MCPL-file
mcpl_file_t mcpl_file;
uint64_t n_sites;
};
#endif
//==============================================================================
//! Wrapper for custom sources that manages opening/closing shared library
//==============================================================================

View file

@ -433,7 +433,7 @@ void read_settings_xml()
#ifdef OPENMC_MCPL
} else if (check_for_node(node, "mcpl")) {
auto path = get_node_value(node, "mcpl", false, true);
model::external_sources.push_back(make_unique<MCPLFileSource>(path));
model::external_sources.push_back(make_unique<FileSource>(mcpl_open_file(path.c_str())));
#endif
} else if (check_for_node(node, "library")) {
// Get shared library path and parameters

View file

@ -277,6 +277,54 @@ FileSource::FileSource(std::string path)
file_close(file_id);
}
#ifdef OPENMC_MCPL
FileSource::FileSource(mcpl_file_t mcpl_file)
{
//do checks on the mcpl_file to see if particles are many enough.
// should model this on the example source shown in the docs.
size_t n_sites=mcpl_hdr_nparticles(mcpl_file);
sites_.resize(n_sites);
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 or a photon. otherwise skip
while ( mcpl_particle->pdgcode!=2112 && mcpl_particle->pdgcode!=22 ) {
mcpl_particle=mcpl_read(mcpl_file);
//should check for file exhaustion This could happen if particles are other than
//neutrons or photons
}
if(mcpl_particle->pdgcode==2112) {
site_.particle=ParticleType::neutron;
} else if (mcpl_particle->pdgcode==22) {
site_.particle=ParticleType::photon;
}
//particle is good, convert to openmc-formalism
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;
site_.delayed_group=0;
sites_[i]=site_;
}
mcpl_close_file(mcpl_file);
}
#endif //OPENMC_MCPL
SourceSite FileSource::sample(uint64_t* seed) const
{
size_t i_site = sites_.size() * prn(seed);
@ -336,86 +384,6 @@ CustomSourceWrapper::~CustomSourceWrapper()
#endif
}
#ifdef OPENMC_MCPL
//===========================================================================
// Read particles from an MCPL-file
//===========================================================================
MCPLFileSource::MCPLFileSource(std::string path)
{
// Check if source file exists
if (!file_exists(path)) {
fatal_error(fmt::format("Source file '{}' does not exist.", path));
}
// Read the source from a binary file instead of sampling from some
// assumed source distribution
write_message(6, "Reading mcpl source file from {}",path);
// Open the mcpl file
mcpl_file = mcpl_open_file(path.c_str());
//do checks on the mcpl_file to see if particles are many enough.
// should model this on the example source shown in the docs.
n_sites=mcpl_hdr_nparticles(mcpl_file);
read_source_bank(sites_);
}
MCPLFileSource::~MCPLFileSource(){
mcpl_close_file(mcpl_file);
}
SourceSite MCPLFileSource::sample(uint64_t* seed) const
{
size_t i_site = sites_.size() * prn(seed);
return sites_[i_site];
}
void MCPLFileSource::read_source_bank(vector<SourceSite> &sites_)
{
sites_.resize(n_sites);
for (int i=0;i<n_sites;i++){
sites_[i]=read_single_particle();
}
}
SourceSite MCPLFileSource::read_single_particle() const
{
SourceSite omc_particle_;
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 ) {
mcpl_particle=mcpl_read(mcpl_file);
//should check for file exhaustion This could happen if particles are other than
//neutrons or photons
}
if(mcpl_particle->pdgcode==2112) {
omc_particle_.particle=ParticleType::neutron;
} else if (mcpl_particle->pdgcode==22) {
omc_particle_.particle=ParticleType::photon;
}
//particle is good, convert to openmc-formalism
omc_particle_.r.x=mcpl_particle->position[0];
omc_particle_.r.y=mcpl_particle->position[1];
omc_particle_.r.z=mcpl_particle->position[2];
omc_particle_.u.x=mcpl_particle->direction[0];
omc_particle_.u.y=mcpl_particle->direction[1];
omc_particle_.u.z=mcpl_particle->direction[2];
//mcpl stores kinetic energy in MeV
omc_particle_.E=mcpl_particle->ekin*1e6;
//mcpl stores time in ms
omc_particle_.time=mcpl_particle->time*1e-3;
omc_particle_.wgt=mcpl_particle->weight;
omc_particle_.delayed_group=0;
return omc_particle_;
}
#endif //OPENMC_MCPL
//==============================================================================
// Non-member functions