wip add MCPL-out to openmc

code for writing a source_bank/point compiles but is untested

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

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)
@@ -598,6 +602,45 @@ void write_source_point(const char* filename, bool surf_source_bank)
     file_close(file_id);
 }
 
+void write_mcpl_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 {VERSION_MAJOR}.{VERSION_MINOR}.{VERSION_RELEASE}-development");
+    } else {
+      line=fmt::format("OpenMC {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);
+  }
+
+}
+
+
 void write_source_bank(hid_t group_id, bool surf_source_bank)
 {
   hid_t banktype = h5banktype();
@@ -714,6 +757,88 @@ void write_source_bank(hid_t group_id, bool surf_source_bank)
   H5Tclose(banktype);
 }
 
+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) {
+    //write particles from the master node
+    //loop over the other nodes and receive data - then write those.
+    for (int i = 0; i < mpi::n_procs; ++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_banke 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;
+
+        p.direction[0]=_site->u.x;
+        p.direction[1]=_site->u.y;
+        p.direction[2]=_site->u.x;
+
+        //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);
+      }
+    }
+  } else {
+#ifdef OPENMC_MPI
+    MPI_Send(source_bank->data(), count_size, mpi::source_site, 0, mpi::rank,
+      mpi::intracomm);
+#endif
+  }
+
+}
+
+
 // Determine member names of a compound HDF5 datatype
 std::string dtype_member_names(hid_t dtype_id)
 {