diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 05512b379f..3ead7d8339 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -95,7 +95,7 @@ extern "C" { int openmc_sphharm_filter_get_cosine(int32_t index, char cosine[]); int openmc_sphharm_filter_set_order(int32_t index, int order); int openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]); - int openmc_statepoint_write(const char filename[]); + int openmc_statepoint_write(const char filename[], bool* write_source); int openmc_tally_allocate(int32_t index, const char* type); int openmc_tally_get_active(int32_t index, bool* active); int openmc_tally_get_estimator(int32_t index, int32_t* estimator); @@ -136,7 +136,6 @@ extern "C" { // Global variables extern char openmc_err_msg[256]; extern int32_t n_cells; - extern int32_t n_filters; extern int32_t n_lattices; extern int32_t n_materials; extern int n_nuclides; diff --git a/include/openmc/container_util.h b/include/openmc/container_util.h new file mode 100644 index 0000000000..1e65a18df6 --- /dev/null +++ b/include/openmc/container_util.h @@ -0,0 +1,17 @@ +#ifndef OPENMC_CONTAINER_UTIL_H +#define OPENMC_CONTAINER_UTIL_H + +#include // for find +#include // for begin, end + +namespace openmc { + +template +inline bool contains(const C& v, const T& x) +{ + return std::end(v) != std::find(std::begin(v), std::end(v), x); +} + +} + +#endif // OPENMC_CONTAINER_UTIL_H diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index cbf0baa2d1..b7a75becbe 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -1,6 +1,10 @@ +//! \file eigenvalue.h +//! \brief Data/functions related to k-eigenvalue calculations + #ifndef OPENMC_EIGENVALUE_H #define OPENMC_EIGENVALUE_H +#include #include // for int64_t #include @@ -15,6 +19,7 @@ namespace openmc { //============================================================================== extern double keff_generation; //!< Single-generation k on each processor +extern std::array k_sum; //!< Used to reduce sum and sum_sq extern std::vector entropy; //!< Shannon entropy at each generation extern xt::xtensor source_frac; //!< Source fraction for UFS @@ -28,6 +33,28 @@ extern "C" int64_t n_bank; //! Collect/normalize the tracklength keff from each process extern "C" void calculate_generation_keff(); +//! Calculate mean/standard deviation of keff during active generations +//! +//! This function sets the global variables keff and keff_std which represent +//! the mean and standard deviation of the mean of k-effective over active +//! generations. It also broadcasts the value from the master process. +extern "C" void calculate_average_keff(); + + +//! Calculates a minimum variance estimate of k-effective +//! +//! The minimum variance estimate is based on a linear combination of the +//! collision, absorption, and tracklength estimates. The theory behind this can +//! be found in M. Halperin, "Almost linearly-optimum combination of unbiased +//! estimates," J. Am. Stat. Assoc., 56, 36-43 (1961), +//! doi:10.1080/01621459.1961.10482088. The implementation here follows that +//! described in T. Urbatsch et al., "Estimation and interpretation of keff +//! confidence intervals in MCNP," Nucl. Technol., 111, 169-182 (1995). +//! +//! \param[out] k_combined Estimate of k-effective and its standard deviation +//! \return Error status +extern "C" int openmc_get_keff(double* k_combined); + //! Sample/redistribute source sites from accumulated fission sites extern "C" void synchronize_bank(); diff --git a/include/openmc/finalize.h b/include/openmc/finalize.h index 5a82388414..f99201abe7 100644 --- a/include/openmc/finalize.h +++ b/include/openmc/finalize.h @@ -3,8 +3,6 @@ namespace openmc { -extern "C" void openmc_free_bank(); - } // namespace openmc #endif // OPENMC_FINALIZE_H diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 7fa5370c76..313b841e20 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -52,7 +52,7 @@ void write_string(hid_t group_id, const char* name, const std::string& buffer, std::vector attribute_shape(hid_t obj_id, const char* name); std::vector dataset_names(hid_t group_id); -void ensure_exists(hid_t group_id, const char* name); +void ensure_exists(hid_t obj_id, const char* name, bool attribute=false); std::vector group_names(hid_t group_id); std::vector object_shape(hid_t obj_id); std::string object_name(hid_t obj_id); @@ -202,7 +202,13 @@ read_attribute(hid_t obj_id, const char* name, std::vector& vec) read_attr_string(obj_id, name, n, buffer[0]); for (int i = 0; i < m; ++i) { - vec.emplace_back(&buffer[i][0], std::min(strlen(buffer[i]), n)); + // Determine proper length of string -- strlen doesn't work because + // buffer[i] might not have any null characters + std::size_t k = 0; + for (; k < n; ++k) if (buffer[i][k] == '\0') break; + + // Create string based on (char*, size_t) constructor + vec.emplace_back(&buffer[i][0], k); } } diff --git a/include/openmc/initialize.h b/include/openmc/initialize.h index 72ccc12f68..d2f4ff0c67 100644 --- a/include/openmc/initialize.h +++ b/include/openmc/initialize.h @@ -5,9 +5,6 @@ #include "mpi.h" #endif -extern "C" void print_usage(); -extern "C" void print_version(); - namespace openmc { int parse_command_line(int argc, char* argv[]); diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h index e0addf5bc0..c813755111 100644 --- a/include/openmc/mgxs_interface.h +++ b/include/openmc/mgxs_interface.h @@ -7,6 +7,7 @@ #include "hdf5_interface.h" #include "mgxs.h" +#include namespace openmc { @@ -17,9 +18,9 @@ namespace openmc { extern std::vector nuclides_MG; extern std::vector macro_xs; extern "C" int num_energy_groups; - -//TODO: When more of the Fortran is converted (input_xml, tallies, etc, also -// bring over energy_bin_avg, energy_bins, etc, as vectors) +extern std::vector energy_bins; +extern std::vector energy_bin_avg; +extern std::vector rev_energy_bins; //============================================================================== // Mgxs data loading interface methods @@ -39,6 +40,8 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[], int n_temps, const double temps[], const double atom_densities[], double tolerance, int& method); +extern "C" void read_mg_cross_sections_header_c(hid_t file_id); + //============================================================================== // Mgxs tracking/transport/tallying interface methods //============================================================================== diff --git a/include/openmc/output.h b/include/openmc/output.h index 1ba5d65991..e396c07d8c 100644 --- a/include/openmc/output.h +++ b/include/openmc/output.h @@ -20,7 +20,7 @@ void header(const char* msg, int level); //! Display information regarding cell overlap checking. //============================================================================== -extern "C" void print_overlap_check(); +void print_overlap_check(); extern "C" void title(); diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 298b2878a0..6f2d58832b 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -146,9 +146,7 @@ extern "C" { //! site may have been produced from an external source, from fission, or //! simply as a secondary particle. //! \param src Source site data - //! \param run_CE Whether continuous-energy data is being used - //! \param energy_bin_avg An array of energy group bin averages - void from_source(const Bank* src, bool run_CE, const double* energy_bin_avg); + void from_source(const Bank* src); //! mark a particle as lost and create a particle restart file //! \param message A warning message to display @@ -174,8 +172,7 @@ extern "C" { void particle_create_secondary(Particle* p, const double* uvw, double E, int type, bool run_CE); void particle_initialize(Particle* p); - void particle_from_source(Particle* p, const Bank* src, bool run_CE, - const double* energy_bin_avg); + void particle_from_source(Particle* p, const Bank* src); void particle_mark_as_lost(Particle* p, const char* message); void particle_write_restart(Particle* p); diff --git a/include/openmc/random_lcg.h b/include/openmc/random_lcg.h index b02deec2c5..8fce803924 100644 --- a/include/openmc/random_lcg.h +++ b/include/openmc/random_lcg.h @@ -17,6 +17,7 @@ extern "C" const int STREAM_SOURCE; extern "C" const int STREAM_URR_PTABLE; extern "C" const int STREAM_VOLUME; extern "C" const int STREAM_PHOTON; +constexpr int64_t DEFAULT_SEED = 1; //============================================================================== //! Generate a pseudo-random number using a linear congruential generator. diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 6565d3e460..aa9f9c0565 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -7,6 +7,8 @@ #include #include #include +#include +#include #include "pugixml.hpp" @@ -66,23 +68,25 @@ extern "C" int32_t gen_per_batch; //!< number of generations per batch extern "C" int64_t n_particles; //!< number of particles per generation extern "C" int electron_treatment; //!< how to treat secondary electrons -extern "C" double energy_cutoff[4]; //!< Energy cutoff in [eV] for each particle type +extern "C" std::array energy_cutoff; //!< Energy cutoff in [eV] for each particle type extern "C" int legendre_to_tabular_points; //!< number of points to convert Legendres extern "C" int max_order; //!< Maximum Legendre order for multigroup data extern "C" int n_log_bins; //!< number of bins for logarithmic energy grid extern "C" int n_max_batches; //!< Maximum number of batches - extern "C" int res_scat_method; //!< resonance upscattering method extern "C" double res_scat_energy_min; //!< Min energy in [eV] for res. upscattering extern "C" double res_scat_energy_max; //!< Max energy in [eV] for res. upscattering extern "C" int run_mode; //!< Run mode (eigenvalue, fixed src, etc.) +extern std::unordered_set sourcepoint_batch; //!< Batches when source should be written +extern std::unordered_set statepoint_batch; //!< Batches when state should be written extern "C" int temperature_method; //!< method for choosing temperatures extern "C" double temperature_tolerance; //!< Tolerance in [K] on choosing temperatures extern "C" double temperature_default; //!< Default T in [K] -extern "C" double temperature_range[2]; //!< Min/max T in [K] over which to load xs +extern "C" std::array temperature_range; //!< Min/max T in [K] over which to load xs extern "C" int trace_batch; //!< Batch to trace particle on extern "C" int trace_gen; //!< Generation to trace particle on extern "C" int64_t trace_particle; //!< Particle ID to enable trace on +extern std::vector> track_identifiers; //!< Particle numbers for writing tracks extern "C" int trigger_batch_interval; //!< Batch interval for triggers extern "C" int verbosity; //!< How verbose to make output extern "C" double weight_cutoff; //!< Weight cutoff for Russian roulette diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index e21ea05d50..2002e95fed 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -4,11 +4,17 @@ #ifndef OPENMC_SIMULATION_H #define OPENMC_SIMULATION_H +#include "openmc/particle.h" + #include #include namespace openmc { +constexpr int STATUS_EXIT_NORMAL {0}; +constexpr int STATUS_EXIT_MAX_BATCH {1}; +constexpr int STATUS_EXIT_ON_TRIGGER {2}; + //============================================================================== // Global variable declarations //============================================================================== @@ -18,6 +24,7 @@ namespace simulation { extern "C" int current_batch; //!< current batch extern "C" int current_gen; //!< current fission generation extern "C" int64_t current_work; //!< index in source back of current particle +extern "C" bool initialized; //!< has simulation been initialized? extern "C" double keff; //!< average k over batches extern "C" double keff_std; //!< standard deviation of average k extern "C" double k_col_abs; //!< sum over batches of k_collision * k_absorption @@ -28,7 +35,6 @@ extern "C" int n_lost_particles; //!< cumulative number of lost particles extern "C" bool need_depletion_rx; //!< need to calculate depletion rx? extern "C" int restart_batch; //!< batch at which a restart job resumed extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied? -extern "C" bool simulation_initialized; //!< has simulation been initialized? extern "C" int total_gen; //!< total number of generations simulated extern "C" int64_t work; //!< number of particles per process @@ -53,18 +59,29 @@ extern "C" int thread_id; //!< ID of a given thread //! Determine number of particles to transport per process void calculate_work(); -//! Initialize simulation -extern "C" void openmc_simulation_init_c(); +//! Initialize a batch +void initialize_batch(); //! Initialize a fission generation -extern "C" void initialize_generation(); +void initialize_generation(); + +void initialize_history(Particle* p, int64_t index_source); + +//! Finalize a batch +//! +//! Handles synchronization and accumulation of tallies, calculation of Shannon +//! entropy, getting single-batch estimate of keff, and turning on tallies when +//! appropriate +void finalize_batch(); + +//! Finalize a fission generation +void finalize_generation(); //! Determine overall generation number extern "C" int overall_generation(); #ifdef OPENMC_MPI -extern "C" void broadcast_results(); -extern "C" void broadcast_triggers(); +void broadcast_results(); #endif } // namespace openmc diff --git a/include/openmc/source.h b/include/openmc/source.h index 2e8c9491a2..efc6f12a61 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -58,6 +58,9 @@ extern "C" void initialize_source(); //! \return Sampled source site extern "C" Bank sample_external_source(); +//! Fill source bank at end of generation for fixed source simulations +void fill_source_bank_fixedsource(); + } // namespace openmc #endif // OPENMC_SOURCE_H diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index 058edc09cf..6b019e7b8f 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -9,10 +9,11 @@ namespace openmc { +void write_source_point(const char* filename); extern "C" void write_source_bank(hid_t group_id, Bank* source_bank); extern "C" void read_source_bank(hid_t group_id, Bank* source_bank); - extern "C" void write_tally_results_nr(hid_t file_id); +extern "C" void restart_set_keff(); } // namespace openmc #endif // OPENMC_STATE_POINT_H diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index c323d3bace..05b5d141f9 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -14,9 +14,12 @@ namespace openmc { extern "C" double total_weight; // Threadprivate variables - extern "C" double global_tally_absorption; -#pragma omp threadprivate(global_tally_absorption) +extern "C" double global_tally_collision; +extern "C" double global_tally_tracklength; +extern "C" double global_tally_leakage; +#pragma omp threadprivate(global_tally_absorption, global_tally_collision, \ + global_tally_tracklength, global_tally_leakage) //============================================================================== // Non-member functions diff --git a/include/openmc/timer.h b/include/openmc/timer.h index 11dfd1adfe..1b4d71a5c7 100644 --- a/include/openmc/timer.h +++ b/include/openmc/timer.h @@ -38,9 +38,22 @@ private: // Global variables //============================================================================== +extern Timer time_active; extern Timer time_bank; extern Timer time_bank_sample; extern Timer time_bank_sendrecv; +extern Timer time_finalize; +extern Timer time_inactive; +extern Timer time_initialize; +extern Timer time_tallies; +extern Timer time_total; +extern Timer time_transport; + +//============================================================================== +// Non-member functions +//============================================================================== + +void reset_timers(); } // namespace openmc diff --git a/src/api.F90 b/src/api.F90 index b9ccdca48c..c5ab06f592 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -2,185 +2,20 @@ module openmc_api use, intrinsic :: ISO_C_BINDING - use bank_header, only: openmc_source_bank - use constants, only: K_BOLTZMANN - use eigenvalue, only: k_sum, openmc_get_keff + use constants use error use geometry, only: find_cell - use geometry_header - use hdf5_interface - use material_header - use math - use message_passing - use nuclide_header - use initialize, only: openmc_init_f use particle_header - use plot, only: openmc_plot_geometry - use random_lcg, only: openmc_get_seed, openmc_set_seed - use settings - use simulation_header - use state_point, only: openmc_statepoint_write - use tally_header - use tally_filter_header - use tally_filter - use tally, only: openmc_tally_allocate - use simulation - use string, only: to_f_string - use timer_header - use volume_calc, only: openmc_calculate_volumes - + use string, only: to_str #ifdef DAGMC use dagmc_header, only: free_memory_dagmc #endif implicit none - private - public :: openmc_calculate_volumes - public :: openmc_cell_get_id - public :: openmc_cell_set_id - public :: openmc_energy_filter_get_bins - public :: openmc_energy_filter_set_bins - public :: openmc_extend_filters - public :: openmc_extend_cells - public :: openmc_extend_materials - public :: openmc_extend_tallies - public :: openmc_filter_get_id - public :: openmc_filter_get_type - public :: openmc_filter_set_id - public :: openmc_filter_set_type - public :: openmc_finalize - public :: openmc_find_cell - public :: openmc_get_cell_index - public :: openmc_get_keff - public :: openmc_get_filter_index - public :: openmc_get_filter_next_id - public :: openmc_get_material_index - public :: openmc_get_nuclide_index - public :: openmc_get_seed - public :: openmc_get_tally_index - public :: openmc_get_tally_next_id - public :: openmc_global_tallies - public :: openmc_hard_reset - public :: openmc_init_f - public :: openmc_load_nuclide - public :: openmc_material_add_nuclide - public :: openmc_material_get_id - public :: openmc_material_get_densities - public :: openmc_material_set_density - public :: openmc_material_set_densities - public :: openmc_material_set_id - public :: openmc_next_batch - public :: openmc_nuclide_name - public :: openmc_plot_geometry - public :: openmc_reset - public :: openmc_set_seed - public :: openmc_simulation_finalize - public :: openmc_simulation_init - public :: openmc_source_bank - public :: openmc_tally_allocate - public :: openmc_tally_get_estimator - public :: openmc_tally_get_id - public :: openmc_tally_get_filters - public :: openmc_tally_get_n_realizations - public :: openmc_tally_get_nuclides - public :: openmc_tally_get_scores - public :: openmc_tally_get_type - public :: openmc_tally_results - public :: openmc_tally_set_estimator - public :: openmc_tally_set_filters - public :: openmc_tally_set_id - public :: openmc_tally_set_nuclides - public :: openmc_tally_set_scores - public :: openmc_tally_set_type contains -!=============================================================================== -! OPENMC_FINALIZE frees up memory by deallocating arrays and resetting global -! variables -!=============================================================================== - - function openmc_finalize() result(err) bind(C) - integer(C_INT) :: err - - interface - subroutine openmc_free_bank() bind(C) - end subroutine openmc_free_bank - end interface - - ! Clear results - err = openmc_reset() - - ! Reset global variables - assume_separate = .false. - check_overlaps = .false. - confidence_intervals = .false. - create_fission_neutrons = .true. - electron_treatment = ELECTRON_LED - energy_cutoff(:) = [ZERO, 1000.0_8, ZERO, ZERO] - energy_max(:) = [INFINITY, INFINITY] - energy_min(:) = [ZERO, ZERO] - entropy_on = .false. - gen_per_batch = 1 - index_entropy_mesh = -1 - index_ufs_mesh = -1 - keff = ONE - legendre_to_tabular = .true. - legendre_to_tabular_points = C_NONE - n_batch_interval = 1 - n_lost_particles = 0 - n_particles = -1 - n_source_points = 0 - n_state_points = 0 - n_tallies = 0 - output_summary = .true. - output_tallies = .true. - particle_restart_run = .false. - photon_transport = .false. - pred_batches = .false. - reduce_tallies = .true. - res_scat_on = .false. - res_scat_method = RES_SCAT_ARES - res_scat_energy_min = 0.01_8 - res_scat_energy_max = 1000.0_8 - restart_run = .false. - root_universe = -1 - run_CE = .true. - run_mode = -1 - dagmc = .false. - satisfy_triggers = .false. - call openmc_set_seed(DEFAULT_SEED) - source_latest = .false. - source_separate = .false. - source_write = .true. - survival_biasing = .false. - temperature_default = 293.6_8 - temperature_method = TEMPERATURE_NEAREST - temperature_multipole = .false. - temperature_range = [ZERO, ZERO] - temperature_tolerance = 10.0_8 - total_gen = 0 - trigger_on = .false. - ufs = .false. - urr_ptables_on = .true. - verbosity = 7 - weight_cutoff = 0.25_8 - weight_survive = ONE - write_all_tracks = .false. - write_initial_source = .false. - - ! Deallocate arrays - call free_memory() - - err = 0 -#ifdef OPENMC_MPI - ! Free all MPI types - call openmc_free_bank() -#endif - - end function openmc_finalize - !=============================================================================== ! OPENMC_FIND_CELL determines what cell contains a given point in space !=============================================================================== @@ -216,85 +51,26 @@ contains end function openmc_find_cell -!=============================================================================== -! OPENMC_HARD_RESET reset tallies and timers as well as the pseudorandom -! generator state -!=============================================================================== - - function openmc_hard_reset() result(err) bind(C) - integer(C_INT) :: err - - ! Reset all tallies and timers - err = openmc_reset() - - ! Reset total generations and keff guess - keff = ONE - total_gen = 0 - - ! Reset the random number generator state - call openmc_set_seed(DEFAULT_SEED) - end function openmc_hard_reset - -!=============================================================================== -! OPENMC_RESET resets tallies and timers -!=============================================================================== - - function openmc_reset() result(err) bind(C) - integer(C_INT) :: err - - integer :: i - - if (allocated(tallies)) then - do i = 1, size(tallies) - associate (t => tallies(i) % obj) - t % n_realizations = 0 - if (allocated(t % results)) then - t % results(:, :, :) = ZERO - end if - end associate - end do - end if - - ! Reset global tallies - n_realizations = 0 - if (allocated(global_tallies)) then - global_tallies(:, :) = ZERO - end if - k_col_abs = ZERO - k_col_tra = ZERO - k_abs_tra = ZERO - k_sum(:) = ZERO - - ! Reset timers - call time_total % reset() - call time_total % reset() - call time_initialize % reset() - call time_read_xs % reset() - call time_unionize % reset() - call time_tallies % reset() - call time_inactive % reset() - call time_active % reset() - call time_transport % reset() - call time_finalize % reset() - call reset_timers() - - err = 0 - end function openmc_reset - !=============================================================================== ! FREE_MEMORY deallocates and clears all global allocatable arrays in the ! program !=============================================================================== - subroutine free_memory() + subroutine free_memory() bind(C) + use bank_header use cmfd_header + use geometry_header + use material_header use photon_header use plot_header use sab_header use settings + use simulation_header use surface_header use tally_derivative_header + use tally_filter_header + use tally_header use trigger_header use volume_header diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index bcd7e9f5c6..681d08c5c1 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -29,7 +29,7 @@ contains ! EXECUTE_CMFD runs the CMFD calculation !============================================================================== - subroutine execute_cmfd() + subroutine execute_cmfd() bind(C) use cmfd_data, only: set_up_cmfd use cmfd_solver, only: cmfd_solver_execute @@ -73,7 +73,7 @@ contains ! CMFD_INIT_BATCH handles cmfd options at the start of every batch !============================================================================== - subroutine cmfd_init_batch() + subroutine cmfd_init_batch() bind(C) ! Check to activate CMFD diffusion and possible feedback ! this guarantees that when cmfd begins at least one batch of tallies are @@ -375,7 +375,7 @@ contains ! CMFD_TALLY_INIT !=============================================================================== - subroutine cmfd_tally_init() + subroutine cmfd_tally_init() bind(C) integer :: i if (cmfd_run) then do i = 1, size(cmfd_tallies) diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index 036e193b0b..f4b14b6f2b 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -143,7 +143,7 @@ module cmfd_header logical, public :: cmfd_run_adjoint = .false. ! CMFD run logicals - logical, public :: cmfd_on = .false. + logical(C_BOOL), public, bind(C) :: cmfd_on = .false. ! CMFD display info character(len=25), public :: cmfd_display = 'balance' diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 92c731eae0..bdbbd1c2bc 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -2,241 +2,21 @@ module eigenvalue use, intrinsic :: ISO_C_BINDING - use algorithm, only: binary_search - use constants, only: ZERO - use error, only: fatal_error, warning - use math, only: t_percentile - use message_passing - use random_lcg, only: prn, set_particle_seed, advance_prn_seed - use settings + use bank_header use simulation_header - use string, only: to_str - use tally_header - use timer_header implicit none - real(8) :: k_sum(2) = ZERO ! Used to reduce sum and sum_sq - interface - subroutine calculate_generation_keff() bind(C) - end subroutine + function openmc_get_keff(k_combined) result(err) bind(C) + import C_INT, C_DOUBLE + real(C_DOUBLE), intent(out) :: k_combined(2) + integer(C_INT) :: err + end function end interface contains -!=============================================================================== -! CALCULATE_AVERAGE_KEFF calculates the mean and standard deviation of the mean -! of k-effective during active generations and broadcasts the mean to all -! processors -!=============================================================================== - - subroutine calculate_average_keff() - - integer :: i ! overall generation within simulation - integer :: n ! number of active generations - real(8) :: alpha ! significance level for CI - real(8) :: t_value ! t-value for confidence intervals - - ! Determine overall generation and number of active generations - i = overall_generation() - if (current_batch > n_inactive) then - n = gen_per_batch*n_realizations + current_gen - else - n = 0 - end if - - if (n <= 0) then - ! For inactive generations, use current generation k as estimate for next - ! generation - keff = k_generation(i) - - else - ! Sample mean of keff - k_sum(1) = k_sum(1) + k_generation(i) - k_sum(2) = k_sum(2) + k_generation(i)**2 - - ! Determine mean - keff = k_sum(1) / n - - if (n > 1) then - if (confidence_intervals) then - ! Calculate t-value for confidence intervals - alpha = ONE - CONFIDENCE_LEVEL - t_value = t_percentile(ONE - alpha/TWO, n - 1) - else - t_value = ONE - end if - - ! Standard deviation of the sample mean of k - keff_std = t_value * sqrt((k_sum(2)/n - keff**2) / (n - 1)) - end if - end if - - end subroutine calculate_average_keff - -!=============================================================================== -! OPENMC_GET_KEFF calculates a minimum variance estimate of k-effective based on -! a linear combination of the collision, absorption, and tracklength -! estimates. The theory behind this can be found in M. Halperin, "Almost -! linearly-optimum combination of unbiased estimates," J. Am. Stat. Assoc., 56, -! 36-43 (1961), doi:10.1080/01621459.1961.10482088. The implementation here -! follows that described in T. Urbatsch et al., "Estimation and interpretation -! of keff confidence intervals in MCNP," Nucl. Technol., 111, 169-182 (1995). -!=============================================================================== - - function openmc_get_keff(k_combined) result(err) bind(C) - real(C_DOUBLE), intent(out) :: k_combined(2) - integer(C_INT) :: err - - integer :: l ! loop index - integer :: i, j, k ! indices referring to collision, absorption, or track - real(8) :: n ! number of realizations - real(8) :: kv(3) ! vector of k-effective estimates - real(8) :: cov(3,3) ! sample covariance matrix - real(8) :: f ! weighting factor - real(8) :: g ! sum of weighting factors - real(8) :: S(3) ! sums used for variance calculation - - k_combined = ZERO - - ! Make sure we have at least four realizations. Notice that at the end, - ! there is a N-3 term in a denominator. - if (n_realizations <= 3) then - err = -1 - return - end if - - ! Initialize variables - n = real(n_realizations, 8) - - ! Copy estimates of k-effective and its variance (not variance of the mean) - kv(1) = global_tallies(RESULT_SUM, K_COLLISION) / n - kv(2) = global_tallies(RESULT_SUM, K_ABSORPTION) / n - kv(3) = global_tallies(RESULT_SUM, K_TRACKLENGTH) / n - cov(1, 1) = (global_tallies(RESULT_SUM_SQ, K_COLLISION) - & - n * kv(1) * kv(1)) / (n - ONE) - cov(2, 2) = (global_tallies(RESULT_SUM_SQ, K_ABSORPTION) - & - n * kv(2) * kv(2)) / (n - ONE) - cov(3, 3) = (global_tallies(RESULT_SUM_SQ, K_TRACKLENGTH) - & - n * kv(3) * kv(3)) / (n - ONE) - - ! Calculate covariances based on sums with Bessel's correction - cov(1, 2) = (k_col_abs - n * kv(1) * kv(2)) / (n - ONE) - cov(1, 3) = (k_col_tra - n * kv(1) * kv(3)) / (n - ONE) - cov(2, 3) = (k_abs_tra - n * kv(2) * kv(3)) / (n - ONE) - cov(2, 1) = cov(1, 2) - cov(3, 1) = cov(1, 3) - cov(3, 2) = cov(2, 3) - - ! Check to see if two estimators are the same; this is guaranteed to happen - ! in MG-mode with survival biasing when the collision and absorption - ! estimators are the same, but can theoretically happen at anytime. - ! If it does, the standard estimators will produce floating-point - ! exceptions and an expression specifically derived for the combination of - ! two estimators (vice three) should be used instead. - - ! First we will identify if there are any matching estimators - if ((abs(kv(1) - kv(2)) / kv(1) < FP_REL_PRECISION) .and. & - (abs(cov(1, 1) - cov(2, 2)) / cov(1, 1) < FP_REL_PRECISION)) then - ! 1 and 2 match, so only use 1 and 3 in our comparisons - i = 1 - j = 3 - - else if ((abs(kv(1) - kv(3)) / kv(1) < FP_REL_PRECISION) .and. & - (abs(cov(1, 1) - cov(3, 3)) / cov(1, 1) < FP_REL_PRECISION)) then - ! 1 and 3 match, so only use 1 and 2 in our comparisons - i = 1 - j = 2 - - else if ((abs(kv(2) - kv(3)) / kv(2) < FP_REL_PRECISION) .and. & - (abs(cov(2, 2) - cov(3, 3)) / cov(2, 2) < FP_REL_PRECISION)) then - ! 2 and 3 match, so only use 1 and 2 in our comparisons - i = 1 - j = 2 - - else - ! No two estimators match, so set i to 0 and this will be the indicator - ! to use all three estimators. - i = 0 - end if - - if (i == 0) then - ! Use three estimators as derived in the paper by Urbatsch - - ! Initialize variables - g = ZERO - S = ZERO - - do l = 1, 3 - ! Permutations of estimates - if (l == 1) then - ! i = collision, j = absorption, k = tracklength - i = 1 - j = 2 - k = 3 - elseif (l == 2) then - ! i = absortion, j = tracklength, k = collision - i = 2 - j = 3 - k = 1 - elseif (l == 3) then - ! i = tracklength, j = collision, k = absorption - i = 3 - j = 1 - k = 2 - end if - - ! Calculate weighting - f = cov(j, j) * (cov(k, k) - cov(i, k)) - cov(k, k) * cov(i, j) + & - cov(j, k) * (cov(i, j) + cov(i, k) - cov(j, k)) - - ! Add to S sums for variance of combined estimate - S(1) = S(1) + f * cov(1, l) - S(2) = S(2) + (cov(j, j) + cov(k, k) - TWO * cov(j, k)) * kv(l) * kv(l) - S(3) = S(3) + (cov(k, k) + cov(i, j) - cov(j, k) - & - cov(i, k)) * kv(l) * kv(j) - - ! Add to sum for combined k-effective - k_combined(1) = k_combined(1) + f * kv(l) - g = g + f - end do - - ! Complete calculations of S sums - S = (n - ONE) * S - S(1) = (n - ONE)**2 * S(1) - - ! Calculate combined estimate of k-effective - k_combined(1) = k_combined(1) / g - - ! Calculate standard deviation of combined estimate - g = (n - ONE)**2 * g - k_combined(2) = sqrt(S(1) / & - (g * n * (n - THREE)) * (ONE + n * ((S(2) - TWO * S(3)) / g))) - - else - ! Use only two estimators - ! These equations are derived analogously to that done in the paper by - ! Urbatsch, but are simpler than for the three estimators case since the - ! block matrices of the three estimator equations reduces to scalars here - - ! Store the commonly used term - f = kv(i) - kv(j) - g = cov(i, i) + cov(j, j) - TWO * cov(i, j) - - ! Calculate combined estimate of k-effective - k_combined(1) = kv(i) - (cov(i, i) - cov(i, j)) / g * f - - ! Calculate standard deviation of combined estimate - k_combined(2) = (cov(i, i) * cov(j, j) - cov(i, j) * cov(i, j)) * & - (g + n * f * f) / (n * (n - TWO) * g * g) - k_combined(2) = sqrt(k_combined(2)) - - end if - err = 0 - - end function openmc_get_keff - #ifdef _OPENMP !=============================================================================== ! JOIN_BANK_FROM_THREADS joins threadprivate fission banks into a single fission @@ -244,7 +24,7 @@ contains ! to preserve the order of the bank when using varying numbers of threads. !=============================================================================== - subroutine join_bank_from_threads() + subroutine join_bank_from_threads() bind(C) integer(8) :: total ! total number of fission bank sites integer :: i ! loop index for threads diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index a9f99a946e..927b014b52 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -1,5 +1,6 @@ #include "openmc/eigenvalue.h" +#include "xtensor/xbuilder.hpp" #include "xtensor/xmath.hpp" #include "xtensor/xtensor.hpp" #include "xtensor/xview.hpp" @@ -8,6 +9,7 @@ #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/hdf5_interface.h" +#include "openmc/math_functions.h" #include "openmc/mesh.h" #include "openmc/message_passing.h" #include "openmc/random_lcg.h" @@ -18,6 +20,8 @@ #include "openmc/tallies/tally.h" #include // for min +#include +#include // for sqrt, abs, pow #include namespace openmc { @@ -27,6 +31,7 @@ namespace openmc { //============================================================================== double keff_generation; +std::array k_sum; std::vector entropy; xt::xtensor source_frac; @@ -129,7 +134,7 @@ void synchronize_bank() // Allocate temporary source bank int64_t index_temp = 0; - Bank temp_sites[3*simulation::work]; + std::vector temp_sites(3*simulation::work); for (int64_t i = 0; i < n_bank; ++i) { // If there are less than n_particles particles banked, automatically add @@ -291,13 +296,200 @@ void synchronize_bank() MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE); #else - std::copy(temp_sites, temp_sites + settings::n_particles, source_bank); + std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles, source_bank); #endif time_bank_sendrecv.stop(); time_bank.stop(); } +void calculate_average_keff() +{ + // Determine overall generation and number of active generations + int i = overall_generation() - 1; + int n; + if (simulation::current_batch > settings::n_inactive) { + n = settings::gen_per_batch*n_realizations + simulation::current_gen; + } else { + n = 0; + } + + if (n <= 0) { + // For inactive generations, use current generation k as estimate for next + // generation + simulation::keff = simulation::k_generation[i]; + } else { + // Sample mean of keff + k_sum[0] += simulation::k_generation[i]; + k_sum[1] += std::pow(simulation::k_generation[i], 2); + + // Determine mean + simulation::keff = k_sum[0] / n; + + if (n > 1) { + double t_value; + if (settings::confidence_intervals) { + // Calculate t-value for confidence intervals + double alpha = 1.0 - CONFIDENCE_LEVEL; + t_value = t_percentile_c(1.0 - alpha/2.0, n - 1); + } else { + t_value = 1.0; + } + + // Standard deviation of the sample mean of k + simulation::keff_std = t_value * std::sqrt((k_sum[1]/n - + std::pow(simulation::keff, 2)) / (n - 1)); + } + } +} + +int openmc_get_keff(double* k_combined) +{ + k_combined[0] = 0.0; + k_combined[1] = 0.0; + + // Make sure we have at least four realizations. Notice that at the end, + // there is a N-3 term in a denominator. + if (n_realizations <= 3) { + return -1; + } + + // Initialize variables + int64_t n = n_realizations; + + // Copy estimates of k-effective and its variance (not variance of the mean) + auto gt = global_tallies(); + + std::array kv {}; + xt::xtensor cov = xt::zeros({3, 3}); + kv[0] = gt(K_COLLISION, RESULT_SUM) / n; + kv[1] = gt(K_ABSORPTION, RESULT_SUM) / n; + kv[2] = gt(K_TRACKLENGTH, RESULT_SUM) / n; + cov(0, 0) = (gt(K_COLLISION, RESULT_SUM_SQ) - n*kv[0]*kv[0]) / (n - 1); + cov(1, 1) = (gt(K_ABSORPTION, RESULT_SUM_SQ) - n*kv[1]*kv[1]) / (n - 1); + cov(2, 2) = (gt(K_TRACKLENGTH, RESULT_SUM_SQ) - n*kv[2]*kv[2]) / (n - 1); + + // Calculate covariances based on sums with Bessel's correction + cov(0, 1) = (simulation::k_col_abs - n * kv[0] * kv[1]) / (n - 1); + cov(0, 2) = (simulation::k_col_tra - n * kv[0] * kv[2]) / (n - 1); + cov(1, 2) = (simulation::k_abs_tra - n * kv[1] * kv[2]) / (n - 1); + cov(1, 0) = cov(0, 1); + cov(2, 0) = cov(0, 2); + cov(2, 1) = cov(1, 2); + + // Check to see if two estimators are the same; this is guaranteed to happen + // in MG-mode with survival biasing when the collision and absorption + // estimators are the same, but can theoretically happen at anytime. + // If it does, the standard estimators will produce floating-point + // exceptions and an expression specifically derived for the combination of + // two estimators (vice three) should be used instead. + + // First we will identify if there are any matching estimators + int i, j, k; + if ((std::abs(kv[0] - kv[1]) / kv[0] < FP_REL_PRECISION) && + (std::abs(cov(0, 0) - cov(1, 1)) / cov(0, 0) < FP_REL_PRECISION)) { + // 0 and 1 match, so only use 0 and 2 in our comparisons + i = 0; + j = 2; + + } else if ((std::abs(kv[0] - kv[2]) / kv[0] < FP_REL_PRECISION) && + (std::abs(cov(0, 0) - cov(2, 2)) / cov(0, 0) < FP_REL_PRECISION)) { + // 0 and 2 match, so only use 0 and 1 in our comparisons + i = 0; + j = 1; + + } else if ((std::abs(kv[1] - kv[2]) / kv[1] < FP_REL_PRECISION) && + (std::abs(cov(1, 1) - cov(2, 2)) / cov(1, 1) < FP_REL_PRECISION)) { + // 1 and 2 match, so only use 0 and 1 in our comparisons + i = 0; + j = 1; + + } else { + // No two estimators match, so set i to -1 and this will be the indicator + // to use all three estimators. + i = -1; + } + + if (i == -1) { + // Use three estimators as derived in the paper by Urbatsch + + // Initialize variables + double g = 0.0; + std::array S {}; + + for (int l = 0; l < 3; ++l) { + // Permutations of estimates + switch (l) { + case 0: + // i = collision, j = absorption, k = tracklength + i = 0; + j = 1; + k = 2; + break; + case 1: + // i = absortion, j = tracklength, k = collision + i = 1; + j = 2; + k = 0; + break; + case 2: + // i = tracklength, j = collision, k = absorption + i = 2; + j = 0; + k = 1; + break; + } + + // Calculate weighting + double f = cov(j, j) * (cov(k, k) - cov(i, k)) - cov(k, k) * cov(i, j) + + cov(j, k) * (cov(i, j) + cov(i, k) - cov(j, k)); + + // Add to S sums for variance of combined estimate + S[0] += f * cov(0, l); + S[1] += (cov(j, j) + cov(k, k) - 2.0 * cov(j, k)) * kv[l] * kv[l]; + S[2] += (cov(k, k) + cov(i, j) - cov(j, k) - cov(i, k)) * kv[l] * kv[j]; + + // Add to sum for combined k-effective + k_combined[0] += f * kv[l]; + g += f; + } + + // Complete calculations of S sums + for (auto& S_i : S) { + S_i *= (n - 1); + } + S[0] *= (n - 1)*(n - 1); + + // Calculate combined estimate of k-effective + k_combined[0] /= g; + + // Calculate standard deviation of combined estimate + g *= (n - 1)*(n - 1); + k_combined[1] = std::sqrt(S[0] / (g*n*(n - 3)) * + (1 + n*((S[1] - 2*S[2]) / g))); + + } else { + // Use only two estimators + // These equations are derived analogously to that done in the paper by + // Urbatsch, but are simpler than for the three estimators case since the + // block matrices of the three estimator equations reduces to scalars here + + // Store the commonly used term + double f = kv[i] - kv[j]; + double g = cov(i, i) + cov(j, j) - 2.0*cov(i, j); + + // Calculate combined estimate of k-effective + k_combined[0] = kv[i] - (cov(i, i) - cov(i, j)) / g * f; + + // Calculate standard deviation of combined estimate + k_combined[1] = (cov(i, i)*cov(j, j) - cov(i, j)*cov(i, j)) * + (g + n*f*f) / (n*(n - 2)*g*g); + k_combined[1] = std::sqrt(k_combined[1]); + + } + return 0; +} + void shannon_entropy() { // Get pointer to entropy mesh @@ -430,6 +622,10 @@ extern "C" void read_eigenvalue_hdf5(hid_t group) read_dataset(group, "k_abs_tra", simulation::k_abs_tra); } +//============================================================================== +// Fortran compatibility +//============================================================================== + extern "C" double entropy_c(int i) { return entropy.at(i - 1); @@ -440,5 +636,6 @@ extern "C" void entropy_clear() entropy.clear(); } +extern "C" void k_sum_reset() { k_sum.fill(0.0); } } // namespace openmc diff --git a/src/finalize.cpp b/src/finalize.cpp index 7f07ddac03..94e3d1a7f2 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -1,14 +1,137 @@ #include "openmc/finalize.h" +#include "openmc/capi.h" +#include "openmc/constants.h" +#include "openmc/eigenvalue.h" +#include "openmc/geometry.h" #include "openmc/message_passing.h" +#include "openmc/nuclide.h" +#include "openmc/random_lcg.h" +#include "openmc/settings.h" +#include "openmc/simulation.h" +#include "openmc/timer.h" +#include "openmc/tallies/tally.h" -namespace openmc { +using namespace openmc; -void openmc_free_bank() +// Functions defined in Fortran +extern "C" void free_memory(); +extern "C" void reset_timers_f(); + +int openmc_finalize() { + // Clear results + openmc_reset(); + + // Reset global variables + settings::assume_separate = false; + settings::check_overlaps = false; + settings::confidence_intervals = false; + settings::create_fission_neutrons = true; + settings::electron_treatment = ELECTRON_LED; + settings::energy_cutoff = {0.0, 1000.0, 0.0, 0.0}; + settings::entropy_on = false; + settings::gen_per_batch = 1; + settings::index_entropy_mesh = -1; + settings::index_ufs_mesh = -1; + settings::legendre_to_tabular = true; + settings::legendre_to_tabular_points = -1; + settings::n_particles = -1; + settings::output_summary = true; + settings::output_tallies = true; + settings::particle_restart_run = false; + settings::photon_transport = false; + settings::reduce_tallies = true; + settings::res_scat_on = false; + settings::res_scat_method = RES_SCAT_ARES; + settings::res_scat_energy_min = 0.01; + settings::res_scat_energy_max = 1000.0; + settings::restart_run = false; + settings::run_CE = true; + settings::run_mode = -1; + settings::dagmc = false; + settings::source_latest = false; + settings::source_separate = false; + settings::source_write = true; + settings::survival_biasing = false; + settings::temperature_default = 293.6; + settings::temperature_method = TEMPERATURE_NEAREST; + settings::temperature_multipole = false; + settings::temperature_range = {0.0, 0.0}; + settings::temperature_tolerance = 10.0; + settings::trigger_on = false; + settings::trigger_predict = false; + settings::trigger_batch_interval = 1; + settings::ufs_on = false; + settings::urr_ptables_on = true; + settings::verbosity = 7; + settings::weight_cutoff = 0.25; + settings::weight_survive = 1.0; + settings::write_all_tracks = false; + settings::write_initial_source = false; + + simulation::keff = 1.0; + simulation::n_lost_particles = 0; + simulation::satisfy_triggers = false; + simulation::total_gen = 0; + + energy_max = {INFTY, INFTY}; + energy_min = {0.0, 0.0}; + n_tallies = 0; + openmc_root_universe = -1; + openmc_set_seed(DEFAULT_SEED); + + // Deallocate arrays + free_memory(); + + // Free all MPI types #ifdef OPENMC_MPI MPI_Type_free(&mpi::bank); #endif + + return 0; } -} // namespace openmc +int openmc_reset() +{ + for (int i = 1; i <= n_tallies; ++i) { + openmc_tally_reset(i); + } + + // Reset global tallies (can't really use global_tallies() right now because + // it doesn't have any information about whether the underlying buffer was + // allocated) + n_realizations = 0; + double* buffer = nullptr; + openmc_global_tallies(&buffer); + if (buffer) { + for (int i = 0; i < 3*N_GLOBAL_TALLIES; ++i) { + buffer[i] = 0.0; + } + } + + simulation::k_col_abs = 0.0; + simulation::k_col_tra = 0.0; + simulation::k_abs_tra = 0.0; + k_sum = {0.0, 0.0}; + + // Reset timers + reset_timers(); + reset_timers_f(); + + return 0; +} + +int openmc_hard_reset() +{ + // Reset all tallies and timers + openmc_reset(); + + // Reset total generations and keff guess + simulation::keff = 1.0; + simulation::total_gen = 0; + + // Reset the random number generator state + openmc_set_seed(DEFAULT_SEED); + return 0; +} diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index d72f7d424d..4b0aedd0ae 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -148,13 +148,22 @@ dataset_typesize(hid_t dset) void -ensure_exists(hid_t group_id, const char* name) +ensure_exists(hid_t obj_id, const char* name, bool attribute) { - if (!object_exists(group_id, name)) { - std::stringstream err_msg; - err_msg << "Object \"" << name << "\" does not exist in group " - << object_name(group_id); - fatal_error(err_msg); + if (attribute) { + if (!attribute_exists(obj_id, name)) { + std::stringstream err_msg; + err_msg << "Attribute \"" << name << "\" does not exist in object " + << object_name(obj_id); + fatal_error(err_msg); + } + } else { + if (!object_exists(obj_id, name)) { + std::stringstream err_msg; + err_msg << "Object \"" << name << "\" does not exist in object " + << object_name(obj_id); + fatal_error(err_msg); + } } } diff --git a/src/initialize.F90 b/src/initialize.F90 index 3d805c97d5..4801a90648 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -2,18 +2,8 @@ module initialize use, intrinsic :: ISO_C_BINDING -#ifdef _OPENMP - use omp_lib -#endif - - use bank_header, only: Bank - use constants - use input_xml, only: read_input_xml - use message_passing - use random_lcg, only: openmc_set_seed use settings - use string, only: ends_with, to_f_string - use timer_header + use string, only: to_f_string implicit none @@ -42,56 +32,11 @@ module initialize contains -!=============================================================================== -! OPENMC_INIT takes care of all initialization tasks, i.e. reading -! from command line, reading xml input files, initializing random -! number seeds, reading cross sections, initializing starting source, -! setting up timers, etc. -!=============================================================================== - - function openmc_init_f() result(err) bind(C) - integer(C_INT) :: err - -#ifdef _OPENMP - character(MAX_WORD_LEN) :: envvar -#endif - - ! Start total and initialization timer - call time_total%start() - call time_initialize%start() - -#ifdef _OPENMP - ! Change schedule of main parallel-do loop if OMP_SCHEDULE is set - call get_environment_variable("OMP_SCHEDULE", envvar) - if (len_trim(envvar) == 0) then - call omp_set_schedule(omp_sched_static, 0) - end if -#endif - - ! Read command line arguments - call read_command_line() - - ! Initialize random number generator -- if the user specifies a seed, it - ! will be re-initialized later - call openmc_set_seed(DEFAULT_SEED) - - ! Read XML input files - call read_input_xml() - - ! Check for particle restart run - if (particle_restart_run) run_mode = MODE_PARTICLE - - ! Stop initialization timer - call time_initialize%stop() - - err = 0 - end function openmc_init_f - !=============================================================================== ! READ_COMMAND_LINE reads all parameters from the command line !=============================================================================== - subroutine read_command_line() + subroutine read_command_line() bind(C) ! Arguments were already read on C++ side (initialize.cpp). Here we just ! convert the C-style strings to Fortran style diff --git a/src/initialize.cpp b/src/initialize.cpp index 26d21bf09b..e04ac2922d 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -1,12 +1,13 @@ #include "openmc/initialize.h" #include +#include // for getenv #include #include #include #ifdef _OPENMP -#include "omp.h" +#include #endif #include "openmc/capi.h" @@ -14,14 +15,17 @@ #include "openmc/error.h" #include "openmc/hdf5_interface.h" #include "openmc/message_passing.h" +#include "openmc/random_lcg.h" #include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/string_utils.h" +#include "openmc/timer.h" // data/functions from Fortran side -extern "C" void openmc_init_f(); extern "C" void print_usage(); extern "C" void print_version(); +extern "C" void read_command_line(); +extern "C" void read_input_xml(); // Paths to various files extern "C" { @@ -30,6 +34,8 @@ extern "C" { int openmc_init(int argc, char* argv[], const void* intracomm) { + using namespace openmc; + #ifdef OPENMC_MPI // Check if intracomm was passed MPI_Comm comm; @@ -40,15 +46,40 @@ int openmc_init(int argc, char* argv[], const void* intracomm) } // Initialize MPI for C++ - openmc::initialize_mpi(comm); + initialize_mpi(comm); #endif // Parse command-line arguments - int err = openmc::parse_command_line(argc, argv); + int err = parse_command_line(argc, argv); if (err) return err; - // Continue with rest of initialization - openmc_init_f(); + // Start total and initialization timer + time_total.start(); + time_initialize.start(); + +#ifdef _OPENMP + // If OMP_SCHEDULE is not set, default to a static schedule + char* envvar = std::getenv("OMP_SCHEDULE"); + if (!envvar) { + omp_set_schedule(omp_sched_static, 0); + } +#endif + + // Read command line arguments + read_command_line(); + + // Initialize random number generator -- if the user specifies a seed, it + // will be re-initialized later + openmc_set_seed(DEFAULT_SEED); + + // Read XML input files + read_input_xml(); + + // Check for particle restart run + if (settings::particle_restart_run) settings::run_mode = RUN_MODE_PARTICLE; + + // Stop initialization timer + time_initialize.stop(); return 0; } diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 81a38b16cd..466bff9e7e 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -28,7 +28,7 @@ module input_xml use plot_header use random_lcg, only: prn use surface_header - use set_header, only: SetChar + use set_header, only: SetChar, SetInt use settings use stl_vector, only: VectorInt, VectorReal, VectorChar use string, only: to_lower, to_str, str_to_int, str_to_real, & @@ -112,7 +112,7 @@ contains ! geometry, materials, and tallies. !=============================================================================== - subroutine read_input_xml() + subroutine read_input_xml() bind(C) type(VectorReal), allocatable :: nuc_temps(:) ! List of T to read for each nuclide type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b) @@ -207,10 +207,7 @@ contains integer :: i integer :: n - integer, allocatable :: temp_int_array(:) - integer :: n_tracks type(XMLNode) :: root - type(XMLNode) :: node_sp type(XMLNode) :: node_res_scat type(XMLNode) :: node_vol type(XMLNode), allocatable :: node_vol_list(:) @@ -218,113 +215,6 @@ contains ! Get proper XMLNode type given pointer root % ptr = root_ptr - if (run_mode == MODE_EIGENVALUE) then - ! Preallocate space for keff and entropy by generation - call k_generation_reserve(n_max_batches*gen_per_batch) - end if - - ! Particle tracks - if (check_for_node(root, "track")) then - ! Make sure that there are three values per particle - n_tracks = node_word_count(root, "track") - if (mod(n_tracks, 3) /= 0) then - call fatal_error("Number of integers specified in 'track' is not & - &divisible by 3. Please provide 3 integers per particle to be & - &tracked.") - end if - - ! Allocate space and get list of tracks - allocate(temp_int_array(n_tracks)) - call get_node_array(root, "track", temp_int_array) - - ! Reshape into track_identifiers - allocate(track_identifiers(3, n_tracks/3)) - track_identifiers = reshape(temp_int_array, [3, n_tracks/3]) - end if - - ! Check if the user has specified to write state points - if (check_for_node(root, "state_point")) then - - ! Get pointer to state_point node - node_sp = root % child("state_point") - - ! Determine number of batches at which to store state points - if (check_for_node(node_sp, "batches")) then - n_state_points = node_word_count(node_sp, "batches") - else - n_state_points = 0 - end if - - if (n_state_points > 0) then - ! User gave specific batches to write state points - allocate(temp_int_array(n_state_points)) - call get_node_array(node_sp, "batches", temp_int_array) - do i = 1, n_state_points - call statepoint_batch % add(temp_int_array(i)) - end do - deallocate(temp_int_array) - else - ! If neither were specified, write state point at last batch - n_state_points = 1 - call statepoint_batch % add(n_batches) - end if - else - ! If no tag was present, by default write state point at - ! last batch only - n_state_points = 1 - call statepoint_batch % add(n_batches) - end if - - ! Check if the user has specified to write source points - if (check_for_node(root, "source_point")) then - - ! Get pointer to source_point node - node_sp = root % child("source_point") - - ! Determine number of batches at which to store source points - if (check_for_node(node_sp, "batches")) then - n_source_points = node_word_count(node_sp, "batches") - else - n_source_points = 0 - end if - - if (n_source_points > 0) then - ! User gave specific batches to write source points - allocate(temp_int_array(n_source_points)) - call get_node_array(node_sp, "batches", temp_int_array) - do i = 1, n_source_points - call sourcepoint_batch % add(temp_int_array(i)) - end do - deallocate(temp_int_array) - else - ! If neither were specified, write source points with state points - n_source_points = n_state_points - do i = 1, n_state_points - call sourcepoint_batch % add(statepoint_batch % get_item(i)) - end do - end if - else - ! If no tag was present, by default we keep source bank in - ! statepoint file and write it out at statepoints intervals - n_source_points = n_state_points - do i = 1, n_state_points - call sourcepoint_batch % add(statepoint_batch % get_item(i)) - end do - end if - - ! If source is not seperate and is to be written out in the statepoint file, - ! make sure that the sourcepoint batch numbers are contained in the - ! statepoint list - if (.not. source_separate) then - do i = 1, n_source_points - if (.not. statepoint_batch % contains(sourcepoint_batch % & - get_item(i))) then - call fatal_error('Sourcepoint batches are not a subset& - & of statepoint batches.') - end if - end do - end if - ! Resonance scattering parameters if (check_for_node(root, "resonance_scattering")) then node_res_scat = root % child("resonance_scattering") @@ -1322,7 +1212,11 @@ contains ! Read derivative attributes. do i = 1, size(node_deriv_list) +!$omp parallel +!$omp critical (ReadTallyDeriv) call tally_derivs(i) % from_xml(node_deriv_list(i)) +!$omp end critical (ReadTallyDeriv) +!$omp end parallel ! Update tally derivative dictionary call tally_deriv_dict % set(tally_derivs(i) % id, i) @@ -2714,6 +2608,13 @@ contains integer(HID_T) :: file_id character(len=MAX_WORD_LEN), allocatable :: names(:) + interface + subroutine read_mg_cross_sections_header_c(file_id) bind(C) + import HID_T + integer(HID_T), value :: file_id + end subroutine + end interface + ! Check if MGXS Library exists inquire(FILE=path_cross_sections, EXIST=file_exists) if (.not. file_exists) then @@ -2761,6 +2662,9 @@ contains energy_bin_avg(i) = HALF * (energy_bins(i) + energy_bins(i + 1)) end do + ! Set up energy bins on C++ side + call read_mg_cross_sections_header_c(file_id) + ! Get the minimum and maximum energies energy_min(NEUTRON) = energy_bins(num_energy_groups + 1) energy_max(NEUTRON) = energy_bins(1) diff --git a/src/mgxs.cpp b/src/mgxs.cpp index 12c73fab4f..1b6e85ffb2 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -26,7 +26,6 @@ namespace openmc { std::vector nuclides_MG; std::vector macro_xs; - //============================================================================== // Mgxs base-class methods //============================================================================== diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index c592d28a28..5cec8811f1 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -8,6 +8,14 @@ namespace openmc { +//============================================================================== +// Global variable definitions +//============================================================================== + +std::vector energy_bins; +std::vector energy_bin_avg; +std::vector rev_energy_bins; + //============================================================================== // Mgxs data loading interface methods //============================================================================== @@ -83,6 +91,26 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[], } } +//============================================================================== + +void read_mg_cross_sections_header_c(hid_t file_id) +{ + ensure_exists(file_id, "energy_groups", true); + read_attribute(file_id, "energy_groups", num_energy_groups); + + ensure_exists(file_id, "group structure", true); + read_attribute(file_id, "group structure", energy_bins); + + // Create reverse energy bins + std::copy(energy_bins.crbegin(), energy_bins.crend(), + std::back_inserter(rev_energy_bins)); + + // Create average energies + for (int i = 0; i < energy_bins.size() - 1; ++i) { + energy_bin_avg.push_back(0.5*(energy_bins[i] + energy_bins[i+1])); + } +} + //============================================================================== // Mgxs tracking/transport/tallying interface methods //============================================================================== diff --git a/src/output.F90 b/src/output.F90 index b1a2f4e77d..ef3d3259cb 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -285,7 +285,7 @@ contains ! below them !=============================================================================== - subroutine print_columns() + subroutine print_columns() bind(C) write(UNIT=ou, FMT='(2X,A9,3X)', ADVANCE='NO') "Bat./Gen." write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') " k " @@ -323,7 +323,7 @@ contains ! PRINT_GENERATION displays information for a generation of neutrons. !=============================================================================== - subroutine print_generation() + subroutine print_generation() bind(C) integer :: i ! overall generation integer :: n ! number of active generations @@ -360,7 +360,7 @@ contains ! multiplication factor as well as the average value if we're in active batches !=============================================================================== - subroutine print_batch_keff() + subroutine print_batch_keff() bind(C) integer :: i ! overall generation integer :: n ! number of active generations @@ -490,7 +490,7 @@ contains ! initialization, for computation, and for intergeneration synchronization. !=============================================================================== - subroutine print_runtime() + subroutine print_runtime() bind(C) integer :: n_active real(8) :: speed_inactive ! # of neutrons/second in inactive batches @@ -501,49 +501,49 @@ contains call header("Timing Statistics", 6) ! display time elapsed for various sections - write(ou,100) "Total time for initialization", time_initialize % elapsed + write(ou,100) "Total time for initialization", time_initialize_elapsed() write(ou,100) " Reading cross sections", time_read_xs % elapsed - write(ou,100) "Total time in simulation", time_inactive % elapsed + & - time_active % elapsed - write(ou,100) " Time in transport only", time_transport % elapsed + write(ou,100) "Total time in simulation", time_inactive_elapsed() + & + time_active_elapsed() + write(ou,100) " Time in transport only", time_transport_elapsed() if (run_mode == MODE_EIGENVALUE) then - write(ou,100) " Time in inactive batches", time_inactive % elapsed + write(ou,100) " Time in inactive batches", time_inactive_elapsed() end if - write(ou,100) " Time in active batches", time_active % elapsed + write(ou,100) " Time in active batches", time_active_elapsed() if (run_mode == MODE_EIGENVALUE) then write(ou,100) " Time synchronizing fission bank", time_bank_elapsed() write(ou,100) " Sampling source sites", time_bank_sample_elapsed() write(ou,100) " SEND/RECV source sites", time_bank_sendrecv_elapsed() end if - write(ou,100) " Time accumulating tallies", time_tallies % elapsed + write(ou,100) " Time accumulating tallies", time_tallies_elapsed() if (cmfd_run) write(ou,100) " Time in CMFD", time_cmfd % elapsed if (cmfd_run) write(ou,100) " Building matrices", & time_cmfdbuild % elapsed if (cmfd_run) write(ou,100) " Solving matrices", & time_cmfdsolve % elapsed - write(ou,100) "Total time for finalization", time_finalize % elapsed - write(ou,100) "Total time elapsed", time_total % elapsed + write(ou,100) "Total time for finalization", time_finalize_elapsed() + write(ou,100) "Total time elapsed", time_total_elapsed() ! Calculate particle rate in active/inactive batches n_active = current_batch - n_inactive if (restart_run) then if (restart_batch < n_inactive) then speed_inactive = real(n_particles * (n_inactive - restart_batch) * & - gen_per_batch) / time_inactive % elapsed + gen_per_batch) / time_inactive_elapsed() speed_active = real(n_particles * n_active * gen_per_batch) / & - time_active % elapsed + time_active_elapsed() else speed_inactive = ZERO speed_active = real(n_particles * (n_batches - restart_batch) * & - gen_per_batch) / time_active % elapsed + gen_per_batch) / time_active_elapsed() end if else if (n_inactive > 0) then speed_inactive = real(n_particles * n_inactive * gen_per_batch) / & - time_inactive % elapsed + time_inactive_elapsed() end if speed_active = real(n_particles * n_active * gen_per_batch) / & - time_active % elapsed + time_active_elapsed() end if ! display calculation rate @@ -566,7 +566,7 @@ contains ! leakage rate. !=============================================================================== - subroutine print_results() + subroutine print_results() bind(C) integer :: n ! number of realizations real(8) :: alpha ! significance level for CI @@ -633,7 +633,7 @@ contains ! tallies and their standard deviations !=============================================================================== - subroutine write_tallies() + subroutine write_tallies() bind(C) integer :: i ! index in tallies array integer :: j ! level in tally hierarchy diff --git a/src/output.cpp b/src/output.cpp index fba6e5af7b..6a3abeed65 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -41,8 +41,7 @@ header(const char* msg, int level) { //============================================================================== -extern "C" void -print_overlap_check() { +void print_overlap_check() { #ifdef OPENMC_MPI std::vector temp(overlap_check_count); int err = MPI_Reduce(temp.data(), overlap_check_count.data(), diff --git a/src/particle.cpp b/src/particle.cpp index d6dfa969c1..c8a5a1b567 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -7,6 +7,7 @@ #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/hdf5_interface.h" +#include "openmc/mgxs_interface.h" #include "openmc/settings.h" #include "openmc/simulation.h" @@ -92,7 +93,7 @@ Particle::initialize() } void -Particle::from_source(const Bank* src, bool run_CE, const double* energy_bin_avg) +Particle::from_source(const Bank* src) { // set defaults initialize(); @@ -106,7 +107,7 @@ Particle::from_source(const Bank* src, bool run_CE, const double* energy_bin_avg std::copy(src->xyz, src->xyz + 3, last_xyz_current); std::copy(src->xyz, src->xyz + 3, last_xyz); std::copy(src->uvw, src->uvw + 3, last_uvw); - if (run_CE) { + if (settings::run_CE) { E = src->E; g = 0; } else { @@ -212,10 +213,9 @@ void particle_create_secondary(Particle* p, const double* uvw, double E, p->create_secondary(uvw, E, type, run_CE); } void particle_initialize(Particle* p) { p->initialize(); } -void particle_from_source(Particle* p, const Bank* src, bool run_CE, - const double* energy_bin_avg) +void particle_from_source(Particle* p, const Bank* src) { - p->from_source(src, run_CE, energy_bin_avg); + p->from_source(src); } void particle_mark_as_lost(Particle* p, const char* message) { diff --git a/src/particle_header.F90 b/src/particle_header.F90 index df5deff8c8..db2317f4b2 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -129,12 +129,10 @@ module particle_header type(Particle), intent(inout) :: p end subroutine particle_initialize - subroutine particle_from_source(p, src, run_CE, energy_bin_avg) bind(C) - import Particle, Bank, C_BOOL, C_DOUBLE + subroutine particle_from_source(p, src) bind(C) + import Particle, Bank, C_DOUBLE type(Particle), intent(inout) :: p type(Bank), intent(in) :: src - logical(C_BOOL), value :: run_CE - real(C_DOUBLE), intent(in) :: energy_bin_avg(*) end subroutine particle_from_source subroutine particle_mark_as_lost_c(p, message) bind(C, name='particle_mark_as_lost') diff --git a/src/settings.F90 b/src/settings.F90 index 5929bdd939..19ec3f7714 100644 --- a/src/settings.F90 +++ b/src/settings.F90 @@ -3,7 +3,6 @@ module settings use, intrinsic :: ISO_C_BINDING use constants - use set_header, only: SetInt implicit none @@ -67,9 +66,6 @@ module settings integer(C_INT32_T), bind(C) :: index_ufs_mesh ! Write source at end of simulation - logical(C_BOOL), bind(C) :: source_separate - logical(C_BOOL), bind(C) :: source_write - logical(C_BOOL), bind(C) :: source_latest ! Variance reduction settins logical(C_BOOL), bind(C) :: survival_biasing @@ -99,25 +95,10 @@ module settings ! Particle tracks logical(C_BOOL), bind(C) :: write_all_tracks - integer, allocatable :: track_identifiers(:,:) - - ! Particle restart run - logical(C_BOOL), bind(C) :: particle_restart_run - - ! Write out initial source - logical(C_BOOL), bind(C) :: write_initial_source ! Whether create fission neutrons or not. Only applied for MODE_FIXEDSOURCE logical(C_BOOL), bind(C) :: create_fission_neutrons - ! Information about state points to be written - integer :: n_state_points = 0 - type(SetInt) :: statepoint_batch - - ! Information about source points to be written - integer :: n_source_points = 0 - type(SetInt) :: sourcepoint_batch - character(MAX_FILE_LEN) :: path_input ! Path to input file character(MAX_FILE_LEN) :: path_cross_sections = '' ! Path to cross_sections.xml character(MAX_FILE_LEN) :: path_multipole ! Path to wmp library @@ -128,7 +109,6 @@ module settings ! Various output options logical(C_BOOL), bind(C) :: output_summary - logical(C_BOOL), bind(C) :: output_tallies ! Resonance scattering settings logical(C_BOOL), bind(C) :: res_scat_on ! is resonance scattering treated? @@ -150,11 +130,14 @@ contains !=============================================================================== subroutine free_memory_settings() - if (allocated(res_scat_nuclides)) deallocate(res_scat_nuclides) - if (allocated(track_identifiers)) deallocate(track_identifiers) + interface + subroutine free_memory_settings_c() bind(C) + end subroutine + end interface - call statepoint_batch % clear() - call sourcepoint_batch % clear() + if (allocated(res_scat_nuclides)) deallocate(res_scat_nuclides) + + call free_memory_settings_c() end subroutine free_memory_settings end module settings diff --git a/src/settings.cpp b/src/settings.cpp index bf97e44b88..ce09259785 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -9,9 +9,11 @@ #include "openmc/capi.h" #include "openmc/constants.h" +#include "openmc/container_util.h" #include "openmc/distribution.h" #include "openmc/distribution_multi.h" #include "openmc/distribution_spatial.h" +#include "openmc/eigenvalue.h" #include "openmc/error.h" #include "openmc/file_utils.h" #include "openmc/mesh.h" @@ -77,7 +79,7 @@ int32_t gen_per_batch {1}; int64_t n_particles {-1}; int electron_treatment {ELECTRON_TTB}; -double energy_cutoff[4] {0.0, 1000.0, 0.0, 0.0}; +std::array energy_cutoff {0.0, 1000.0, 0.0, 0.0}; int legendre_to_tabular_points {C_NONE}; int max_order {0}; int n_log_bins {8000}; @@ -86,13 +88,16 @@ int res_scat_method {RES_SCAT_ARES}; double res_scat_energy_min {0.01}; double res_scat_energy_max {1000.0}; int run_mode {-1}; +std::unordered_set sourcepoint_batch; +std::unordered_set statepoint_batch; int temperature_method {TEMPERATURE_NEAREST}; double temperature_tolerance {10.0}; double temperature_default {293.6}; -double temperature_range[2] {0.0, 0.0}; +std::array temperature_range {0.0, 0.0}; int trace_batch; int trace_gen; int64_t trace_particle; +std::vector> track_identifiers; int trigger_batch_interval {1}; int verbosity {7}; double weight_cutoff {0.25}; @@ -141,7 +146,10 @@ void get_run_parameters(pugi::xml_node node_base) gen_per_batch = std::stoi(get_node_value(node_base, "generations_per_batch")); } - // TODO: Preallocate space for keff and entropy by generation + // Preallocate space for keff and entropy by generation + int m = settings::n_max_batches * settings::gen_per_batch; + simulation::k_generation.reserve(m); + entropy.reserve(m); // Get the trigger information for keff if (check_for_node(node_base, "keff_trigger")) { @@ -491,7 +499,7 @@ void read_settings_xml() // Particle tracks if (check_for_node(root, "track")) { // Get values and make sure there are three per particle - auto temp = get_node_array(root, "track"); + auto temp = get_node_array(root, "track"); if (temp.size() % 3 != 0) { fatal_error("Number of integers specified in 'track' is not " "divisible by 3. Please provide 3 integers per particle to be " @@ -499,8 +507,11 @@ void read_settings_xml() } // Reshape into track_identifiers - //allocate(track_identifiers(3, n_tracks/3)) - //track_identifiers = reshape(temp_int_array, [3, n_tracks/3]) + int n_tracks = temp.size() / 3; + for (int i = 0; i < n_tracks; ++i) { + track_identifiers.push_back({temp[3*i], temp[3*i + 1], + temp[3*i + 2]}); + } } // Read meshes @@ -539,7 +550,7 @@ void read_settings_xml() // If the user did not specify how many mesh cells are to be used in // each direction, we automatically determine an appropriate number of // cells - int n = std::ceil(std::pow(settings::n_particles / 20.0, 1.0/3.0)); + int n = std::ceil(std::pow(n_particles / 20.0, 1.0/3.0)); m.shape_ = {n, n, n}; m.n_dimension_ = 3; @@ -548,7 +559,7 @@ void read_settings_xml() } // Turn on Shannon entropy calculation - settings::entropy_on = true; + entropy_on = true; } // Uniform fission source weighting mesh @@ -581,18 +592,48 @@ void read_settings_xml() if (index_ufs_mesh >= 0) { // Turn on uniform fission source weighting - settings::ufs_on = true; + ufs_on = true; } - // TODO: Read + // Check if the user has specified to write state points + if (check_for_node(root, "state_point")) { + // Get pointer to state_point node + auto node_sp = root.child("state_point"); + + // Determine number of batches at which to store state points + if (check_for_node(node_sp, "batches")) { + // User gave specific batches to write state points + auto temp = get_node_array(node_sp, "batches"); + for (const auto& b : temp) { + statepoint_batch.insert(b); + } + } else { + // If neither were specified, write state point at last batch + statepoint_batch.insert(n_batches); + } + } else { + // If no tag was present, by default write state point at + // last batch only + statepoint_batch.insert(n_batches); + } // Check if the user has specified to write source points if (check_for_node(root, "source_point")) { // Get source_point node xml_node node_sp = root.child("source_point"); - // TODO: Read source point batches + // Determine batches at which to store source points + if (check_for_node(node_sp, "batches")) { + // User gave specific batches to write source points + auto temp = get_node_array(node_sp, "batches"); + for (const auto& b : temp) { + sourcepoint_batch.insert(b); + } + } else { + // If neither were specified, write source points with state points + sourcepoint_batch = statepoint_batch; + } // Check if the user has specified to write binary source file if (check_for_node(node_sp, "separate")) { @@ -609,10 +650,19 @@ void read_settings_xml() // If no tag was present, by default we keep source bank in // statepoint file and write it out at statepoints intervals source_separate = false; - // TODO: add defaults + sourcepoint_batch = statepoint_batch; } - // TODO: Check source points are subset + // If source is not seperate and is to be written out in the statepoint file, + // make sure that the sourcepoint batch numbers are contained in the + // statepoint list + if (!source_separate) { + for (const auto& b : sourcepoint_batch) { + if (!contains(statepoint_batch, b)) { + fatal_error("Sourcepoint batches are not a subset of statepoint batches."); + } + } + } // Check if the user has specified to not reduce tallies at the end of every // batch @@ -778,6 +828,11 @@ extern "C" { const char* openmc_path_particle_restart() { return settings::path_particle_restart.c_str(); } + + void free_memory_settings_c() { + settings::statepoint_batch.clear(); + settings::sourcepoint_batch.clear(); + } } } // namespace openmc diff --git a/src/simulation.F90 b/src/simulation.F90 index e744f1ea82..3eabbfecb4 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -7,392 +7,30 @@ module simulation #endif use bank_header, only: source_bank - use cmfd_execute, only: cmfd_init_batch, cmfd_tally_init, execute_cmfd - use cmfd_header, only: cmfd_on use constants, only: ZERO - use eigenvalue, only: calculate_average_keff, calculate_generation_keff, & - k_sum -#ifdef _OPENMP - use eigenvalue, only: join_bank_from_threads -#endif - use error, only: fatal_error, write_message - use geometry_header, only: n_cells + use error, only: fatal_error use material_header, only: n_materials, materials use message_passing - use mgxs_interface, only: energy_bins, energy_bin_avg use nuclide_header, only: micro_xs, n_nuclides - use output, only: header, print_columns, & - print_batch_keff, print_generation, print_runtime, & - print_results, write_tallies - use particle_header use photon_header, only: micro_photon_xs, n_elements - use random_lcg, only: set_particle_seed use settings use simulation_header - use state_point, only: openmc_statepoint_write, write_source_point, load_state_point - use string, only: to_str - use tally, only: accumulate_tallies, setup_active_tallies, & - init_tally_routines use tally_header use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer - use tally_derivative_header, only: tally_derivs - use timer_header - use trigger, only: check_triggers - use tracking, only: transport implicit none private - public :: openmc_next_batch - public :: openmc_simulation_init - public :: openmc_simulation_finalize - - integer(C_INT), parameter :: STATUS_EXIT_NORMAL = 0 - integer(C_INT), parameter :: STATUS_EXIT_MAX_BATCH = 1 - integer(C_INT), parameter :: STATUS_EXIT_ON_TRIGGER = 2 - - interface - subroutine openmc_simulation_init_c() bind(C) - end subroutine - - subroutine initialize_source() bind(C) - end subroutine - - subroutine initialize_generation() bind(C) - end subroutine - - function sample_external_source() result(site) bind(C) - import Bank - type(Bank) :: site - end function - end interface contains -!=============================================================================== -! OPENMC_NEXT_BATCH -!=============================================================================== - - function openmc_next_batch(status) result(err) bind(C) - integer(C_INT), intent(out), optional :: status - integer(C_INT) :: err - - type(Particle) :: p - integer(8) :: i_work - - err = 0 - - ! Make sure simulation has been initialized - if (.not. simulation_initialized) then - err = E_ALLOCATE - call set_errmsg("Simulation has not been initialized yet.") - return - end if - - call initialize_batch() - - ! ======================================================================= - ! LOOP OVER GENERATIONS - GENERATION_LOOP: do current_gen = 1, gen_per_batch - - call initialize_generation() - - ! Start timer for transport - call time_transport % start() - - ! ==================================================================== - ! LOOP OVER PARTICLES -!$omp parallel do schedule(runtime) firstprivate(p) copyin(tally_derivs) - PARTICLE_LOOP: do i_work = 1, work - current_work = i_work - - ! grab source particle from bank - call initialize_history(p, current_work) - - ! transport particle - call transport(p) - - end do PARTICLE_LOOP -!$omp end parallel do - - ! Accumulate time for transport - call time_transport % stop() - - call finalize_generation() - - end do GENERATION_LOOP - - call finalize_batch() - - ! Check simulation ending criteria - if (present(status)) then - if (current_batch == n_max_batches) then - status = STATUS_EXIT_MAX_BATCH - elseif (satisfy_triggers) then - status = STATUS_EXIT_ON_TRIGGER - else - status = STATUS_EXIT_NORMAL - end if - end if - - end function openmc_next_batch - -!=============================================================================== -! INITIALIZE_HISTORY -!=============================================================================== - - subroutine initialize_history(p, index_source) - - type(Particle), intent(inout) :: p - integer(8), intent(in) :: index_source - - integer(8) :: particle_seed ! unique index for particle - integer :: i - - ! set defaults - call particle_from_source(p, source_bank(index_source), run_CE, & - energy_bin_avg) - - ! set identifier for particle - p % id = work_index(rank) + index_source - - ! set random number seed - particle_seed = (total_gen + overall_generation() - 1)*n_particles + p % id - call set_particle_seed(particle_seed) - - ! set particle trace - trace = .false. - if (current_batch == trace_batch .and. current_gen == trace_gen .and. & - p % id == trace_particle) trace = .true. - - ! Set particle track. - p % write_track = .false. - if (write_all_tracks) then - p % write_track = .true. - else if (allocated(track_identifiers)) then - do i=1, size(track_identifiers(1,:)) - if (current_batch == track_identifiers(1,i) .and. & - ¤t_gen == track_identifiers(2,i) .and. & - &p % id == track_identifiers(3,i)) then - p % write_track = .true. - exit - end if - end do - end if - - end subroutine initialize_history - -!=============================================================================== -! INITIALIZE_BATCH -!=============================================================================== - - subroutine initialize_batch() - - integer :: i - - ! Increment current batch - current_batch = current_batch + 1 - - if (run_mode == MODE_FIXEDSOURCE) then - call write_message("Simulating batch " // trim(to_str(current_batch)) & - // "...", 6) - end if - - ! Reset total starting particle weight used for normalizing tallies - total_weight = ZERO - - if ((n_inactive > 0 .and. current_batch == 1) .or. & - (restart_run .and. restart_batch < n_inactive .and. current_batch == restart_batch + 1)) then - ! Turn on inactive timer - call time_inactive % start() - elseif ((current_batch == n_inactive + 1) .or. & - (restart_run .and. restart_batch > n_inactive .and. current_batch == restart_batch + 1)) then - ! Switch from inactive batch timer to active batch timer - call time_inactive % stop() - call time_active % start() - - do i = 1, n_tallies - tallies(i) % obj % active = .true. - end do - end if - - ! check CMFD initialize batch - if (run_mode == MODE_EIGENVALUE) then - if (cmfd_run) call cmfd_init_batch() - end if - - ! Add user tallies to active tallies list - call setup_active_tallies() - - end subroutine initialize_batch - -!=============================================================================== -! FINALIZE_GENERATION -!=============================================================================== - - subroutine finalize_generation() - - interface - subroutine fill_source_bank_fixedsource() bind(C) - end subroutine - - subroutine shannon_entropy() bind(C) - end subroutine - - subroutine synchronize_bank() bind(C) - end subroutine - end interface - - ! Update global tallies with the omp private accumulation variables -!$omp parallel -!$omp critical - if (run_mode == MODE_EIGENVALUE) then - global_tallies(RESULT_VALUE, K_COLLISION) = & - global_tallies(RESULT_VALUE, K_COLLISION) + global_tally_collision - global_tallies(RESULT_VALUE, K_ABSORPTION) = & - global_tallies(RESULT_VALUE, K_ABSORPTION) + global_tally_absorption - global_tallies(RESULT_VALUE, K_TRACKLENGTH) = & - global_tallies(RESULT_VALUE, K_TRACKLENGTH) + global_tally_tracklength - end if - global_tallies(RESULT_VALUE, LEAKAGE) = & - global_tallies(RESULT_VALUE, LEAKAGE) + global_tally_leakage -!$omp end critical - - ! reset private tallies - if (run_mode == MODE_EIGENVALUE) then - global_tally_collision = ZERO - global_tally_absorption = ZERO - global_tally_tracklength = ZERO - end if - global_tally_leakage = ZERO -!$omp end parallel - - if (run_mode == MODE_EIGENVALUE) then -#ifdef _OPENMP - ! Join the fission bank from each thread into one global fission bank - call join_bank_from_threads() -#endif - - ! Distribute fission bank across processors evenly - call synchronize_bank() - - ! Calculate shannon entropy - if (entropy_on) call shannon_entropy() - - ! Collect results and statistics - call calculate_generation_keff() - call calculate_average_keff() - - ! Write generation output - if (master .and. verbosity >= 7) then - if (current_gen /= gen_per_batch) then - call print_generation() - end if - end if - - elseif (run_mode == MODE_FIXEDSOURCE) then - ! For fixed-source mode, we need to sample the external source - call fill_source_bank_fixedsource() - end if - - end subroutine finalize_generation - -!=============================================================================== -! FINALIZE_BATCH handles synchronization and accumulation of tallies, -! calculation of Shannon entropy, getting single-batch estimate of keff, and -! turning on tallies when appropriate -!=============================================================================== - - subroutine finalize_batch() - - integer(C_INT) :: err - character(MAX_FILE_LEN) :: filename - - interface - subroutine broadcast_triggers() bind(C) - end subroutine broadcast_triggers - end interface - - ! Reduce tallies onto master process and accumulate - call time_tallies % start() - call accumulate_tallies() - call time_tallies % stop() - - ! Reset global tally results - if (current_batch <= n_inactive) then - global_tallies(:,:) = ZERO - n_realizations = 0 - end if - - if (run_mode == MODE_EIGENVALUE) then - ! Perform CMFD calculation if on - if (cmfd_on) call execute_cmfd() - ! Write batch output - if (master .and. verbosity >= 7) call print_batch_keff() - end if - - ! Check_triggers - if (master) call check_triggers() -#ifdef OPENMC_MPI - call broadcast_triggers() -#endif - if (satisfy_triggers .or. & - (trigger_on .and. current_batch == n_max_batches)) then - call statepoint_batch % add(current_batch) - end if - - ! Write out state point if it's been specified for this batch - if (statepoint_batch % contains(current_batch)) then - if (sourcepoint_batch % contains(current_batch) .and. source_write & - .and. .not. source_separate) then - err = openmc_statepoint_write(write_source=.true._C_BOOL) - else - err = openmc_statepoint_write(write_source=.false._C_BOOL) - end if - end if - - ! Write out a separate source point if it's been specified for this batch - if (sourcepoint_batch % contains(current_batch) .and. source_write & - .and. source_separate) call write_source_point() - - ! Write a continously-overwritten source point if requested. - if (source_latest) then - filename = trim(path_output) // 'source' // '.h5' - call write_source_point(filename) - end if - - end subroutine finalize_batch - - !=============================================================================== ! INITIALIZE_SIMULATION !=============================================================================== - function openmc_simulation_init() result(err) bind(C) - integer(C_INT) :: err + subroutine simulation_init_f() bind(C) integer :: i - err = 0 - - ! Skip if simulation has already been initialized - if (simulation_initialized) return - - ! Call initialization on C++ side - call openmc_simulation_init_c() - - ! Set up tally procedure pointers - call init_tally_routines() - - ! Allocate source bank, and for eigenvalue simulations also allocate the - ! fission bank - call allocate_banks() - - ! Allocate tally results arrays if they're not allocated yet - call configure_tallies() - - ! Activate the CMFD tallies - call cmfd_tally_init() - ! Set up material nuclide index mapping do i = 1, n_materials call materials(i) % init_nuclide_index() @@ -410,66 +48,16 @@ contains end do !$omp end parallel - ! Reset global variables -- this is done before loading state point (as that - ! will potentially populate k_generation and entropy) - current_batch = 0 - call k_generation_clear() - call entropy_clear() - need_depletion_rx = .false. - - ! If this is a restart run, load the state point data and binary source - ! file - if (restart_run) then - call load_state_point() - call write_message("Resuming simulation...", 6) - else - call initialize_source() - end if - - ! Display header - if (master) then - if (run_mode == MODE_FIXEDSOURCE) then - call header("FIXED SOURCE TRANSPORT SIMULATION", 3) - elseif (run_mode == MODE_EIGENVALUE) then - call header("K EIGENVALUE SIMULATION", 3) - if (verbosity >= 7) call print_columns() - end if - end if - - ! Set flag indicating initialization is done - simulation_initialized = .true. - - end function openmc_simulation_init + end subroutine !=============================================================================== ! FINALIZE_SIMULATION calculates tally statistics, writes tallies, and displays ! execution time and results !=============================================================================== - function openmc_simulation_finalize() result(err) bind(C) - integer(C_INT) :: err + subroutine simulation_finalize_f() bind(C) - integer :: i ! loop index - - interface - subroutine openmc_simulation_finalize_c() bind(C) - end subroutine openmc_simulation_finalize_c - - subroutine print_overlap_check() bind(C) - end subroutine print_overlap_check - - subroutine broadcast_results() bind(C) - end subroutine broadcast_results - end interface - - err = 0 - - ! Skip if simulation was never run - if (.not. simulation_initialized) return - - ! Stop active batch timer and start finalization timer - call time_active % stop() - call time_finalize % start() + integer :: i ! loop index ! Free up simulation-specific memory do i = 1, n_materials @@ -479,45 +67,13 @@ contains deallocate(micro_xs, micro_photon_xs, filter_matches) !$omp end parallel - ! Increment total number of generations - total_gen = total_gen + current_batch*gen_per_batch - -#ifdef OPENMC_MPI - call broadcast_results() -#endif - - ! Write tally results to tallies.out - if (output_tallies .and. master) call write_tallies() - - ! Deactivate all tallies - if (allocated(tallies)) then - do i = 1, n_tallies - tallies(i) % obj % active = .false. - end do - end if - - call openmc_simulation_finalize_c() - - ! Stop timers and show timing statistics - call time_finalize%stop() - call time_total%stop() - if (master) then - if (verbosity >= 6) call print_runtime() - if (verbosity >= 4) call print_results() - end if - if (check_overlaps) call print_overlap_check() - - ! Reset flags - need_depletion_rx = .false. - simulation_initialized = .false. - - end function openmc_simulation_finalize + end subroutine !=============================================================================== ! ALLOCATE_BANKS allocates memory for the fission and source banks !=============================================================================== - subroutine allocate_banks() + subroutine allocate_banks() bind(C) integer :: alloc_err ! allocation error code diff --git a/src/simulation.cpp b/src/simulation.cpp index 36f9811094..84ebb01921 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1,13 +1,54 @@ #include "openmc/simulation.h" #include "openmc/capi.h" +#include "openmc/container_util.h" #include "openmc/eigenvalue.h" +#include "openmc/error.h" #include "openmc/message_passing.h" +#include "openmc/output.h" +#include "openmc/particle.h" +#include "openmc/random_lcg.h" #include "openmc/settings.h" +#include "openmc/source.h" +#include "openmc/state_point.h" +#include "openmc/timer.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/tally.h" #include +#include + +namespace openmc { + +// data/functions from Fortran side +extern "C" bool cmfd_on; + +extern "C" void accumulate_tallies(); +extern "C" void allocate_banks(); +extern "C" void check_triggers(); +extern "C" void cmfd_init_batch(); +extern "C" void cmfd_tally_init(); +extern "C" void configure_tallies(); +extern "C" void execute_cmfd(); +extern "C" void init_tally_routines(); +extern "C" void join_bank_from_threads(); +extern "C" void load_state_point(); +extern "C" void print_batch_keff(); +extern "C" void print_columns(); +extern "C" void print_generation(); +extern "C" void print_results(); +extern "C" void print_runtime(); +extern "C" void setup_active_tallies(); +extern "C" void simulation_init_f(); +extern "C" void simulation_finalize_f(); +extern "C" void transport(Particle* p); +extern "C" void write_tallies(); + +} // namespace openmc + +//============================================================================== +// C API functions +//============================================================================== // OPENMC_RUN encompasses all the main logic where iterations are performed // over the batches, generations, and histories in a fixed source or k-eigenvalue @@ -27,6 +68,176 @@ int openmc_run() return err; } +int openmc_simulation_init() +{ + using namespace openmc; + + // Skip if simulation has already been initialized + if (simulation::initialized) return 0; + + // Determine how much work each process should do + calculate_work(); + + // Allocate array for matching filter bins +#pragma omp parallel + { + filter_matches.resize(n_filters); + } + + // Set up tally procedure pointers + init_tally_routines(); + + // Allocate source bank, and for eigenvalue simulations also allocate the + // fission bank + allocate_banks(); + + // Allocate tally results arrays if they're not allocated yet + configure_tallies(); + + // Activate the CMFD tallies + cmfd_tally_init(); + + // Call Fortran initialization + simulation_init_f(); + + // Reset global variables -- this is done before loading state point (as that + // will potentially populate k_generation and entropy) + simulation::current_batch = 0; + simulation::k_generation.clear(); + entropy.clear(); + simulation::need_depletion_rx = false; + + // If this is a restart run, load the state point data and binary source + // file + if (settings::restart_run) { + load_state_point(); + write_message("Resuming simulation...", 6); + } else { + initialize_source(); + } + + // Display header + if (mpi::master) { + if (settings::run_mode == RUN_MODE_FIXEDSOURCE) { + header("FIXED SOURCE TRANSPORT SIMULATION", 3); + } else if (settings::run_mode == RUN_MODE_EIGENVALUE) { + header("K EIGENVALUE SIMULATION", 3); + if (settings::verbosity >= 7) print_columns(); + } + } + + // Set flag indicating initialization is done + simulation::initialized = true; + return 0; +} + +int openmc_simulation_finalize() +{ + using namespace openmc; + + // Skip if simulation was never run + if (!simulation::initialized) return 0; + + // Stop active batch timer and start finalization timer + time_active.stop(); + time_finalize.start(); + +#pragma omp parallel + { + filter_matches.clear(); + } + + // Deallocate Fortran variables, set tallies to inactive + simulation_finalize_f(); + + // Increment total number of generations + simulation::total_gen += simulation::current_batch*settings::gen_per_batch; + +#ifdef OPENMC_MPI + broadcast_results(); +#endif + + // Write tally results to tallies.out + if (settings::output_tallies && mpi::master) write_tallies(); + + // Deactivate all tallies + for (int i = 1; i <= n_tallies; ++i) { + openmc_tally_set_active(i, false); + } + + // Stop timers and show timing statistics + time_finalize.stop(); + time_total.stop(); + if (mpi::master) { + if (settings::verbosity >= 6) print_runtime(); + if (settings::verbosity >= 4) print_results(); + } + if (settings::check_overlaps) print_overlap_check(); + + // Reset flags + simulation::need_depletion_rx = false; + simulation::initialized = false; + return 0; +} + +int openmc_next_batch(int* status) +{ + using namespace openmc; + using openmc::simulation::current_gen; + + // Make sure simulation has been initialized + if (!simulation::initialized) { + set_errmsg("Simulation has not been initialized yet."); + return OPENMC_E_ALLOCATE; + } + + initialize_batch(); + + // ======================================================================= + // LOOP OVER GENERATIONS + for (current_gen = 1; current_gen <= settings::gen_per_batch; ++current_gen) { + + initialize_generation(); + + // Start timer for transport + time_transport.start(); + + // ==================================================================== + // LOOP OVER PARTICLES + +#pragma omp parallel for schedule(runtime) + for (int64_t i_work = 1; i_work <= simulation::work; ++i_work) { + simulation::current_work = i_work; + + // grab source particle from bank + Particle p; + initialize_history(&p, simulation::current_work); + + // transport particle + transport(&p); + } + + // Accumulate time for transport + time_transport.stop(); + + finalize_generation(); + } + + finalize_batch(); + + // Check simulation ending criteria + if (status) { + if (simulation::current_batch == settings::n_max_batches) { + *status = STATUS_EXIT_MAX_BATCH; + } else if (simulation::satisfy_triggers) { + *status = STATUS_EXIT_ON_TRIGGER; + } else { + *status = STATUS_EXIT_NORMAL; + } + } + return 0; +} + namespace openmc { //============================================================================== @@ -38,6 +249,7 @@ namespace simulation { int current_batch; int current_gen; int64_t current_work; +bool initialized {false}; double keff {1.0}; double keff_std; double k_col_abs {0.0}; @@ -48,7 +260,6 @@ int n_lost_particles {0}; bool need_depletion_rx {false}; int restart_batch; bool satisfy_triggers {false}; -bool simulation_initialized {false}; int total_gen {0}; int64_t work; @@ -68,15 +279,104 @@ int thread_id; //!< ID of a given thread // Non-member functions //============================================================================== -void openmc_simulation_init_c() +void initialize_batch() { - // Determine how much work each process should do - calculate_work(); + // Increment current batch + ++simulation::current_batch; - // Allocate array for matching filter bins - #pragma omp parallel - { - filter_matches.resize(n_filters); + if (settings::run_mode == RUN_MODE_FIXEDSOURCE) { + int b = simulation::current_batch; + write_message("Simulating batch " + std::to_string(b), 6); + } + + // Reset total starting particle weight used for normalizing tallies + total_weight = 0.0; + + // Determine if this batch is the first inactive or active batch. + bool first_inactive = false; + bool first_active = false; + if (!settings::restart_run) { + first_inactive = settings::n_inactive > 0 && simulation::current_batch == 1; + first_active = simulation::current_batch == settings::n_inactive + 1; + } else if (simulation::current_batch == simulation::restart_batch + 1){ + first_inactive = simulation::restart_batch < settings::n_inactive; + first_active = !first_inactive; + } + + // Manage active/inactive timers and activate tallies if necessary. + if (first_inactive) { + time_inactive.start(); + } else if (first_active) { + time_inactive.stop(); + time_active.start(); + for (int i = 1; i <= n_tallies; ++i) { + // TODO: change one-based index + openmc_tally_set_active(i, true); + } + } + + // check CMFD initialize batch + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + if (settings::cmfd_run) cmfd_init_batch(); + } + + // Add user tallies to active tallies list + setup_active_tallies(); +} + +void finalize_batch() +{ + // Reduce tallies onto master process and accumulate + time_tallies.start(); + accumulate_tallies(); + time_tallies.stop(); + + // Reset global tally results + if (simulation::current_batch <= settings::n_inactive) { + auto gt = global_tallies(); + std::fill(gt.begin(), gt.end(), 0.0); + n_realizations = 0; + } + + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + // Perform CMFD calculation if on + if (cmfd_on) execute_cmfd(); + // Write batch output + if (mpi::master && settings::verbosity >= 7) print_batch_keff(); + } + + // Check_triggers + if (mpi::master) check_triggers(); +#ifdef OPENMC_MPI + MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm); +#endif + if (simulation::satisfy_triggers || (settings::trigger_on && + simulation::current_batch == settings::n_max_batches)) { + settings::statepoint_batch.insert(simulation::current_batch); + } + + // Write out state point if it's been specified for this batch + if (contains(settings::statepoint_batch, simulation::current_batch)) { + if (contains(settings::sourcepoint_batch, simulation::current_batch) + && settings::source_write && !settings::source_separate) { + bool b = true; + openmc_statepoint_write(nullptr, &b); + } else { + bool b = false; + openmc_statepoint_write(nullptr, &b); + } + } + + // 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()); } } @@ -94,12 +394,99 @@ void initialize_generation() } } -extern "C" void -openmc_simulation_finalize_c() +void finalize_generation() { - #pragma omp parallel + auto gt = global_tallies(); + + // Update global tallies with the omp private accumulation variables +#pragma omp parallel { - filter_matches.clear(); +#pragma omp critical(increment_global_tallies) + { + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + gt(K_COLLISION, RESULT_VALUE) += global_tally_collision; + gt(K_ABSORPTION, RESULT_VALUE) += global_tally_absorption; + gt(K_TRACKLENGTH, RESULT_VALUE) += global_tally_tracklength; + } + gt(LEAKAGE, RESULT_VALUE) += global_tally_leakage; + } + + // reset threadprivate tallies + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + global_tally_collision = 0.0; + global_tally_absorption = 0.0; + global_tally_tracklength = 0.0; + } + global_tally_leakage = 0.0; + } + + + if (settings::run_mode == RUN_MODE_EIGENVALUE) { +#ifdef _OPENMP + // Join the fission bank from each thread into one global fission bank + join_bank_from_threads(); +#endif + + // Distribute fission bank across processors evenly + synchronize_bank(); + + // Calculate shannon entropy + if (settings::entropy_on) shannon_entropy(); + + // Collect results and statistics + calculate_generation_keff(); + calculate_average_keff(); + + // Write generation output + if (mpi::master && settings::verbosity >= 7) { + if (simulation::current_gen != settings::gen_per_batch) { + print_generation(); + } + } + + } else if (settings::run_mode == RUN_MODE_FIXEDSOURCE) { + // For fixed-source mode, we need to sample the external source + fill_source_bank_fixedsource(); + } +} + +void initialize_history(Particle* p, int64_t index_source) +{ + // Get pointer to source bank + Bank* source_bank; + int64_t n; + openmc_source_bank(&source_bank, &n); + + // set defaults + p->from_source(&source_bank[index_source - 1]); + + // set identifier for particle + p->id = simulation::work_index[mpi::rank] + index_source; + + // set random number seed + int64_t particle_seed = (simulation::total_gen + overall_generation() - 1) + * settings::n_particles + p->id; + set_particle_seed(particle_seed); + + // set particle trace + simulation::trace = false; + if (simulation::current_batch == settings::trace_batch && + simulation::current_gen == settings::trace_gen && + p->id == settings::trace_particle) simulation::trace = true; + + // Set particle track. + p->write_track = false; + if (settings::write_all_tracks) { + p->write_track = true; + } else if (settings::track_identifiers.size() > 0) { + for (const auto& t : settings::track_identifiers) { + if (simulation::current_batch == t[0] && + simulation::current_gen == t[1] && + p->id == t[2]) { + p->write_track = true; + break; + } + } } } @@ -165,10 +552,6 @@ void broadcast_results() { simulation::k_abs_tra = temp[2]; } -void broadcast_triggers() -{ - MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm); -} #endif //============================================================================== diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 1844c8da6e..1df4b0a4e3 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -23,7 +23,6 @@ module simulation_header integer(C_INT), bind(C) :: current_batch ! current batch integer(C_INT), bind(C) :: current_gen ! current generation within a batch integer(C_INT), bind(C) :: total_gen ! total number of generations simulated - logical(C_BOOL), bind(C) :: simulation_initialized logical(C_BOOL), bind(C) :: need_depletion_rx ! need to calculate depletion reaction rx? ! ============================================================================ @@ -32,7 +31,6 @@ module simulation_header logical(C_BOOL), bind(C) :: satisfy_triggers ! whether triggers are satisfied integer(C_INT64_T), bind(C) :: work ! number of particles per processor - integer(C_INT64_T), bind(C) :: current_work ! index in source bank of current history simulated ! ============================================================================ ! K-EIGENVALUE SIMULATION VARIABLES @@ -59,7 +57,7 @@ module simulation_header logical(C_BOOL), bind(C) :: trace -!$omp threadprivate(trace, thread_id, current_work) +!$omp threadprivate(trace, thread_id) interface subroutine entropy_clear() bind(C) @@ -84,11 +82,6 @@ module simulation_header subroutine k_generation_clear() bind(C) end subroutine - subroutine k_generation_reserve(i) bind(C) - import C_INT - integer(C_INT), value :: i - end subroutine - function work_index(rank) result(i) bind(C) import C_INT, C_INT64_T integer(C_INT), value :: rank diff --git a/src/source.cpp b/src/source.cpp index 3427579ee4..bfed41b6a3 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -352,11 +352,7 @@ extern "C" double total_source_strength() return strength; } -// Needed in fill_source_bank_fixedsource -extern "C" int overall_generation(); - -//! Fill source bank at end of generation for fixed source simulations -extern "C" void fill_source_bank_fixedsource() +void fill_source_bank_fixedsource() { if (settings::path_source.empty()) { // Get pointer to source bank diff --git a/src/state_point.F90 b/src/state_point.F90 index 52b01557d4..1286bb63d4 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -16,7 +16,6 @@ module state_point use bank_header, only: Bank use cmfd_header use constants - use eigenvalue, only: openmc_get_keff, k_sum use endf, only: reaction_name use error, only: fatal_error, warning, write_message use hdf5_interface @@ -391,19 +390,19 @@ contains ! Write out the runtime metrics. runtime_group = create_group(file_id, "runtime") call write_dataset(runtime_group, "total initialization", & - time_initialize % get_value()) + time_initialize_elapsed()) call write_dataset(runtime_group, "reading cross sections", & time_read_xs % get_value()) call write_dataset(runtime_group, "simulation", & - time_inactive % get_value() + time_active % get_value()) + time_inactive_elapsed() + time_active_elapsed()) call write_dataset(runtime_group, "transport", & - time_transport % get_value()) + time_transport_elapsed()) if (run_mode == MODE_EIGENVALUE) then call write_dataset(runtime_group, "inactive batches", & - time_inactive % get_value()) + time_inactive_elapsed()) end if call write_dataset(runtime_group, "active batches", & - time_active % get_value()) + time_active_elapsed()) if (run_mode == MODE_EIGENVALUE) then call write_dataset(runtime_group, "synchronizing fission bank", & time_bank_elapsed()) @@ -413,7 +412,7 @@ contains time_bank_sendrecv_elapsed()) end if call write_dataset(runtime_group, "accumulating tallies", & - time_tallies % get_value()) + time_tallies_elapsed()) if (cmfd_run) then call write_dataset(runtime_group, "CMFD", time_cmfd % get_value()) call write_dataset(runtime_group, "CMFD building matrices", & @@ -421,7 +420,7 @@ contains call write_dataset(runtime_group, "CMFD solving matrices", & time_cmfdsolve % get_value()) end if - call write_dataset(runtime_group, "total", time_total % get_value()) + call write_dataset(runtime_group, "total", time_total_elapsed()) call close_group(runtime_group) call file_close(file_id) @@ -443,52 +442,13 @@ contains end if end function openmc_statepoint_write -!=============================================================================== -! WRITE_SOURCE_POINT -!=============================================================================== - - subroutine write_source_point(filename) - character(MAX_FILE_LEN), intent(in), optional :: filename - - logical :: parallel - integer(HID_T) :: file_id - character(MAX_FILE_LEN) :: filename_ - - ! When using parallel HDF5, the file is written to collectively by all - ! processes. With MPI-only, the file is opened and written by the master - ! (note that the call to write_source_bank is by all processes since slave - ! processes need to send source bank data to the master. -#ifdef PHDF5 - parallel = .true. -#else - parallel = .false. -#endif - - if (present(filename)) then - filename_ = filename - else - filename_ = trim(path_output) // 'source.' // & - & zero_padded(current_batch, count_digits(n_max_batches)) - filename_ = trim(filename_) // '.h5' - end if - - if (master .or. parallel) then - file_id = file_open(filename_, 'w', parallel=.true.) - call write_attribute(file_id, "filetype", 'source') - end if - call write_source_bank(file_id, source_bank) - if (master .or. parallel) call file_close(file_id) - - end subroutine write_source_point - !=============================================================================== ! LOAD_STATE_POINT !=============================================================================== - subroutine load_state_point() + subroutine load_state_point() bind(C) integer :: i - integer :: n integer :: int_array(3) integer, allocatable :: array(:) integer(C_INT64_T) :: seed @@ -504,6 +464,9 @@ contains import HID_T integer(HID_T), value :: group end subroutine + + subroutine restart_set_keff() bind(C) + end subroutine end interface ! Write message @@ -612,17 +575,7 @@ contains ! Set k_sum, keff, and current_batch based on whether restart file is part ! of active cycle or inactive cycle - if (restart_batch > n_inactive) then - do i = n_inactive + 1, restart_batch - k_sum(1) = k_sum(1) + k_generation(i) - k_sum(2) = k_sum(2) + k_generation(i)**2 - end do - n = gen_per_batch*n_realizations - keff = k_sum(1) / n - else - n = k_generation_size() - keff = k_generation(n) - end if + call restart_set_keff() current_batch = restart_batch ! Check to make sure source bank is present diff --git a/src/state_point.cpp b/src/state_point.cpp index 63c8264a36..620b5f71a7 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -1,6 +1,7 @@ #include "openmc/state_point.h" #include +#include // for setfill, setw #include #include @@ -9,6 +10,7 @@ #include "openmc/capi.h" #include "openmc/constants.h" +#include "openmc/eigenvalue.h" #include "openmc/error.h" #include "openmc/hdf5_interface.h" #include "openmc/message_passing.h" @@ -18,7 +20,6 @@ namespace openmc { - hid_t h5banktype() { // Create type for array of 3 reals hsize_t dims[] {3}; @@ -36,6 +37,46 @@ hid_t h5banktype() { return banktype; } +void +write_source_point(const char* filename) +{ + // When using parallel HDF5, the file is written to collectively by all + // processes. With MPI-only, the file is opened and written by the master + // (note that the call to write_source_bank is by all processes since slave + // processes need to send source bank data to the master. +#ifdef PHDF5 + bool parallel = true; +#else + bool parallel = false; +#endif + + std::string filename_; + if (filename) { + filename_ = filename; + } else { + // Determine width for zero padding + int w = std::to_string(settings::n_max_batches).size(); + + std::stringstream s; + s << settings::path_output << "source." << std::setfill('0') + << std::setw(w) << simulation::current_batch << ".h5"; + filename_ = s.str(); + } + + hid_t file_id; + if (mpi::master || parallel) { + file_id = file_open(filename_, 'w', true); + write_attribute(file_id, "filetype", "source"); + } + + // Get pointer to source bank and write to file + Bank* source_bank; + int64_t n; + openmc_source_bank(&source_bank, &n); + write_source_bank(file_id, source_bank); + + if (mpi::master || parallel) file_close(file_id); +} void write_source_bank(hid_t group_id, Bank* source_bank) @@ -275,4 +316,18 @@ 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) { + k_sum[0] += simulation::k_generation[i]; + k_sum[1] += std::pow(simulation::k_generation[i], 2); + } + int n = settings::gen_per_batch*n_realizations; + simulation::keff = k_sum[0] / n; + } else { + simulation::keff = simulation::k_generation.back(); + } +} + } // namespace openmc diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index c15325418d..ba152353e7 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -52,7 +52,7 @@ contains ! with the CE and MG modes. !=============================================================================== - subroutine init_tally_routines() + subroutine init_tally_routines() bind(C) if (run_CE) then score_general => score_general_ce score_analog_tally => score_analog_tally_ce @@ -3769,7 +3769,7 @@ contains ! within the batch to a new random variable !=============================================================================== - subroutine accumulate_tallies() + subroutine accumulate_tallies() bind(C) integer :: i real(C_DOUBLE) :: k_col ! Copy of batch collision estimate of keff @@ -3830,7 +3830,7 @@ contains ! SETUP_ACTIVE_TALLIES !=============================================================================== - subroutine setup_active_tallies() + subroutine setup_active_tallies() bind(C) integer :: i ! loop counter diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 481f05edc6..7f64832ec3 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -13,6 +13,15 @@ namespace openmc { +//============================================================================== +// Global variable definitions +//============================================================================== + +double global_tally_absorption; +double global_tally_collision; +double global_tally_tracklength; +double global_tally_leakage; + //============================================================================== // Non-member functions //============================================================================== diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index bc9a9bf9bc..0b80255d50 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -136,10 +136,10 @@ module tally_header ! global tally, it can cause a higher cache miss rate due to ! invalidation. Thus, we use threadprivate variables to accumulate global ! tallies and then reduce at the end of a generation. - real(C_DOUBLE), public :: global_tally_collision = ZERO - real(C_DOUBLE), public, bind(C) :: global_tally_absorption = ZERO - real(C_DOUBLE), public :: global_tally_tracklength = ZERO - real(C_DOUBLE), public :: global_tally_leakage = ZERO + real(C_DOUBLE), public, bind(C) :: global_tally_collision + real(C_DOUBLE), public, bind(C) :: global_tally_absorption + real(C_DOUBLE), public, bind(C) :: global_tally_tracklength + real(C_DOUBLE), public, bind(C) :: global_tally_leakage !$omp threadprivate(global_tally_collision, global_tally_absorption, & !$omp& global_tally_tracklength, global_tally_leakage) @@ -392,7 +392,7 @@ contains ! tallies.xml file. !=============================================================================== - subroutine configure_tallies() + subroutine configure_tallies() bind(C) integer :: i diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index fbecaecd93..f5efd04378 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -29,7 +29,7 @@ contains ! and predicts the number of remainining batches to convergence. !=============================================================================== - subroutine check_triggers() + subroutine check_triggers() bind(C) implicit none diff --git a/src/timer.cpp b/src/timer.cpp index 85b1eb7717..46a84a6f1e 100644 --- a/src/timer.cpp +++ b/src/timer.cpp @@ -38,23 +38,48 @@ double Timer::elapsed() // Global variables //============================================================================== +Timer time_active; Timer time_bank; Timer time_bank_sample; Timer time_bank_sendrecv; +Timer time_finalize; +Timer time_inactive; +Timer time_initialize; +Timer time_tallies; +Timer time_total; +Timer time_transport; //============================================================================== // Fortran compatibility //============================================================================== +extern "C" double time_active_elapsed() { return time_active.elapsed(); } extern "C" double time_bank_elapsed() { return time_bank.elapsed(); } extern "C" double time_bank_sample_elapsed() { return time_bank_sample.elapsed(); } extern "C" double time_bank_sendrecv_elapsed() { return time_bank_sendrecv.elapsed(); } +extern "C" double time_finalize_elapsed() { return time_finalize.elapsed(); } +extern "C" double time_inactive_elapsed() { return time_inactive.elapsed(); } +extern "C" double time_initialize_elapsed() { return time_initialize.elapsed(); } +extern "C" double time_tallies_elapsed() { return time_tallies.elapsed(); } +extern "C" double time_total_elapsed() { return time_total.elapsed(); } +extern "C" double time_transport_elapsed() { return time_transport.elapsed(); } -extern "C" void reset_timers() +//============================================================================== +// Non-member functions +//============================================================================== + +void reset_timers() { + time_active.reset(); time_bank.reset(); time_bank_sample.reset(); time_bank_sendrecv.reset(); + time_finalize.reset(); + time_inactive.reset(); + time_initialize.reset(); + time_tallies.reset(); + time_total.reset(); + time_transport.reset(); } } // namespace openmc diff --git a/src/timer_header.F90 b/src/timer_header.F90 index 1691546964..2064f7ae66 100644 --- a/src/timer_header.F90 +++ b/src/timer_header.F90 @@ -7,6 +7,10 @@ module timer_header implicit none interface + function time_active_elapsed() result(t) bind(C) + import C_DOUBLE + real(C_DOUBLE) :: t + end function function time_bank_elapsed() result(t) bind(C) import C_DOUBLE real(C_DOUBLE) :: t @@ -19,6 +23,30 @@ module timer_header import C_DOUBLE real(C_DOUBLE) :: t end function + function time_finalize_elapsed() result(t) bind(C) + import C_DOUBLE + real(C_DOUBLE) :: t + end function + function time_inactive_elapsed() result(t) bind(C) + import C_DOUBLE + real(C_DOUBLE) :: t + end function + function time_initialize_elapsed() result(t) bind(C) + import C_DOUBLE + real(C_DOUBLE) :: t + end function + function time_tallies_elapsed() result(t) bind(C) + import C_DOUBLE + real(C_DOUBLE) :: t + end function + function time_total_elapsed() result(t) bind(C) + import C_DOUBLE + real(C_DOUBLE) :: t + end function + function time_transport_elapsed() result(t) bind(C) + import C_DOUBLE + real(C_DOUBLE) :: t + end function subroutine reset_timers() bind(C) end subroutine end interface @@ -44,15 +72,7 @@ module timer_header ! ============================================================================ ! TIMING VARIABLES - type(Timer) :: time_total ! timer for total run - type(Timer) :: time_initialize ! timer for initialization type(Timer) :: time_read_xs ! timer for reading cross sections - type(Timer) :: time_unionize ! timer for material xs-energy grid union - type(Timer) :: time_tallies ! timer for accumulate tallies - type(Timer) :: time_inactive ! timer for inactive batches - type(Timer) :: time_active ! timer for active batches - type(Timer) :: time_transport ! timer for transport only - type(Timer) :: time_finalize ! timer for finalization contains @@ -117,4 +137,12 @@ contains self % elapsed = ZERO end subroutine timer_reset +!=============================================================================== +! RESET_TIMERS resets timers on the Fortran side +!=============================================================================== + + subroutine reset_timers_f() bind(C) + call time_read_xs % reset() + end subroutine + end module timer_header diff --git a/src/tracking.F90 b/src/tracking.F90 index dc0c008d61..038c7fdb54 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -48,7 +48,7 @@ contains ! TRANSPORT encompasses the main logic for moving a particle through geometry. !=============================================================================== - subroutine transport(p) + subroutine transport(p) bind(C) type(Particle), intent(inout) :: p @@ -287,8 +287,7 @@ contains ! Check for secondary particles if this particle is dead if (.not. p % alive) then if (p % n_secondary > 0) then - call particle_from_source(p, p % secondary_bank(p % n_secondary), & - run_CE, energy_bin_avg) + call particle_from_source(p, p % secondary_bank(p % n_secondary)) p % n_secondary = p % n_secondary - 1 n_event = 0