Merge pull request #2 from davidjohnlong/fix-source-dev

Review of fixed source bank model developments

include/openmc/particle.h

@@ -181,6 +181,15 @@ public:
     int delayed_group; //!< particle delayed group
   };
 
+  //! Saved ("banked") state of primary source particle, for particle restart
+  struct RestartBank {
+    Position r;
+    Direction u;
+    double E;
+    double wgt;
+  };
+
+
   //==========================================================================
   // Constructors
 
@@ -361,6 +370,8 @@ public:
   #endif
 
   int64_t n_progeny_ {0}; // Number of progeny produced by this particle
+
+  RestartBank primary_source_; // primary source particle data
 };
 
 } // namespace openmc

include/openmc/source.h

@@ -65,8 +65,11 @@ extern "C" void initialize_source();
 //! \return Sampled source site
 Particle::Bank sample_external_source(uint64_t* seed);
 
-//! Fill source bank at end of generation for fixed source simulations
-void fill_source_bank_fixedsource();
+//! Sample a site from custom source library
+Particle::Bank sample_custom_source_library(uint64_t* seed);
+
+void load_custom_source_library();
+void close_custom_source_library();
 
 //! Fill source bank at the end of a generation for dlopen based source simulation
 void fill_source_bank_custom_source();

src/output.cpp

@@ -441,8 +441,6 @@ void print_runtime()
     show_time("Time synchronizing fission bank", time_bank.elapsed(), 1);
     show_time("Sampling source sites", time_bank_sample.elapsed(), 2);
     show_time("SEND/RECV source sites", time_bank_sendrecv.elapsed(), 2);
-  } else {
-    show_time("Time sampling source", time_sample_source.elapsed(), 1);
   }
   show_time("Time accumulating tallies", time_tallies.elapsed(), 1);
   show_time("Total time for finalization", time_finalize.elapsed());

src/particle.cpp

@@ -128,6 +128,13 @@ Particle::from_source(const Bank* src)
     E_ = data::mg.energy_bin_avg_[g_];
   }
   E_last_ = E_;
+
+  // Saved copy of primary source attributes
+  primary_source_.r = src->r;
+  primary_source_.u = src->u;
+  primary_source_.E = E_last_;
+  primary_source_.wgt = wgt_;
+
 }
 
 void
@@ -681,11 +688,10 @@ Particle::write_restart() const
     write_dataset(file_id, "id", id_);
     write_dataset(file_id, "type", static_cast<int>(type_));
 
-    int64_t i = current_work_;
-    write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt);
-    write_dataset(file_id, "energy", simulation::source_bank[i-1].E);
-    write_dataset(file_id, "xyz", simulation::source_bank[i-1].r);
-    write_dataset(file_id, "uvw", simulation::source_bank[i-1].u);
+    write_dataset(file_id, "weight", primary_source_.wgt);
+    write_dataset(file_id, "energy", primary_source_.E);
+    write_dataset(file_id, "xyz", primary_source_.r);
+    write_dataset(file_id, "uvw", primary_source_.u);
 
     // Close file
     file_close(file_id);

src/simulation.cpp

@@ -69,9 +69,10 @@ int openmc_simulation_init()
   // Determine how much work each process should do
   calculate_work();
 
-  // Allocate source bank, and for eigenvalue simulations also allocate the
-  // fission bank
-  allocate_banks();
+  // Allocate source and fission banks for eigenvalue simulations
+  if (settings::run_mode == RunMode::EIGENVALUE) {
+    allocate_banks();
+  }
 
   // If doing an event-based simulation, intialize the particle buffer
   // and event queues
@@ -105,7 +106,16 @@ int openmc_simulation_init()
     load_state_point();
     write_message("Resuming simulation...", 6);
   } else {
-    initialize_source();
+    //Only initialize primary source bank for eigenvalue simulations
+    if (settings::run_mode == RunMode::EIGENVALUE) {
+      initialize_source();
+    }
+  }
+
+  // If fixed source and using custom source library then need to load
+  if (settings::run_mode == RunMode::FIXED_SOURCE &&
+      !settings::path_source_library.empty()) {
+    load_custom_source_library();
   }
 
   // Display header
@@ -154,6 +164,12 @@ int openmc_simulation_finalize()
     t->active_ = false;
   }
 
+  // If fixed source and using custom source library then need to close
+  if (settings::run_mode == RunMode::FIXED_SOURCE &&
+      !settings::path_source_library.empty()) {
+    close_custom_source_library();
+  }
+
   // Stop timers and show timing statistics
   simulation::time_finalize.stop();
   simulation::time_total.stop();
@@ -271,10 +287,8 @@ void allocate_banks()
 {
   // Allocate source bank
   simulation::source_bank.resize(simulation::work_per_rank);
-
-  if (settings::run_mode == RunMode::EIGENVALUE) {
-    init_fission_bank(3*simulation::work_per_rank);
-  }
+  // Allocate fission bank
+  init_fission_bank(3*simulation::work_per_rank);
 }
 
 void initialize_batch()
@@ -350,7 +364,7 @@ void finalize_batch()
       && !settings::cmfd_run) {
     if (contains(settings::sourcepoint_batch, simulation::current_batch)
         && settings::source_write && !settings::source_separate) {
-      bool b = true;
+      bool b = (settings::run_mode == RunMode::EIGENVALUE) ? true : false;
       openmc_statepoint_write(nullptr, &b);
     } else {
       bool b = false;
@@ -358,16 +372,18 @@ void finalize_batch()
     }
   }
 
-  // Write out a separate source point if it's been specified for this batch
-  if (contains(settings::sourcepoint_batch, simulation::current_batch)
-      && settings::source_write && settings::source_separate) {
-    write_source_point(nullptr);
-  }
+  if (settings::run_mode == RunMode::EIGENVALUE) {
+    // Write out a separate source point if it's been specified for this batch
+    if (contains(settings::sourcepoint_batch, simulation::current_batch)
+        && settings::source_write && settings::source_separate) {
+      write_source_point(nullptr);
+    }
 
-  // Write a continously-overwritten source point if requested.
-  if (settings::source_latest) {
-    auto filename = settings::path_output + "source.h5";
-    write_source_point(filename.c_str());
+    // Write a continously-overwritten source point if requested.
+    if (settings::source_latest) {
+      auto filename = settings::path_output + "source.h5";
+      write_source_point(filename.c_str());
+    }
   }
 }
 
@@ -429,18 +445,29 @@ void finalize_generation()
       }
     }
 
-  } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
-    // For fixed-source mode, we need to sample the external source
-    simulation::time_sample_source.start();
-    fill_source_bank_fixedsource();
-    simulation::time_sample_source.stop();
   }
 }
 
 void initialize_history(Particle* p, int64_t index_source)
 {
   // set defaults
-  p->from_source(&simulation::source_bank[index_source - 1]);
+  if(settings::run_mode == RunMode::FIXED_SOURCE) {
+    // initialize random number seed
+    int64_t id = (simulation::total_gen + overall_generation() - 1)*settings::n_particles +
+      simulation::work_index[mpi::rank] + index_source;
+    uint64_t seed = init_seed(id, STREAM_SOURCE);
+    // sample from external source distribution or custom library then set
+    Particle::Bank site;
+    if (!settings::path_source_library.empty()) {
+      site = sample_custom_source_library(&seed);
+    } else {
+      site = sample_external_source(&seed);
+    }
+    p->from_source(&site);
+  } else if (settings::run_mode == RunMode::EIGENVALUE) {
+    // set defaults for eigenvalue simulations from primary bank
+    p->from_source(&simulation::source_bank[index_source - 1]);
+  }
   p->current_work_ = index_source;
 
   // set identifier for particle

src/source.cpp

@@ -38,6 +38,10 @@ namespace openmc {
 
 namespace model {
 
+typedef Particle::Bank (*sample_t)(uint64_t  &seed);
+sample_t sample_source;
+void* source_library;
+
 std::vector<SourceDistribution> external_sources;
 
 }
@@ -276,25 +280,10 @@ void initialize_source()
     file_close(file_id);
   } else if (!settings::path_source_library.empty()) {
 
-  #ifdef HAS_DYNAMIC_LINKING
-    // Get the source from a library object
     write_message(fmt::format("Sampling from library source {}...",
       settings::path_source), 6);
 
-    // Open the library
-    auto source_library = dlopen(settings::path_source_library.c_str(), RTLD_LAZY);
-    if (!source_library) {
-      fatal_error("Couldn't open source library " + settings::path_source_library);
-    }
-
-    // reset errors
-    dlerror();
-
     fill_source_bank_custom_source();
-  #else
-    fatal_error("Custom source libraries have not yet been implemented for "
-      "non-POSIX systems");
-  #endif
 
   } else {
     // Generation source sites from specified distribution in user input
@@ -355,13 +344,14 @@ void free_memory_source()
   model::external_sources.clear();
 }
 
-// fill the source bank from the external source
-void fill_source_bank_custom_source()
+//Load custom source library
+void load_custom_source_library()
 {
 #ifdef HAS_DYNAMIC_LINKING
+
   // Open the library
-  auto source_library = dlopen(settings::path_source_library.c_str(), RTLD_LAZY);
-  if (!source_library) {
+  model::source_library = dlopen(settings::path_source_library.c_str(), RTLD_LAZY);
+  if (!model::source_library) {
     fatal_error("Couldn't open source library " + settings::path_source_library);
   }
 
@@ -369,16 +359,42 @@ void fill_source_bank_custom_source()
   dlerror();
 
   // get the function from the library
-  using sample_t = Particle::Bank (*)(uint64_t* seed);
-  auto sample_source = reinterpret_cast<sample_t>(dlsym(source_library, "sample_source"));
+  //using sample_t = Particle::Bank (*)(uint64_t* seed);
+  model::sample_source = reinterpret_cast<model::sample_t>(dlsym(model::source_library, "sample_source"));
 
   // check for any dlsym errors
   auto dlsym_error = dlerror();
   if (dlsym_error) {
-    dlclose(source_library);
+    dlclose(model::source_library);
     fatal_error(fmt::format("Couldn't open the sample_source symbol: {}", dlsym_error));
   }
 
+#else
+  fatal_error("Custom source libraries have not yet been implemented for "
+    "non-POSIX systems");
+#endif
+
+}
+
+//Release custom source library
+void close_custom_source_library()
+{
+  dlclose(model::source_library);
+}
+
+//Sample source particle from custom library
+Particle::Bank sample_custom_source_library(uint64_t* seed)
+{
+  Particle::Bank site = model::sample_source(*seed);
+  return site;
+}
+
+// fill the source bank from the external source
+void fill_source_bank_custom_source()
+{
+  // Load the custom library
+  load_custom_source_library();
+
   // Generation source sites from specified distribution in the
   // library source
   for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
@@ -387,30 +403,12 @@ void fill_source_bank_custom_source()
       settings::n_particles + simulation::work_index[mpi::rank] + i + 1;
      uint64_t seed = init_seed(id, STREAM_SOURCE);
 
-    // sample external source distribution
-    simulation::source_bank[i] = sample_source(&seed);
+    // sample custom library source
+    simulation::source_bank[i] = sample_custom_source_library(&seed);
   }
 
   // release the library
-  dlclose(source_library);
-#endif
-}
-
-void fill_source_bank_fixedsource()
-{
-  if (settings::path_source.empty() && settings::path_source_library.empty()) {
-    for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
-      // initialize random number seed
-      int64_t id = (simulation::total_gen + overall_generation()) *
-        settings::n_particles + simulation::work_index[mpi::rank] + i + 1;
-      uint64_t seed = init_seed(id, STREAM_SOURCE);
-
-      // sample external source distribution
-      simulation::source_bank[i] = sample_external_source(&seed);
-    }
-  } else if (settings::path_source.empty() && !settings::path_source_library.empty()) {
-     fill_source_bank_custom_source();
-  }
+  close_custom_source_library();
 }
 
 } // namespace openmc

src/state_point.cpp

@@ -380,15 +380,15 @@ void load_state_point()
   // Read batch number to restart at
   read_dataset(file_id, "current_batch", simulation::restart_batch);
 
-  // Check for source in statepoint if needed
-  bool source_present;
-  read_attribute(file_id, "source_present", source_present);
-
   if (simulation::restart_batch > settings::n_batches) {
     fatal_error("The number batches specified in settings.xml is fewer "
       " than the number of batches in the given statepoint file.");
   }
 
+  // Logical flag for source present in statepoint file
+  bool source_present;
+  read_attribute(file_id, "source_present", source_present);
+
   // Read information specific to eigenvalue run
   if (settings::run_mode == RunMode::EIGENVALUE) {
     read_dataset(file_id, "n_inactive", temp);
@@ -396,6 +396,12 @@ void load_state_point()
 
     // Take maximum of statepoint n_inactive and input n_inactive
     settings::n_inactive = std::max(settings::n_inactive, temp);
+
+    // Check to make sure source bank is present
+    if (settings::path_sourcepoint == settings::path_statepoint &&
+        !source_present) {
+      fatal_error("Source bank must be contained in statepoint restart file");
+    }
   }
 
   // Read number of realizations for global tallies
@@ -403,15 +409,13 @@ void load_state_point()
 
   // Set k_sum, keff, and current_batch based on whether restart file is part
   // of active cycle or inactive cycle
-  restart_set_keff();
-  simulation::current_batch = simulation::restart_batch;
-
-  // Check to make sure source bank is present
-  if (settings::path_sourcepoint == settings::path_statepoint &&
-      !source_present) {
-    fatal_error("Source bank must be contained in statepoint restart file");
+  if (settings::run_mode == RunMode::EIGENVALUE) {
+    restart_set_keff();
   }
 
+  // Set current batch number
+  simulation::current_batch = simulation::restart_batch;
+
   // Read tallies to master. If we are using Parallel HDF5, all processes
   // need to be included in the HDF5 calls.
 #ifdef PHDF5