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

CMakeLists.txt

@@ -36,7 +36,7 @@ option(OPENMC_ENABLE_COVERAGE "Compile with coverage analysis flags"
 option(OPENMC_USE_DAGMC       "Enable support for DAGMC (CAD) geometry"              OFF)
 option(OPENMC_USE_LIBMESH     "Enable support for libMesh unstructured mesh tallies" OFF)
 option(OPENMC_USE_MPI         "Enable MPI"                                           OFF)
-option(OPENMC_USE_MCPL        "Enable MPCPL"                                         OFF)
+option(OPENMC_USE_MCPL        "Enable MCPL"                                          OFF)
 
 #===============================================================================
 # Set a default build configuration if not explicitly specified

include/openmc/source.h

@@ -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
 //==============================================================================

src/settings.cpp

@@ -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

src/source.cpp

@@ -283,6 +283,65 @@ 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, 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
+    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);
@@ -342,6 +401,7 @@ CustomSourceWrapper::~CustomSourceWrapper()
 #endif
 }
 
+<<<<<<< HEAD
 #ifdef OPENMC_MCPL
 //===========================================================================
 // Read particles from an MCPL-file
@@ -433,6 +493,8 @@ SourceSite MCPLFileSource::read_single_particle() const
   return omc_particle_;
 }
 #endif //OPENMC_MCPL
+=======
+>>>>>>> 380cfc150b449ef56debb5261caa010a52c42bdd
 
 //==============================================================================
 // Non-member functions