From 1c7ebf198eb8eff4ff891bbdd0513ed09470b685 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 16 Oct 2018 14:32:37 -0500 Subject: [PATCH] Convert openmc_get_keff to C++ --- include/openmc/eigenvalue.h | 15 ++++ src/eigenvalue.F90 | 169 ++---------------------------------- src/eigenvalue.cpp | 151 +++++++++++++++++++++++++++++++- 3 files changed, 171 insertions(+), 164 deletions(-) diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index c2ecca30e..ab1c82265 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -37,6 +37,21 @@ extern "C" void calculate_generation_keff(); //! 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/src/eigenvalue.F90 b/src/eigenvalue.F90 index ed9dce3c6..fbf4c65ae 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -2,7 +2,6 @@ module eigenvalue use, intrinsic :: ISO_C_BINDING - use constants, only: ZERO use message_passing use settings use simulation_header @@ -16,172 +15,16 @@ module eigenvalue subroutine calculate_average_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 -!=============================================================================== -! 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 diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 0dec2066b..c60961850 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" @@ -19,7 +20,8 @@ #include "openmc/tallies/tally.h" #include // for min -#include // for sqrt +#include +#include // for sqrt, abs, pow #include namespace openmc { @@ -341,6 +343,153 @@ void calculate_average_keff() } } +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