From 05ecbb153af4639b0aa669797b18e61c930c96a6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 16 Oct 2018 13:40:12 -0500 Subject: [PATCH] Convert calculate_average_keff to C++ --- include/openmc/eigenvalue.h | 9 ++++++ include/openmc/state_point.h | 2 +- src/api.F90 | 8 +++-- src/eigenvalue.F90 | 61 ++---------------------------------- src/eigenvalue.cpp | 48 ++++++++++++++++++++++++++++ src/simulation.F90 | 3 +- src/state_point.F90 | 18 +++-------- src/state_point.cpp | 15 +++++++++ 8 files changed, 88 insertions(+), 76 deletions(-) diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index cbf0baa2d..c2ecca30e 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -1,6 +1,7 @@ #ifndef OPENMC_EIGENVALUE_H #define OPENMC_EIGENVALUE_H +#include #include // for int64_t #include @@ -15,6 +16,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 +30,13 @@ extern "C" int64_t n_bank; //! Collect/normalize the tracklength keff from each process extern "C" void calculate_generation_keff(); +//! Calcaulte 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(); + //! Sample/redistribute source sites from accumulated fission sites extern "C" void synchronize_bank(); diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index 058edc09c..8f98fccd8 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -11,8 +11,8 @@ namespace openmc { 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/src/api.F90 b/src/api.F90 index 18162e600..4a1b7313e 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -4,7 +4,7 @@ module openmc_api use bank_header, only: openmc_source_bank use constants, only: K_BOLTZMANN - use eigenvalue, only: k_sum, openmc_get_keff + use eigenvalue, only: openmc_get_keff use error use geometry, only: find_cell use geometry_header @@ -248,6 +248,10 @@ contains integer(C_INT) :: err integer :: i + interface + subroutine k_sum_reset() bind(C) + end subroutine + end interface if (allocated(tallies)) then do i = 1, size(tallies) @@ -268,7 +272,7 @@ contains k_col_abs = ZERO k_col_tra = ZERO k_abs_tra = ZERO - k_sum(:) = ZERO + call k_sum_reset() ! Reset timers call time_total % reset() diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 92c731eae..ed9dce3c6 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -2,79 +2,24 @@ 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 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 + + subroutine calculate_average_keff() bind(C) + end subroutine 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 diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index a9f99a946..0dec2066b 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -8,6 +8,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 +19,7 @@ #include "openmc/tallies/tally.h" #include // for min +#include // for sqrt #include namespace openmc { @@ -27,6 +29,7 @@ namespace openmc { //============================================================================== double keff_generation; +std::array k_sum; std::vector entropy; xt::xtensor source_frac; @@ -298,6 +301,46 @@ void synchronize_bank() 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)); + } + } +} + void shannon_entropy() { // Get pointer to entropy mesh @@ -430,6 +473,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 +487,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/simulation.F90 b/src/simulation.F90 index 0d6c534f1..e3f8cd3b4 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -10,8 +10,7 @@ module simulation 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 + use eigenvalue, only: calculate_average_keff, calculate_generation_keff #ifdef _OPENMP use eigenvalue, only: join_bank_from_threads #endif diff --git a/src/state_point.F90 b/src/state_point.F90 index 52b01557d..4943ba134 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -16,7 +16,7 @@ module state_point use bank_header, only: Bank use cmfd_header use constants - use eigenvalue, only: openmc_get_keff, k_sum + use eigenvalue, only: openmc_get_keff use endf, only: reaction_name use error, only: fatal_error, warning, write_message use hdf5_interface @@ -488,7 +488,6 @@ contains subroutine load_state_point() integer :: i - integer :: n integer :: int_array(3) integer, allocatable :: array(:) integer(C_INT64_T) :: seed @@ -504,6 +503,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 +614,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 63c8264a3..ed679836b 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -9,6 +9,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" @@ -275,4 +276,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