Merge branch 'mcpl_output_dev' of github.com:ebknudsen/openmc into mcpl_output_dev

This commit is contained in:
Erik B Knudsen 2022-10-05 11:01:34 +02:00
commit 0b8bd9b43c

View file

@ -297,17 +297,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:
omc_particle_.particle=ParticleType::neutron;
break;
case 22:
omc_particle_.particle=ParticleType::photon;
break;
case 11:
omc_particle_.particle=ParticleType::electron;
break;
case -11:
omc_particle_.particle=ParticleType::positron;
break;
}
//particle is good, convert to openmc-formalism
@ -390,6 +401,100 @@ CustomSourceWrapper::~CustomSourceWrapper()
#endif
}
<<<<<<< HEAD
<<<<<<< HEAD
#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_);
}
<<<<<<< HEAD
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, 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;
//should check for file exhaustion. This could happen if particles are other than
//neutrons, photons, electrons, or positrons.
}
switch(pdg){
case 2112:
omc_particle_.particle=ParticleType::neutron;
break;
case 22:
omc_particle_.particle=ParticleType::photon;
break;
case 11:
omc_particle_.particle=ParticleType::electron;
break;
case -11:
omc_particle_.particle=ParticleType::positron;
break;
}
//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