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

docs/source/io_formats/settings.rst

@@ -767,6 +767,16 @@ certain surfaces and write out the source bank in a separate file called
 
     *Default*: None
 
+  :mcpl:
+    An optional boolean which indicates if the banked particles should
+    be written to a file in the MCPL-format (documented in mcpl_).
+    instead of the native hdf5-based format.If activated the output
+    output file name is altered to ``surface_source.mcpl``
+
+    *Default*: false
+
+    .. _mcpl: https://mctools.github.io/mcpl/mcpl.pdf
+
 ------------------------------
 ``<survival_biasing>`` Element
 ------------------------------

openmc/settings.py

@@ -1341,13 +1341,15 @@ class Settings:
     def _surf_source_write_from_xml_element(self, root):
         elem = root.find('surf_source_write')
         if elem is not None:
-            for key in ('surface_ids', 'max_particles'):
+            for key in ('surface_ids', 'max_particles','mcpl'):
                 value = get_text(elem, key)
                 if value is not None:
                     if key == 'surface_ids':
                         value = [int(x) for x in value.split()]
                     elif key in ('max_particles'):
                         value = int(value)
+                    elif key == 'mcpl':
+                        value = True
                     self.surf_source_write[key] = value
 
     def _confidence_intervals_from_xml_element(self, root):

src/settings.cpp

@@ -431,12 +431,13 @@ 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);
-      model::external_sources.push_back(make_unique<FileSource>(path));
 #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<FileSource>(mcpl_open_file(path.c_str())));
+      if ( (path.size() >= 5 && path.compare(path.size()-5,5,".mcpl")==0 ) ||
+          (path.size() >= 8 && path.compare(path.size()-8,8,".mcpl.gz")==0 ) ) {
+        model::external_sources.push_back(make_unique<FileSource>(mcpl_open_file(path.c_str())));
+      } else
 #endif
+        model::external_sources.push_back(make_unique<FileSource>(path));
     } else if (check_for_node(node, "library")) {
       // Get shared library path and parameters
       auto path = get_node_value(node, "library", false, true);

src/source.cpp

@@ -234,7 +234,8 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
             auto id = (domain_type_ == DomainType::CELL)
                         ? model::cells[coord.cell]->id_
                         : model::universes[coord.universe]->id_;
-            if (found = contains(domain_ids_, id)) break;
+            if (found = contains(domain_ids_, id))
+              break;
           }
         }
       }
@@ -332,61 +333,57 @@ FileSource::FileSource(std::string path)
 #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);
+  size_t n_sites = mcpl_hdr_nparticles(mcpl_file);
 
   sites_.resize(n_sites);
-  for (int i=0;i<n_sites;i++){
+  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);
+    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.
+    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;
+    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
-    site_.r.x=mcpl_particle->position[0];
-    site_.r.y=mcpl_particle->position[1];
-    site_.r.z=mcpl_particle->position[2];
+    // 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];
+    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 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
+#endif // OPENMC_MCPL
 
 SourceSite FileSource::sample(uint64_t* seed) const
 {
@@ -447,7 +444,6 @@ CustomSourceWrapper::~CustomSourceWrapper()
 #endif
 }
 
-
 //==============================================================================
 // Non-member functions
 //==============================================================================

src/state_point.cpp

@@ -621,10 +621,7 @@ void write_mcpl_source_point(const char *filename, bool surf_source_bank)
 
   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 {

tests/regression_tests/source_mcpl_file/gen_dummy_mcpl.c

@@ -1,98 +0,0 @@
-#include <mcpl.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-
-const int batches=10;
-const int particles=10000;
-
-double rand01(){
-  double r=rand()/((double) RAND_MAX);
-  return r;
-}
-
-
-int main(int argc, char **argv){
-  char outfilename[128];
-  snprintf(outfilename,127,"source.%d.mcpl",batches);
-
-  /*generate an mcpl_header_of sorts*/
-  mcpl_outfile_t outfile=mcpl_create_outfile(outfilename);
-  mcpl_particle_t *particle, Particle;
-
-  char line[256];
-  snprintf(line,255,"gen_dummy_mcpl.c");
-  mcpl_hdr_set_srcname(outfile,line);
-  mcpl_enable_universal_pdgcode(outfile,2112);/*all particles are neutrons*/
-  snprintf(line,255,"Dummy MCPL-file for testing OpenMC mcpl input");
-  mcpl_hdr_add_comment(outfile,line);
-
-  /*also add the instrument file and the command line as blobs*/
-  FILE *fp;
-  if( (fp=fopen("gen_dummy_mcpl.c","rb"))!=NULL){
-    unsigned char *buffer;
-    int size,status;
-    /*find the file size by seeking to end, "tell" the position, and then go back again*/
-    fseek(fp, 0L, SEEK_END);
-    size = ftell(fp); // get current file pointer
-    fseek(fp, 0L, SEEK_SET); // seek back to beginning of file
-    if ( size && (buffer=malloc(size))!=NULL){
-      if (size!=(fread(buffer,1,size,fp))){
-        fprintf(stderr,"Warning: c source generator file not read cleanly\n");
-      }
-      mcpl_hdr_add_data(outfile, "mcpl_point_source_file_generator", size, buffer);
-      free(buffer);
-    }
-    fclose(fp);
-  } else {
-    fprintf(stderr,"Warning: could not open c source generator file, hence not embedded.\n");
-  }
-
-
-  /*the main particle loop*/
-  particle=&Particle;
-  int i;
-  srand(1234);
-  for (i=0;i<batches*particles; i++){
-    if(i%particles==0){
-      printf("MCPL_test: generating batch number %d\n",i/particles);
-    }
-    /*generate a point source - i.e. x,y,z are 0*/
-    particle->position[0]=0;
-    particle->position[1]=0;
-    particle->position[2]=0;
-
-    /*generate a random direction on unit sphere*/
-    double nrm=2.0;
-    double vx,vy,vz;
-    int iter=0;
-    do {
-      if(iter>100){
-        printf("Warning: Exceeded max random vector attempts: at loop %d -> %g   %d\n",i,nrm,iter);
-      }
-      vx=rand01()*2.0-1.0;
-      vy=rand01()*2.0-1.0;
-      vz=rand01()*2.0-1.0;
-      nrm=(vx*vx + vy*vy + vz*vz);
-      iter++;
-    } while (nrm>1.0);
-    vx/=sqrt(nrm);
-    vy/=sqrt(nrm);
-    vz/=sqrt(nrm);
-    particle->direction[0]=vx;
-    particle->direction[1]=vy;
-    particle->direction[2]=vz;
-
-    /*generate the kinetic energy (in MeV) of the particle*/
-    particle->ekin=1e-3;
-
-    /*set time=0*/
-    particle->time=0;
-    /*set weight*/
-    particle->weight=1;
-
-    particle->userflags=0;
-    mcpl_add_particle(outfile,particle);
-  }
-  mcpl_close_outfile(outfile);
-}