From 8ebd26b160cca2b9f0382b41f0a4ef188aa4aae6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 12 Oct 2018 08:01:23 -0500 Subject: [PATCH] Move code related to k_generation to C++ --- include/openmc/constants.h | 8 +++---- include/openmc/eigenvalue.h | 4 ++++ include/openmc/simulation.h | 3 +++ src/eigenvalue.F90 | 43 +++++++---------------------------- src/eigenvalue.cpp | 45 +++++++++++++++++++++++++++++++++++-- src/input_xml.F90 | 2 +- src/output.F90 | 4 ++-- src/simulation.F90 | 31 +++++-------------------- src/simulation.cpp | 26 ++++++++++++++++++++- src/simulation_header.F90 | 23 ++++++++++++++++--- src/state_point.F90 | 39 +++++++++----------------------- 11 files changed, 126 insertions(+), 102 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 6ead514d7..361cfc158 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -408,10 +408,10 @@ constexpr int STANDARD_DEVIATION {3}; // Global tally parameters constexpr int N_GLOBAL_TALLIES {4}; -constexpr int K_COLLISION {1}; -constexpr int K_ABSORPTION {2}; -constexpr int K_TRACKLENGTH {3}; -constexpr int LEAKAGE {4}; +constexpr int K_COLLISION {0}; +constexpr int K_ABSORPTION {1}; +constexpr int K_TRACKLENGTH {2}; +constexpr int LEAKAGE {3}; // Differential tally independent variables constexpr int DIFF_DENSITY {1}; diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index 59b448e29..cbf0baa2d 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -14,6 +14,7 @@ namespace openmc { // Global variables //============================================================================== +extern double keff_generation; //!< Single-generation k on each processor extern std::vector entropy; //!< Shannon entropy at each generation extern xt::xtensor source_frac; //!< Source fraction for UFS @@ -24,6 +25,9 @@ extern "C" int64_t n_bank; // Non-member functions //============================================================================== +//! Collect/normalize the tracklength keff from each process +extern "C" void calculate_generation_keff(); + //! Sample/redistribute source sites from accumulated fission sites extern "C" void synchronize_bank(); diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index ceb5eb92a..e21ea05d5 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -56,6 +56,9 @@ void calculate_work(); //! Initialize simulation extern "C" void openmc_simulation_init_c(); +//! Initialize a fission generation +extern "C" void initialize_generation(); + //! Determine overall generation number extern "C" int overall_generation(); diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index df108f6d7..92c731eae 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -16,42 +16,15 @@ module eigenvalue implicit none - real(8) :: keff_generation ! Single-generation k on each processor real(8) :: k_sum(2) = ZERO ! Used to reduce sum and sum_sq + interface + subroutine calculate_generation_keff() bind(C) + end subroutine + end interface + contains -!=============================================================================== -! CALCULATE_GENERATION_KEFF collects the single-processor tracklength k's onto -! the master processor and normalizes them. This should work whether or not the -! no-reduce method is being used. -!=============================================================================== - - subroutine calculate_generation_keff() - - real(8) :: keff_reduced -#ifdef OPENMC_MPI - integer :: mpi_err ! MPI error code -#endif - - ! Get keff for this generation by subtracting off the starting value - keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - keff_generation - -#ifdef OPENMC_MPI - ! Combine values across all processors - call MPI_ALLREDUCE(keff_generation, keff_reduced, 1, MPI_REAL8, & - MPI_SUM, mpi_intracomm, mpi_err) -#else - keff_reduced = keff_generation -#endif - - ! Normalize single batch estimate of k - ! TODO: This should be normalized by total_weight, not by n_particles - keff_reduced = keff_reduced / n_particles - call k_generation % push_back(keff_reduced) - - end subroutine calculate_generation_keff - !=============================================================================== ! CALCULATE_AVERAGE_KEFF calculates the mean and standard deviation of the mean ! of k-effective during active generations and broadcasts the mean to all @@ -76,12 +49,12 @@ contains if (n <= 0) then ! For inactive generations, use current generation k as estimate for next ! generation - keff = k_generation % data(i) + keff = k_generation(i) else ! Sample mean of keff - k_sum(1) = k_sum(1) + k_generation % data(i) - k_sum(2) = k_sum(2) + k_generation % data(i)**2 + 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 diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index ba3ebfae5..a9f99a946 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -5,6 +5,7 @@ #include "xtensor/xview.hpp" #include "openmc/capi.h" +#include "openmc/constants.h" #include "openmc/error.h" #include "openmc/hdf5_interface.h" #include "openmc/mesh.h" @@ -14,6 +15,7 @@ #include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/timer.h" +#include "openmc/tallies/tally.h" #include // for min #include @@ -24,6 +26,7 @@ namespace openmc { // Global variables //============================================================================== +double keff_generation; std::vector entropy; xt::xtensor source_frac; @@ -31,6 +34,28 @@ xt::xtensor source_frac; // Non-member functions //============================================================================== +void calculate_generation_keff() +{ + auto gt = global_tallies(); + + // Get keff for this generation by subtracting off the starting value + keff_generation = gt(K_TRACKLENGTH, RESULT_VALUE) - keff_generation; + + double keff_reduced; +#ifdef OPENMC_MPI + // Combine values across all processors + MPI_Allreduce(&keff_generation, &keff_reduced, 1, MPI_DOUBLE, + MPI_SUM, mpi::intracomm); +#else + keff_reduced = keff_generation; +#endif + + // Normalize single batch estimate of k + // TODO: This should be normalized by total_weight, not by n_particles + keff_reduced /= settings::n_particles; + simulation::k_generation.push_back(keff_reduced); +} + void synchronize_bank() { time_bank.start(); @@ -375,18 +400,34 @@ double ufs_get_weight(const Particle* p) } } -extern "C" void entropy_to_hdf5(hid_t group) +extern "C" void write_eigenvalue_hdf5(hid_t group) { + write_dataset(group, "n_inactive", settings::n_inactive); + write_dataset(group, "generations_per_batch", settings::gen_per_batch); + write_dataset(group, "k_generation", simulation::k_generation); if (settings::entropy_on) { write_dataset(group, "entropy", entropy); } + write_dataset(group, "k_col_abs", simulation::k_col_abs); + write_dataset(group, "k_col_tra", simulation::k_col_tra); + write_dataset(group, "k_abs_tra", simulation::k_abs_tra); + std::array k_combined; + openmc_get_keff(k_combined.data()); + write_dataset(group, "k_combined", k_combined); } -extern "C" void entropy_from_hdf5(hid_t group) +extern "C" void read_eigenvalue_hdf5(hid_t group) { + read_dataset(group, "generations_per_batch", settings::gen_per_batch); + int n = simulation::restart_batch*settings::gen_per_batch; + simulation::k_generation.resize(n); + read_dataset(group, "k_generation", simulation::k_generation); if (settings::entropy_on) { read_dataset(group, "entropy", entropy); } + read_dataset(group, "k_col_abs", simulation::k_col_abs); + read_dataset(group, "k_col_tra", simulation::k_col_tra); + read_dataset(group, "k_abs_tra", simulation::k_abs_tra); } extern "C" double entropy_c(int i) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 989f7b357..81f323351 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -223,7 +223,7 @@ contains if (run_mode == MODE_EIGENVALUE) then ! Preallocate space for keff and entropy by generation - call k_generation % reserve(n_max_batches*gen_per_batch) + call k_generation_reserve(n_max_batches*gen_per_batch) end if ! Particle tracks diff --git a/src/output.F90 b/src/output.F90 index 051167511..d775a2da7 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -339,7 +339,7 @@ contains ! write out information about batch and generation write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') & trim(to_str(current_batch)) // "/" // trim(to_str(current_gen)) - write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation % data(i) + write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation(i) ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & @@ -373,7 +373,7 @@ contains write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') & trim(to_str(current_batch)) // "/" // trim(to_str(gen_per_batch)) write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') & - k_generation % data(i) + k_generation(i) ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & diff --git a/src/simulation.F90 b/src/simulation.F90 index 4bbf94ed1..b94166f87 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -11,7 +11,7 @@ module simulation use cmfd_header, only: cmfd_on use constants, only: ZERO use eigenvalue, only: calculate_average_keff, calculate_generation_keff, & - keff_generation, k_sum + k_sum #ifdef _OPENMP use eigenvalue, only: join_bank_from_threads #endif @@ -57,6 +57,9 @@ module simulation 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 @@ -222,30 +225,6 @@ contains end subroutine initialize_batch -!=============================================================================== -! INITIALIZE_GENERATION -!=============================================================================== - - subroutine initialize_generation() - - interface - subroutine ufs_count_sites() bind(C) - end subroutine - end interface - - if (run_mode == MODE_EIGENVALUE) then - ! Reset number of fission bank sites - n_bank = 0 - - ! Count source sites if using uniform fission source weighting - if (ufs) call ufs_count_sites() - - ! Store current value of tracklength k - keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - end if - - end subroutine initialize_generation - !=============================================================================== ! FINALIZE_GENERATION !=============================================================================== @@ -438,7 +417,7 @@ contains ! 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 k_generation_clear() call entropy_clear() need_depletion_rx = .false. diff --git a/src/simulation.cpp b/src/simulation.cpp index 2935d71ac..0d94966ed 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1,6 +1,7 @@ #include "openmc/simulation.h" #include "openmc/capi.h" +#include "openmc/eigenvalue.h" #include "openmc/message_passing.h" #include "openmc/settings.h" #include "openmc/tallies/tally.h" @@ -62,7 +63,7 @@ int thread_id; //!< ID of a given thread } // namespace simulation //============================================================================== -// Functions +// Non-member functions //============================================================================== void openmc_simulation_init_c() @@ -71,6 +72,20 @@ void openmc_simulation_init_c() calculate_work(); } +void initialize_generation() +{ + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + // Reset number of fission bank sites + n_bank = 0; + + // Count source sites if using uniform fission source weighting + if (settings::ufs_on) ufs_count_sites(); + + // Store current value of tracklength k + keff_generation = global_tallies()(K_TRACKLENGTH, RESULT_VALUE); + } +} + int overall_generation() { using namespace simulation; @@ -139,4 +154,13 @@ void broadcast_triggers() } #endif +//============================================================================== +// Fortran compatibility +//============================================================================== + +extern "C" double k_generation(int i) { return simulation::k_generation.at(i - 1); } +extern "C" int k_generation_size() { return simulation::k_generation.size(); } +extern "C" void k_generation_clear() { simulation::k_generation.clear(); } +extern "C" void k_generation_reserve(int i) { simulation::k_generation.reserve(i); } + } // namespace openmc diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 29ee34ccf..65a36654f 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -39,7 +39,6 @@ module simulation_header ! K-EIGENVALUE SIMULATION VARIABLES ! Temporary k-effective values - type(VectorReal) :: k_generation ! single-generation estimates of k real(C_DOUBLE), bind(C) :: keff ! average k over active batches real(C_DOUBLE), bind(C) :: keff_std ! standard deviation of average k real(C_DOUBLE), bind(C) :: k_col_abs ! sum over batches of k_collision * k_absorption @@ -71,6 +70,25 @@ module simulation_header import C_INT integer(C_INT) :: gen end function overall_generation + + function k_generation(i) result(k) bind(C) + import C_DOUBLE, C_INT + integer(C_INT), value :: i + real(C_DOUBLE) :: k + end function + + function k_generation_size() result(sz) bind(C) + import C_INT + integer(C_INT) :: sz + end function + + 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 end interface contains @@ -83,8 +101,7 @@ contains if (allocated(work_index)) deallocate(work_index) - call k_generation % clear() - call k_generation % shrink_to_fit() + call k_generation_clear() call entropy_clear() end subroutine free_memory_simulation diff --git a/src/state_point.F90 b/src/state_point.F90 index 548db766c..f861f871a 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -63,15 +63,13 @@ contains integer(C_INT) :: err logical :: write_source_ - integer :: i, j, k + integer :: i, j integer :: i_xs integer, allocatable :: id_array(:) integer(HID_T) :: file_id integer(HID_T) :: cmfd_group, tallies_group, tally_group, & filters_group, filter_group, derivs_group, & deriv_group, runtime_group - integer(C_INT) :: ignored_err - real(C_DOUBLE) :: k_combined(2) character(MAX_WORD_LEN), allocatable :: str_array(:) character(C_CHAR), pointer :: string(:) character(len=:, kind=C_CHAR), allocatable :: filename_ @@ -83,7 +81,7 @@ contains import HID_T integer(HID_T), value :: group end subroutine - subroutine entropy_to_hdf5(group) bind(C) + subroutine write_eigenvalue_hdf5(group) bind(C) import HID_T integer(HID_T), value :: group end subroutine @@ -173,16 +171,7 @@ contains ! Write out information for eigenvalue run if (run_mode == MODE_EIGENVALUE) then - call write_dataset(file_id, "n_inactive", n_inactive) - call write_dataset(file_id, "generations_per_batch", gen_per_batch) - k = k_generation % size() - call write_dataset(file_id, "k_generation", k_generation % data(1:k)) - call entropy_to_hdf5(file_id) - call write_dataset(file_id, "k_col_abs", k_col_abs) - call write_dataset(file_id, "k_col_tra", k_col_tra) - call write_dataset(file_id, "k_abs_tra", k_abs_tra) - ignored_err = openmc_get_keff(k_combined) - call write_dataset(file_id, "k_combined", k_combined) + call write_eigenvalue_hdf5(file_id) ! Write out CMFD info if (cmfd_on) then @@ -513,7 +502,9 @@ contains character(MAX_WORD_LEN) :: word interface - subroutine entropy_from_hdf5() bind(C) + subroutine read_eigenvalue_hdf5(group) bind(C) + import HID_T + integer(HID_T), value :: group end subroutine end interface @@ -592,16 +583,7 @@ contains ! Read information specific to eigenvalue run if (run_mode == MODE_EIGENVALUE) then call read_dataset(int_array(1), file_id, "n_inactive") - call read_dataset(gen_per_batch, file_id, "generations_per_batch") - - n = restart_batch*gen_per_batch - call k_generation % resize(n) - call read_dataset(k_generation % data(1:n), file_id, "k_generation") - - call entropy_from_hdf5() - call read_dataset(k_col_abs, file_id, "k_col_abs") - call read_dataset(k_col_tra, file_id, "k_col_tra") - call read_dataset(k_abs_tra, file_id, "k_abs_tra") + call read_eigenvalue_hdf5(file_id) ! Take maximum of statepoint n_inactive and input n_inactive n_inactive = max(n_inactive, int_array(1)) @@ -634,13 +616,14 @@ contains ! 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 % data(i) - k_sum(2) = k_sum(2) + k_generation % data(i)**2 + 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 - keff = k_generation % data(n) + n = k_generation_size() + keff = k_generation(n) end if current_batch = restart_batch