Merge branch 'mcpl_output_dev' into mcpl_input

CMakeLists.txt

@@ -159,6 +159,13 @@ if(${HDF5_VERSION} VERSION_GREATER_EQUAL 1.12.0)
   list(APPEND cxxflags -DH5Oget_info_by_idx_vers=1 -DH5O_info_t_vers=1)
 endif()
 
+#===============================================================================
+# MCPL
+#===============================================================================
+if (OPENMC_USE_MCPL)
+  find_package(MCPL REQUIRED)
+endif()
+
 #===============================================================================
 # Set compile/link flags based on which compiler is being used
 #===============================================================================
@@ -186,10 +193,6 @@ if(OPENMC_ENABLE_COVERAGE)
   list(APPEND ldflags --coverage)
 endif()
 
-if(OPENMC_USE_MCPL)
-  list(APPEND ldflags -lmcpl)
-endif()
-
 # Show flags being used
 message(STATUS "OpenMC C++ flags: ${cxxflags}")
 message(STATUS "OpenMC Linker flags: ${ldflags}")
@@ -444,10 +447,6 @@ endif()
 if (OPENMC_USE_MPI)
   target_compile_definitions(libopenmc PUBLIC -DOPENMC_MPI)
 endif()
-if(OPENMC_USE_MCPL)
-  target_compile_definitions(libopenmc PUBLIC -DOPENMC_MCPL)
-endif()
-
 
 # Set git SHA1 hash as a compile definition
 if(GIT_FOUND)
@@ -491,6 +490,11 @@ if (OPENMC_USE_MPI)
   target_link_libraries(libopenmc MPI::MPI_CXX)
 endif()
 
+if (OPENMC_USE_MCPL)
+  target_compile_definitions(libopenmc PUBLIC OPENMC_MCPL)
+  target_link_libraries(libopenmc MCPL::mcpl)
+endif()
+
 #===============================================================================
 # openmc executable
 #===============================================================================

include/openmc/settings.h

@@ -48,6 +48,7 @@ extern bool source_latest;         //!< write latest source at each batch?
 extern bool source_separate;       //!< write source to separate file?
 extern bool source_write;          //!< write source in HDF5 files?
 extern bool surf_source_write;     //!< write surface source file?
+extern bool surf_mcpl_write;       //!< write surface mcpl file?
 extern bool surf_source_read;      //!< read surface source file?
 extern bool survival_biasing;      //!< use survival biasing?
 extern bool temperature_multipole; //!< use multipole data?

include/openmc/state_point.h

@@ -9,6 +9,10 @@
 #include "openmc/particle.h"
 #include "openmc/vector.h"
 
+#ifdef OPENMC_MCPL
+#include <mcpl.h>
+#endif
+
 namespace openmc {
 
 void load_state_point();
@@ -21,5 +25,10 @@ void write_tally_results_nr(hid_t file_id);
 void restart_set_keff();
 void write_unstructured_mesh_results();
 
+#ifdef OPENMC_MCPL
+void write_mcpl_source_point(const char *filename_, bool surf_source_bank = false);
+void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank);
+#endif
+
 } // namespace openmc
 #endif // OPENMC_STATE_POINT_H

openmc/settings.py

@@ -171,6 +171,7 @@ class Settings:
                    banked (int)
         :max_particles: Maximum number of particles to be banked on
                    surfaces per process (int)
+        :mcpl: Output in the form of an MCPL-file (bool)
     survival_biasing : bool
         Indicate whether survival biasing is to be used
     tabular_legendre : dict
@@ -652,7 +653,7 @@ class Settings:
         cv.check_type('surface source writing options', surf_source_write, Mapping)
         for key, value in surf_source_write.items():
             cv.check_value('surface source writing key', key,
-                           ('surface_ids', 'max_particles'))
+                           ('surface_ids', 'max_particles', 'mcpl'))
             if key == 'surface_ids':
                 cv.check_type('surface ids for source banking', value,
                               Iterable, Integral)
@@ -664,6 +665,9 @@ class Settings:
                               value, Integral)
                 cv.check_greater_than('maximum particle banks on surfaces per process',
                                       value, 0)
+            elif key == 'mcpl':
+                cv.check_type('write to an MCPL-format file', value, bool)
+
         self._surf_source_write = surf_source_write
 
     @confidence_intervals.setter
@@ -1032,6 +1036,9 @@ class Settings:
             if 'max_particles' in self._surf_source_write:
                 subelement = ET.SubElement(element, "max_particles")
                 subelement.text = str(self._surf_source_write['max_particles'])
+            if 'mcpl' in self._surf_source_write:
+                subelement = ET.SubElement(element, "mcpl")
+                subelement.text = str(self._surf_source_write['mcpl']).lower()
 
     def _create_confidence_intervals(self, root):
         if self._confidence_intervals is not None:

openmc/source.py

@@ -223,7 +223,10 @@ class Source:
         if self.particle != 'neutron':
             element.set("particle", self.particle)
         if self.file is not None:
-            element.set("file", self.file)
+            if (self.file.endswith('.mcpl')):
+              element.set("mcpl", self.file)
+            else:
+              element.set("file", self.file)
         if self.library is not None:
             element.set("library", self.library)
         if self.parameters is not None:

src/settings.cpp

@@ -61,6 +61,7 @@ bool source_latest {false};
 bool source_separate {false};
 bool source_write {true};
 bool surf_source_write {false};
+bool surf_mcpl_write {false};
 bool surf_source_read {false};
 bool survival_biasing {false};
 bool temperature_multipole {false};
@@ -682,9 +683,14 @@ void read_settings_xml()
       max_surface_particles =
         std::stoll(get_node_value(node_ssw, "max_particles"));
     }
+#ifdef OPENMC_MCPL
+    if(check_for_node(node_ssw, "mcpl")){
+      surf_mcpl_write=true;
+    }
+#endif
   }
 
-  // If source is not seperate and is to be written out in the statepoint file,
+  // If source is not separate and is to be written out in the statepoint file,
   // make sure that the sourcepoint batch numbers are contained in the
   // statepoint list
   if (!source_separate) {

src/simulation.cpp

@@ -407,6 +407,13 @@ void finalize_batch()
     auto filename = settings::path_output + "surface_source.h5";
     write_source_point(filename.c_str(), true);
   }
+#ifdef OPENMC_MCPL
+  if(settings::surf_mcpl_write && simulation::current_batch == settings::n_batches){
+    auto filename = settings::path_output + "surface_source.mcpl";
+    write_mcpl_source_point(filename.c_str(), true);
+  }
+#endif
+
 }
 
 void initialize_generation()

src/source.cpp

@@ -343,17 +343,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:
+        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

src/state_point.cpp

@@ -28,6 +28,10 @@
 #include "openmc/timer.h"
 #include "openmc/vector.h"
 
+#ifdef OPENMC_MCPL
+#include <mcpl.h>
+#endif
+
 namespace openmc {
 
 extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
@@ -599,6 +603,46 @@ void write_source_point(const char* filename, bool surf_source_bank)
     file_close(file_id);
 }
 
+#ifdef OPENMC_MCPL
+void write_mcpl_source_point(const char *filename, bool surf_source_bank)
+{
+  std::string filename_;
+  if (filename) {
+    filename_ = filename;
+  } else {
+    // Determine width for zero padding
+    int w = std::to_string(settings::n_max_batches).size();
+
+    filename_ = fmt::format("{0}source.{1:0{2}}.mcpl", settings::path_output,
+      simulation::current_batch, w);
+  }
+
+  mcpl_outfile_t file_id;
+
+  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 {
+      line=fmt::format("OpenMC {0}.{1}.{2}",VERSION_MAJOR,VERSION_MINOR,VERSION_RELEASE);
+    }
+    mcpl_hdr_set_srcname(file_id,line.c_str());
+  }
+
+  write_mcpl_source_bank(file_id, surf_source_bank);
+
+  if (mpi::master) {
+    //change this - this is h5 specific
+    mcpl_closeandgzip_outfile(file_id);
+  }
+
+}
+#endif
+
 void write_source_bank(hid_t group_id, bool surf_source_bank)
 {
   hid_t banktype = h5banktype();
@@ -715,6 +759,105 @@ void write_source_bank(hid_t group_id, bool surf_source_bank)
   H5Tclose(banktype);
 }
 
+#ifdef OPENMC_MCPL
+void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank)
+{
+  int64_t dims_size = settings::n_particles;
+  int64_t count_size = simulation::work_per_rank;
+
+  // Set vectors for source bank and starting bank index of each process
+  vector<int64_t>* bank_index = &simulation::work_index;
+  vector<SourceSite>* source_bank = &simulation::source_bank;
+  vector<int64_t> surf_source_index_vector;
+  vector<SourceSite> surf_source_bank_vector;
+
+  if(surf_source_bank) {
+    surf_source_index_vector = calculate_surf_source_size();
+    dims_size = surf_source_index_vector[mpi::n_procs];
+    count_size = simulation::surf_source_bank.size();
+
+    bank_index = &surf_source_index_vector;
+
+    // Copy data in a SharedArray into a vector.
+    surf_source_bank_vector.resize(count_size);
+    surf_source_bank_vector.assign(simulation::surf_source_bank.data(),
+      simulation::surf_source_bank.data() + count_size);
+    source_bank = &surf_source_bank_vector;
+  }
+
+  if (mpi::master) {
+    // Particles are writeen to disk from the master node only
+
+    // Save source bank sites since the array is overwritten below
+#ifdef OPENMC_MPI
+    vector<SourceSite> temp_source {source_bank->begin(), source_bank->end()};
+#endif
+
+    //loop over the other nodes and receive data - then write those.
+    for (int i = 0; i < mpi::n_procs; ++i) {
+      // number of particles for node node i
+      size_t count[] {
+        static_cast<size_t>((*bank_index)[i + 1] - (*bank_index)[i])};
+
+#ifdef OPENMC_MPI
+      if (i>0)
+        MPI_Recv(source_bank->data(), count[0], mpi::source_site, i, i,
+          mpi::intracomm, MPI_STATUS_IGNORE);
+#endif
+      // now write the source_bank data again.
+      for (vector<SourceSite>::iterator _site=source_bank->begin(); _site!=source_bank->end();_site++){
+        // particle is now at the iterator
+        // write it to the mcpl-file
+        mcpl_particle_t p;
+        p.position[0]=_site->r.x;
+        p.position[1]=_site->r.y;
+        p.position[2]=_site->r.z;
+
+        // mcpl requires that the direction vector is unit length
+        // which is also the case in openmc
+        p.direction[0]=_site->u.x;
+        p.direction[1]=_site->u.y;
+        p.direction[2]=_site->u.z;
+
+        // mcpl stores kinetic energy in MeV
+        p.ekin=_site->E*1e-6;
+
+        p.time=_site->time*1e3;
+
+        p.weight=_site->wgt;
+
+        switch(_site->particle){
+          case ParticleType::neutron:
+            p.pdgcode=2112;
+            break;
+          case ParticleType::photon:
+            p.pdgcode=22;
+            break;
+          case ParticleType::electron:
+            p.pdgcode=11;
+            break;
+          case ParticleType::positron:
+            p.pdgcode=-11;
+            break;
+        }
+
+        mcpl_add_particle(file_id,&p);
+      }
+    }
+#ifdef OPENMC_MPI
+    // Restore state of source bank
+    std::copy(temp_source.begin(), temp_source.end(), source_bank->begin());
+#endif
+  } else {
+#ifdef OPENMC_MPI
+    MPI_Send(source_bank->data(), count_size, mpi::source_site, 0, mpi::rank,
+      mpi::intracomm);
+#endif
+  }
+
+}
+#endif
+
 // Determine member names of a compound HDF5 datatype
 std::string dtype_member_names(hid_t dtype_id)
 {