diff --git a/include/openmc/source.h b/include/openmc/source.h index 4db17f241a..97b18d2286 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -65,8 +65,14 @@ 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); + +//! Load custom source library +void load_custom_source_library(); + +//! Release 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(); diff --git a/src/output.cpp b/src/output.cpp index 18b94be9be..8da32de55a 100644 --- a/src/output.cpp +++ b/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()); diff --git a/src/particle.cpp b/src/particle.cpp index 5f9713941e..96cf3fdcf8 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -21,6 +21,7 @@ #include "openmc/physics_mg.h" #include "openmc/random_lcg.h" #include "openmc/settings.h" +#include "openmc/source.h" #include "openmc/surface.h" #include "openmc/simulation.h" #include "openmc/tallies/derivative.h" @@ -682,10 +683,24 @@ Particle::write_restart() const write_dataset(file_id, "type", static_cast(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); + if (settings::run_mode == RunMode::EIGENVALUE) { + //take source data from primary bank for eigenvalue simulation + 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); + } else if (settings::run_mode == RunMode::FIXED_SOURCE) { + // re-sample using rng random number seed used to generate source particle + int64_t id = (simulation::total_gen + overall_generation() - 1)*settings::n_particles + + simulation::work_index[mpi::rank] + i; + uint64_t seed = init_seed(id, STREAM_SOURCE); + // re-sample source site + auto site = sample_external_source(&seed); + write_dataset(file_id, "weight", site.wgt); + write_dataset(file_id, "energy", site.E); + write_dataset(file_id, "xyz", site.r); + write_dataset(file_id, "uvw", site.u); + } // Close file file_close(file_id); diff --git a/src/simulation.cpp b/src/simulation.cpp index 086d740593..64a7d0e806 100644 --- a/src/simulation.cpp +++ b/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); 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,24 @@ 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 + auto 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 diff --git a/src/source.cpp b/src/source.cpp index 96175804ce..21b12ed0dc 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -42,6 +42,15 @@ std::vector external_sources; } +namespace { + +using sample_t = Particle::Bank (*)(uint64_t* seed); +sample_t custom_source_function; +void* custom_source_library; + +} + + //============================================================================== // SourceDistribution implementation //============================================================================== @@ -276,25 +285,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 @@ -321,6 +315,11 @@ void initialize_source() Particle::Bank sample_external_source(uint64_t* seed) { + // return values from custom source if using + if (!settings::path_source_library.empty()) { + return sample_custom_source_library(seed); + } + // Determine total source strength double total_strength = 0.0; for (auto& s : model::external_sources) @@ -355,13 +354,13 @@ void free_memory_source() model::external_sources.clear(); } -// fill the source bank from the external source -void fill_source_bank_custom_source() +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) { + custom_source_library = dlopen(settings::path_source_library.c_str(), RTLD_LAZY); + if (!custom_source_library) { fatal_error("Couldn't open source library " + settings::path_source_library); } @@ -370,15 +369,36 @@ void fill_source_bank_custom_source() // get the function from the library using sample_t = Particle::Bank (*)(uint64_t* seed); - auto sample_source = reinterpret_cast(dlsym(source_library, "sample_source")); + custom_source_function = reinterpret_cast(dlsym(custom_source_library, "sample_source")); // check for any dlsym errors auto dlsym_error = dlerror(); if (dlsym_error) { - dlclose(source_library); + dlclose(custom_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 +} + +void close_custom_source_library() +{ + dlclose(custom_source_library); +} + +Particle::Bank sample_custom_source_library(uint64_t* seed) +{ + return custom_source_function(seed); +} + +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 +407,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 diff --git a/src/state_point.cpp b/src/state_point.cpp index d077590365..c8c3c87bd8 100644 --- a/src/state_point.cpp +++ b/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