Store sites from source file in SourceDistribution

include/openmc/settings.h

@@ -59,7 +59,6 @@ extern std::string path_cross_sections;   //!< path to cross_sections.xml
 extern std::string path_input;            //!< directory where main .xml files resides
 extern std::string path_output;           //!< directory where output files are written
 extern std::string path_particle_restart; //!< path to a particle restart file
-extern std::string path_source;
 extern std::string path_source_library;   //!< path to the source shared object
 extern std::string path_sourcepoint;      //!< path to a source file
 extern "C" std::string path_statepoint;   //!< path to a statepoint file

include/openmc/source.h

@@ -57,6 +57,7 @@ private:
   UPtrSpace space_; //!< Spatial distribution
   UPtrAngle angle_; //!< Angular distribution
   UPtrDist energy_; //!< Energy distribution
+  std::vector<Particle::Bank> sites_; //!< Source sites from a file
 };
 
 class CustomSource {

include/openmc/state_point.h

@@ -2,17 +2,19 @@
 #define OPENMC_STATE_POINT_H
 
 #include <cstdint>
+#include <vector>
 
 #include "hdf5.h"
 
 #include "openmc/capi.h"
+#include "openmc/particle.h"
 
 namespace openmc {
 
 void load_state_point();
 void write_source_point(const char* filename);
 void write_source_bank(hid_t group_id);
-void read_source_bank(hid_t group_id);
+void read_source_bank(hid_t group_id, std::vector<Particle::Bank>& sites);
 void write_tally_results_nr(hid_t file_id);
 void restart_set_keff();
 

src/settings.cpp

@@ -73,7 +73,6 @@ std::string path_cross_sections;
 std::string path_input;
 std::string path_output;
 std::string path_particle_restart;
-std::string path_source;
 std::string path_source_library;
 std::string path_sourcepoint;
 std::string path_statepoint;

src/source.cpp

@@ -82,13 +82,35 @@ SourceDistribution::SourceDistribution(pugi::xml_node node)
   // Check for external source file
   if (check_for_node(node, "file")) {
     // Copy path of source file
-    settings::path_source = get_node_value(node, "file", false, true);
+    auto path_source = get_node_value(node, "file", false, true);
 
     // Check if source file exists
-    if (!file_exists(settings::path_source)) {
-      fatal_error(fmt::format("Source file '{}' does not exist.",
-        settings::path_source));
+    if (!file_exists(path_source)) {
+      fatal_error(fmt::format("Source file '{}' does not exist.", path_source));
     }
+
+    // Read the source from a binary file instead of sampling from some
+    // assumed source distribution
+    write_message(6, "Reading source file from {}...", path_source);
+
+    // Open the binary file
+    hid_t file_id = file_open(path_source, 'r', true);
+
+    // Read the file type
+    std::string filetype;
+    read_attribute(file_id, "filetype", filetype);
+
+    // Check to make sure this is a source file
+    if (filetype != "source" && filetype != "statepoint") {
+      fatal_error("Specified starting source file not a source file type.");
+    }
+
+    // Read in the source bank
+    read_source_bank(file_id, sites_);
+
+    // Close file
+    file_close(file_id);
+
   } else if (check_for_node(node, "library")) {
     settings::path_source_library = get_node_value(node, "library", false, true);
     if (!file_exists(settings::path_source_library)) {
@@ -180,15 +202,20 @@ Particle::Bank SourceDistribution::sample(uint64_t* seed) const
   int n_reject = 0;
   static int n_accept = 0;
   while (!found) {
-    // Set particle type
-    site.particle = particle_;
+    if (!sites_.empty()) {
+      size_t i_site = sites_.size()*prn(seed);
+      site = sites_[i_site];
+    } else {
+      // Set particle type
+      site.particle = particle_;
 
-    // Sample spatial distribution
-    site.r = space_->sample(seed);
-    double xyz[] {site.r.x, site.r.y, site.r.z};
+      // Sample spatial distribution
+      site.r = space_->sample(seed);
+    }
 
     // Now search to see if location exists in geometry
     int32_t cell_index, instance;
+    double xyz[] {site.r.x, site.r.y, site.r.z};
     int err = openmc_find_cell(xyz, &cell_index, &instance);
     found = (err != OPENMC_E_GEOMETRY);
 
@@ -225,6 +252,8 @@ Particle::Bank SourceDistribution::sample(uint64_t* seed) const
   // Increment number of accepted samples
   ++n_accept;
 
+  if (!sites_.empty()) return site;
+
   // Sample angle
   site.u = angle_->sample(seed);
 
@@ -264,33 +293,8 @@ void initialize_source()
 {
   write_message("Initializing source particles...", 5);
 
-  if (!settings::path_source.empty()) {
-    // Read the source from a binary file instead of sampling from some
-    // assumed source distribution
-
-    write_message(6, "Reading source file from {}...", settings::path_source);
-
-    // Open the binary file
-    hid_t file_id = file_open(settings::path_source, 'r', true);
-
-    // Read the file type
-    std::string filetype;
-    read_attribute(file_id, "filetype", filetype);
-
-    // Check to make sure this is a source file
-    if (filetype != "source" && filetype != "statepoint") {
-      fatal_error("Specified starting source file not a source file type.");
-    }
-
-    // Read in the source bank
-    read_source_bank(file_id);
-
-    // Close file
-    file_close(file_id);
-  } else if (!settings::path_source_library.empty()) {
-
-    write_message(6, "Sampling library source {}...", settings::path_source);
-
+  if (!settings::path_source_library.empty()) {
+    write_message(6, "Sampling library source {}...", settings::path_source_library);
     fill_source_bank_custom_source();
 
   } else {

src/state_point.cpp

@@ -485,11 +485,11 @@ void load_state_point()
         + "...", 5);
 
       // Open source file
-      file_id = file_open(settings::path_source.c_str(), 'r', true);
+      file_id = file_open(settings::path_sourcepoint.c_str(), 'r', true);
     }
 
     // Read source
-    read_source_bank(file_id);
+    read_source_bank(file_id, simulation::source_bank);
 
   }
 
@@ -653,7 +653,7 @@ write_source_bank(hid_t group_id)
 }
 
 
-void read_source_bank(hid_t group_id)
+void read_source_bank(hid_t group_id, std::vector<Particle::Bank>& sites)
 {
   hid_t banktype = h5banktype();
 
@@ -666,11 +666,15 @@ void read_source_bank(hid_t group_id)
 
   // Make sure source bank is big enough
   hid_t dspace = H5Dget_space(dset);
-  hsize_t dims_all[1];
-  H5Sget_simple_extent_dims(dspace, dims_all, nullptr);
-  if (simulation::work_index[mpi::n_procs] > dims_all[0]) {
-    fatal_error("Number of source sites in source file is less "
-                "than number of source particles per generation.");
+  hsize_t dims_all;
+  H5Sget_simple_extent_dims(dspace, &dims_all, nullptr);
+  if (&sites == &simulation::source_bank) {
+    if (simulation::work_index[mpi::n_procs] > dims_all) {
+      fatal_error("Number of source sites in source file is less "
+                  "than number of source particles per generation.");
+    }
+  } else {
+    sites.resize(dims_all);
   }
 
   // Select hyperslab for each process
@@ -681,10 +685,10 @@ void read_source_bank(hid_t group_id)
     // Read data in parallel
   hid_t plist = H5Pcreate(H5P_DATASET_XFER);
   H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
-  H5Dread(dset, banktype, memspace, dspace, plist, simulation::source_bank.data());
+  H5Dread(dset, banktype, memspace, dspace, plist, sites.data());
   H5Pclose(plist);
 #else
-  H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, simulation::source_bank.data());
+  H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, sites.data());
 #endif
 
   // Close all ids