diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index f89bdc8bf..d29f404b9 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -83,6 +83,10 @@ extern "C" double ufs_get_weight(const Particle* p); //! \param[in] group HDF5 group void write_eigenvalue_hdf5(hid_t group); +//! Read data related to k-eigenvalue from statepoint +//! \param[in] group HDF5 group +void read_eigenvalue_hdf5(hid_t group); + } // namespace openmc #endif // OPENMC_EIGENVALUE_H diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 6881cac77..398f03b09 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -229,6 +229,19 @@ void read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false) read_dataset(obj_id, name, H5TypeMap::type_id, &buffer, indep); } +// overload for std::string +inline void +read_dataset(hid_t obj_id, const char* name, std::string& str, bool indep=false) +{ + // Create buffer to read data into + auto n = attribute_typesize(obj_id, name); + char buffer[n]; + + // Read attribute and set string + read_string(obj_id, name, n, buffer, indep); + str = std::string{buffer, n}; +} + template void read_dataset(hid_t dset, std::vector& vec, bool indep=false) { diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index 0c1db86f3..e44c57f1a 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -9,11 +9,12 @@ namespace openmc { +void load_state_point(); void write_source_point(const char* filename); -extern "C" void write_source_bank(hid_t group_id); -extern "C" void read_source_bank(hid_t group_id); -extern "C" void write_tally_results_nr(hid_t file_id); -extern "C" void restart_set_keff(); +void write_source_bank(hid_t group_id); +void read_source_bank(hid_t group_id); +void write_tally_results_nr(hid_t file_id); +void restart_set_keff(); } // namespace openmc #endif // OPENMC_STATE_POINT_H diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 446052fbb..0f52649d8 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -638,7 +638,7 @@ void write_eigenvalue_hdf5(hid_t group) write_dataset(group, "k_combined", k_combined); } -extern "C" void read_eigenvalue_hdf5(hid_t group) +void read_eigenvalue_hdf5(hid_t group) { read_dataset(group, "generations_per_batch", settings::gen_per_batch); int n = simulation::restart_batch*settings::gen_per_batch; diff --git a/src/simulation.cpp b/src/simulation.cpp index cf6000398..7eba7f628 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -30,7 +30,6 @@ extern "C" void accumulate_tallies(); extern "C" void allocate_tally_results(); extern "C" void check_triggers(); extern "C" void init_tally_routines(); -extern "C" void load_state_point(); extern "C" void setup_active_tallies(); extern "C" void simulation_init_f(); extern "C" void simulation_finalize_f(); diff --git a/src/state_point.F90 b/src/state_point.F90 index 04e354142..113ce1749 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -13,7 +13,6 @@ module state_point use, intrinsic :: ISO_C_BINDING - use bank_header, only: Bank use constants use endf, only: reaction_name use error, only: fatal_error, warning, write_message @@ -21,11 +20,9 @@ module state_point use message_passing use mgxs_interface use nuclide_header, only: nuclides - use output, only: time_stamp - use random_lcg, only: openmc_get_seed, openmc_set_seed use settings use simulation_header - use string, only: to_str, count_digits, zero_padded, to_f_string + use string, only: to_str, zero_padded use tally_header use tally_filter_header use tally_derivative_header, only: tally_derivs @@ -33,18 +30,6 @@ module state_point implicit none - interface - subroutine write_source_bank(group_id) bind(C) - import HID_T - integer(HID_T), value :: group_id - end subroutine write_source_bank - - subroutine read_source_bank(group_id) bind(C) - import HID_T - integer(HID_T), value :: group_id - end subroutine read_source_bank - end interface - contains !=============================================================================== @@ -248,181 +233,39 @@ contains ! LOAD_STATE_POINT !=============================================================================== - subroutine load_state_point() bind(C) + subroutine load_state_point_f(file_id) bind(C) + integer(HID_T), value :: file_id integer :: i - integer :: int_array(3) - integer, allocatable :: array(:) - integer(C_INT64_T) :: seed - integer(HID_T) :: file_id + integer :: temp integer(HID_T) :: tallies_group integer(HID_T) :: tally_group - logical :: source_present - character(MAX_WORD_LEN) :: word - interface - subroutine read_eigenvalue_hdf5(group) bind(C) - import HID_T - integer(HID_T), value :: group - end subroutine + ! Read global tally data + call read_dataset(global_tallies, file_id, "global_tallies") - subroutine restart_set_keff() bind(C) - end subroutine - end interface + ! Check if tally results are present + call read_attribute(temp, file_id, "tallies_present") - ! Write message - call write_message("Loading state point " // trim(path_state_point) & - // "...", 5) + ! Read in sum and sum squared + if (temp == 1) then + tallies_group = open_group(file_id, "tallies") - ! Open file for reading - file_id = file_open(path_state_point, 'r', parallel=.true.) + TALLY_RESULTS: do i = 1, n_tallies + associate (t => tallies(i) % obj) + ! Read sum, sum_sq, and N for each bin + tally_group = open_group(tallies_group, "tally " // & + trim(to_str(t % id))) + call t % read_results_hdf5(tally_group) + call read_dataset(t % n_realizations, tally_group, & + "n_realizations") + call close_group(tally_group) + end associate + end do TALLY_RESULTS - ! Read filetype - call read_attribute(word, file_id, "filetype") - if (word /= 'statepoint') then - call fatal_error("OpenMC tried to restart from a non-statepoint file.") + call close_group(tallies_group) end if - ! Read revision number for state point file and make sure it matches with - ! current version - call read_attribute(array, file_id, "version") - if (any(array /= VERSION_STATEPOINT)) then - call fatal_error("State point version does not match current version & - &in OpenMC.") - end if - - ! Read and overwrite random number seed - call read_dataset(seed, file_id, "seed") - call openmc_set_seed(seed) - - ! It is not impossible for a state point to be generated from a CE run but - ! to be loaded in to an MG run (or vice versa), check to prevent that. - call read_dataset(word, file_id, "energy_mode") - if (word == "multi-group" .and. run_CE) then - call fatal_error("State point file is from multi-group run but & - & current run is continous-energy!") - else if (word == "continuous-energy" .and. .not. run_CE) then - call fatal_error("State point file is from continuous-energy run but & - & current run is multi-group!") - end if - - ! Read and overwrite run information except number of batches - call read_dataset(word, file_id, "run_mode") - select case(word) - case ('fixed source') - run_mode = MODE_FIXEDSOURCE - case ('eigenvalue') - run_mode = MODE_EIGENVALUE - end select - call read_attribute(int_array(1), file_id, "photon_transport") - if (int_array(1) == 1) then - photon_transport = .true. - else - photon_transport = .false. - end if - call read_dataset(n_particles, file_id, "n_particles") - call read_dataset(int_array(1), file_id, "n_batches") - - ! Take maximum of statepoint n_batches and input n_batches - n_batches = max(n_batches, int_array(1)) - - ! Read batch number to restart at - call read_dataset(restart_batch, file_id, "current_batch") - - ! Check for source in statepoint if needed - call read_attribute(int_array(1), file_id, "source_present") - if (int_array(1) == 1) then - source_present = .true. - else - source_present = .false. - end if - - if (restart_batch > n_batches) then - call fatal_error("The number batches specified in settings.xml is fewer & - & than the number of batches in the given statepoint file.") - end if - - ! Read information specific to eigenvalue run - if (run_mode == MODE_EIGENVALUE) then - call read_dataset(int_array(1), file_id, "n_inactive") - call read_eigenvalue_hdf5(file_id) - - ! Take maximum of statepoint n_inactive and input n_inactive - n_inactive = max(n_inactive, int_array(1)) - end if - - ! Read number of realizations for global tallies - call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.) - - ! Set k_sum, keff, and current_batch based on whether restart file is part - ! of active cycle or inactive cycle - call restart_set_keff() - current_batch = restart_batch - - ! Check to make sure source bank is present - if (path_source_point == path_state_point .and. .not. source_present) then - call fatal_error("Source bank must be contained in statepoint restart & - &file") - end if - - ! Read tallies to master. If we are using Parallel HDF5, all processes - ! need to be included in the HDF5 calls. -#ifdef PHDF5 - if (.true.) then -#else - if (master) then -#endif - ! Read global tally data - call read_dataset(global_tallies, file_id, "global_tallies") - - ! Check if tally results are present - call read_attribute(int_array(1), file_id, "tallies_present") - - ! Read in sum and sum squared - if (int_array(1) == 1) then - tallies_group = open_group(file_id, "tallies") - - TALLY_RESULTS: do i = 1, n_tallies - associate (t => tallies(i) % obj) - ! Read sum, sum_sq, and N for each bin - tally_group = open_group(tallies_group, "tally " // & - trim(to_str(t % id))) - call t % read_results_hdf5(tally_group) - call read_dataset(t % n_realizations, tally_group, & - "n_realizations") - call close_group(tally_group) - end associate - end do TALLY_RESULTS - - call close_group(tallies_group) - end if - end if - - ! Read source if in eigenvalue mode - if (run_mode == MODE_EIGENVALUE) then - - ! Check if source was written out separately - if (.not. source_present) then - - ! Close statepoint file - call file_close(file_id) - - ! Write message - call write_message("Loading source file " // trim(path_source_point) & - // "...", 5) - - ! Open source file - file_id = file_open(path_source_point, 'r', parallel=.true.) - end if - - ! Read source - call read_source_bank(file_id) - - end if - - ! Close file - call file_close(file_id) - - end subroutine load_state_point + end subroutine load_state_point_f end module state_point diff --git a/src/state_point.cpp b/src/state_point.cpp index 3668b7b7c..6dc8eafdd 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -1,6 +1,7 @@ #include "openmc/state_point.h" #include +#include // for int64_t #include // for setfill, setw #include #include @@ -26,6 +27,7 @@ namespace openmc { extern "C" void statepoint_write_f(hid_t file_id); +extern "C" void load_state_point_f(hid_t file_id); extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) @@ -160,6 +162,147 @@ openmc_statepoint_write(const char* filename, bool* write_source) return 0; } +void restart_set_keff() +{ + if (simulation::restart_batch > settings::n_inactive) { + for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) { + simulation::k_sum[0] += simulation::k_generation[i]; + simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2); + } + int n = settings::gen_per_batch*n_realizations; + simulation::keff = simulation::k_sum[0] / n; + } else { + simulation::keff = simulation::k_generation.back(); + } +} + +void load_state_point() +{ + // Write message + write_message("Loading state point " + settings::path_statepoint + "...", 5); + + // Open file for reading + hid_t file_id = file_open(settings::path_statepoint.c_str(), 'r', true); + + // Read filetype + std::string word; + read_attribute(file_id, "filetype", word); + if (word != "statepoint") { + fatal_error("OpenMC tried to restart from a non-statepoint file."); + } + + // Read revision number for state point file and make sure it matches with + // current version + std::array array; + read_attribute(file_id, "version", array); + if (array != VERSION_STATEPOINT) { + fatal_error("State point version does not match current version in OpenMC."); + } + + // Read and overwrite random number seed + int64_t seed; + read_dataset(file_id, "seed", seed); + openmc_set_seed(seed); + + // It is not impossible for a state point to be generated from a CE run but + // to be loaded in to an MG run (or vice versa), check to prevent that. + read_dataset(file_id, "energy_mode", word); + if (word == "multi-group" && settings::run_CE) { + fatal_error("State point file is from multigroup run but current run is " + "continous energy."); + } else if (word == "continuous-energy" && !settings::run_CE) { + fatal_error("State point file is from continuous-energy run but current " + "run is multigroup!"); + } + + // Read and overwrite run information except number of batches + read_dataset(file_id, "run_mode", word); + if (word == "fixed source") { + settings::run_mode = RUN_MODE_FIXEDSOURCE; + } else if (word == "eigenvalue") { + settings::run_mode = RUN_MODE_EIGENVALUE; + } + read_attribute(file_id, "photon_transport", settings::photon_transport); + read_dataset(file_id, "n_particles", settings::n_particles); + int temp; + read_dataset(file_id, "n_batches", temp); + + // Take maximum of statepoint n_batches and input n_batches + settings::n_batches = std::max(settings::n_batches, temp); + + // 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."); + } + + // Read information specific to eigenvalue run + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + read_dataset(file_id, "n_inactive", temp); + read_eigenvalue_hdf5(file_id); + + // Take maximum of statepoint n_inactive and input n_inactive + settings::n_inactive = std::max(settings::n_inactive, temp); + } + + // Read number of realizations for global tallies + read_dataset(file_id, "n_realizations", n_realizations); + + // 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"); + } + + // Read tallies to master. If we are using Parallel HDF5, all processes + // need to be included in the HDF5 calls. +#ifdef PHDF5 + if (true) { +#else + if (mpi::master) { +#endif + // Read tally results + load_state_point_f(file_id); + } + + // Read source if in eigenvalue mode + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + + // Check if source was written out separately + if (!source_present) { + + // Close statepoint file + file_close(file_id); + + // Write message + write_message("Loading source file " + settings::path_sourcepoint + + "...", 5); + + // Open source file + file_id = file_open(settings::path_source.c_str(), 'r', true); + } + + // Read source + read_source_bank(file_id); + + } + + // Close file + file_close(file_id); +} + + hid_t h5banktype() { // Create type for array of 3 reals hsize_t dims[] {3}; @@ -455,18 +598,4 @@ void write_tally_results_nr(hid_t file_id) } } -void restart_set_keff() -{ - if (simulation::restart_batch > settings::n_inactive) { - for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) { - simulation::k_sum[0] += simulation::k_generation[i]; - simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2); - } - int n = settings::gen_per_batch*n_realizations; - simulation::keff = simulation::k_sum[0] / n; - } else { - simulation::keff = simulation::k_generation.back(); - } -} - } // namespace openmc diff --git a/src/string.F90 b/src/string.F90 index 7e1ca41ba..6a2b3ae90 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -247,25 +247,6 @@ contains end function ends_with -!=============================================================================== -! COUNT_DIGITS returns the number of digits needed to represent the input -! integer. -!=============================================================================== - - pure function count_digits(num) result(n_digits) - integer, intent(in) :: num - integer :: n_digits - - n_digits = 1 - do while (num / 10**(n_digits) /= 0 .and. abs(num / 10 **(n_digits-1)) /= 1& - &.and. n_digits /= 10) - ! Note that 10 digits is the maximum needed to represent an integer(4) so - ! the loop automatically exits when n_digits = 10. - n_digits = n_digits + 1 - end do - - end function count_digits - !=============================================================================== ! INT4_TO_STR converts an integer(4) to a string. !===============================================================================