diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index d3bccca7b5d..400f30f701b 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -8,8 +8,10 @@ #include #include +#include "openmc/array.h" #include "openmc/position.h" #include "openmc/search.h" +#include "openmc/tensor.h" namespace openmc { @@ -259,5 +261,44 @@ double standard_normal_cdf(double z); //============================================================================== bool isclose(double a, double b, double rel_tol, double abs_tol); +//============================================================================== +//! Combine three correlated estimates of the same quantity +//! +//! Returns the linear combination of the estimates, with weights summing to +//! one, that has the smallest variance. 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), whose +//! expression for the standard deviation accounts for the weights having been +//! estimated from the same realizations as the estimates themselves. +//! +//! If two of the estimates coincide the three-estimate expression is singular, +//! and an expression derived for a combination of two estimates is used +//! instead. +//! +//! \p n must be at least MIN_REALIZATIONS_TO_COMBINE; below that the +//! covariance is singular and the expressions are undefined. +//! +//! \param[in] estimates The three estimates +//! \param[in] cov Covariance of the three estimates over a single +//! realization, not of the mean +//! \param[in] n Number of realizations each estimate was formed from +//! \param[out] combined The combination and the standard deviation of its mean +//============================================================================== + +//! Fewest realizations from which a combination can be formed +//! +//! A k by k sample covariance built from n realizations has rank at most +//! n - 1, so it is singular unless n exceeds k. The n-3 and n-2 factors in the +//! expressions for the standard deviation are the residual degrees of freedom +//! and vanish at the same point. +constexpr int64_t MIN_REALIZATIONS_TO_COMBINE {4}; + +void combine_estimates(const array& estimates, + const tensor::StaticTensor2D& cov, int64_t n, + array& combined); + } // namespace openmc #endif // OPENMC_MATH_FUNCTIONS_H diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 71ae1707ca6..8b5e60c8469 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -427,29 +427,32 @@ void calculate_average_keff() int openmc_get_keff(double* k_combined) { - k_combined[0] = 0.0; - k_combined[1] = 0.0; + int64_t n = simulation::n_realizations; - // Special case for n <=3. Notice that at the end, - // there is a N-3 term in a denominator. - if (simulation::n_realizations <= 3 || - settings::solver_type == SolverType::RANDOM_RAY) { + // Random ray computes a single estimate of k from the scalar flux rather + // than three independent ones, and a combination is not defined below + // MIN_REALIZATIONS_TO_COMBINE realizations. In both cases report the average + // over generations, which is the only estimate of k defined at every point + // in a run: during inactive generations it holds the most recent generation + // estimate, and thereafter the average over active ones. For random ray it + // is not a substitute at all, but the only estimate there is. + if (settings::solver_type == SolverType::RANDOM_RAY || + n < MIN_REALIZATIONS_TO_COMBINE) { k_combined[0] = simulation::keff; - k_combined[1] = simulation::keff_std; - if (simulation::n_realizations <= 1) { - k_combined[1] = std::numeric_limits::infinity(); - } + + // keff_std is only assigned once there is more than one active generation + // to take a spread over, so it carries no meaning below that + k_combined[1] = + n > 1 ? simulation::keff_std : std::numeric_limits::infinity(); return 0; } - // Initialize variables - int64_t n = simulation::n_realizations; - // Copy estimates of k-effective and its variance (not variance of the mean) const auto& gt = simulation::global_tallies; array kv {}; - tensor::Tensor cov = tensor::zeros({3, 3}); + tensor::StaticTensor2D cov; + cov.fill(0.0); kv[0] = gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n; kv[1] = gt(GlobalTally::K_ABSORPTION, TallyResult::SUM) / n; kv[2] = gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM) / n; @@ -471,116 +474,15 @@ int openmc_get_keff(double* k_combined) 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; - bool use_three = false; - 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; + // In multi-group mode with survival biasing the collision and absorption + // estimators are identical, which combine_estimates() detects and handles + // with its two-estimate expression. Whatever it produces is reported as it + // stands, including for a degenerate covariance. + array result; + combine_estimates(kv, cov, n, result); + k_combined[0] = result[0]; + k_combined[1] = result[1]; - } else { - // No two estimators match, so set boolean to use all three estimators. - use_three = true; - } - - if (use_three) { - // Use three estimators as derived in the paper by Urbatsch - - // Initialize variables - double g = 0.0; - array S {}; - - for (int l = 0; l < 3; ++l) { - // Permutations of estimates - int k; - 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; } diff --git a/src/math_functions.cpp b/src/math_functions.cpp index ddacc2bd9b6..92562c4e70d 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -1,10 +1,14 @@ #include "openmc/math_functions.h" +#include // for abs, sqrt #include // for numeric_limits +#include // for to_string #include "openmc/external/Faddeeva.hh" +#include "openmc/array.h" #include "openmc/constants.h" +#include "openmc/error.h" #include "openmc/random_lcg.h" namespace openmc { @@ -1009,4 +1013,133 @@ bool isclose(double a, double b, double rel_tol, double abs_tol) std::max(rel_tol * std::max(std::abs(a), std::abs(b)), abs_tol); } +void combine_estimates(const array& estimates, + const tensor::StaticTensor2D& cov, int64_t n, + array& combined) +{ + combined[0] = 0.0; + combined[1] = 0.0; + + if (n < MIN_REALIZATIONS_TO_COMBINE) { + fatal_error("combine_estimates() requires at least " + + std::to_string(MIN_REALIZATIONS_TO_COMBINE) + + " realizations; the covariance is singular below that."); + } + + // Check to see if two estimates are the same. If they are, the three + // estimate expressions are singular and will produce floating-point + // exceptions, so an expression specifically derived for the combination of + // two estimates (vice three) is used instead. + + // First we will identify if there are any matching estimates + int i, j; + bool use_three = false; + if ((std::abs(estimates[0] - estimates[1]) / estimates[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(estimates[0] - estimates[2]) / estimates[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(estimates[1] - estimates[2]) / estimates[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 estimates match, so set boolean to use all three estimates. + use_three = true; + } + + if (use_three) { + // Use three estimates as derived in the paper by Urbatsch + + // Initialize variables + double g = 0.0; + array S {}; + + for (int l = 0; l < 3; ++l) { + // Permutations of the three estimates + int k; + switch (l) { + case 0: + i = 0; + j = 1; + k = 2; + break; + case 1: + i = 1; + j = 2; + k = 0; + break; + case 2: + 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)) * estimates[l] * estimates[l]; + S[2] += (cov(k, k) + cov(i, j) - cov(j, k) - cov(i, k)) * estimates[l] * + estimates[j]; + + // Add to sum for the combination + combined[0] += f * estimates[l]; + g += f; + } + + // Complete calculations of S sums + for (auto& S_i : S) { + S_i *= (n - 1); + } + S[0] *= (n - 1) * (n - 1); + + // Calculate the combination + combined[0] /= g; + + // Calculate standard deviation of the combination + g *= (n - 1) * (n - 1); + combined[1] = + std::sqrt(S[0] / (g * n * (n - 3)) * (1 + n * ((S[1] - 2 * S[2]) / g))); + + } else { + // Use only two estimates + // These equations are derived analogously to that done in the paper by + // Urbatsch, but are simpler than for the three estimate case since the + // block matrices of the three estimate equations reduces to scalars here + + // Store the commonly used term + double f = estimates[i] - estimates[j]; + double g = cov(i, i) + cov(j, j) - 2.0 * cov(i, j); + + // Calculate the combination + combined[0] = estimates[i] - (cov(i, i) - cov(i, j)) / g * f; + + // Calculate standard deviation of the combination. Urbatsch's Eq. 40 is + // written in terms of the matrix S rather than the sample covariance + // Sigma = S / (n - 1). The factor cancels in the combination itself but + // not here, and omitting it understates the standard deviation by up to + // sqrt(n - 1). + combined[1] = (cov(i, i) * cov(j, j) - cov(i, j) * cov(i, j)) * + ((n - 1) * g + n * f * f) / (n * (n - 2) * g * g); + combined[1] = std::sqrt(combined[1]); + } +} + } // namespace openmc diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index 991f219f528..e61f8d7eb75 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -1,4 +1,5 @@ set(TEST_NAMES + test_combine_estimates test_distribution test_file_utils test_tally diff --git a/tests/cpp_unit_tests/test_combine_estimates.cpp b/tests/cpp_unit_tests/test_combine_estimates.cpp new file mode 100644 index 00000000000..bada38f1b28 --- /dev/null +++ b/tests/cpp_unit_tests/test_combine_estimates.cpp @@ -0,0 +1,137 @@ +#include + +#include +#include + +#include "openmc/array.h" +#include "openmc/math_functions.h" +#include "openmc/tensor.h" + +using Catch::Matchers::WithinAbs; +using Catch::Matchers::WithinRel; +using openmc::array; +using openmc::combine_estimates; + +namespace { + +//! Covariance in which no two diagonal entries match, so that the +//! three-estimate expression is used +openmc::tensor::StaticTensor2D distinct_cov() +{ + openmc::tensor::StaticTensor2D cov; + cov.fill(0.0); + cov(0, 0) = 9.0e-8; + cov(1, 1) = 1.6e-7; + cov(2, 2) = 4.0e-8; + cov(0, 1) = cov(1, 0) = 1.1e-7; + cov(0, 2) = cov(2, 0) = 5.5e-8; + cov(1, 2) = cov(2, 1) = 7.5e-8; + return cov; +} + +//! Covariance in which the first two estimators are identical, as in +//! multi-group mode with survival biasing +openmc::tensor::StaticTensor2D coincident_cov() +{ + openmc::tensor::StaticTensor2D cov; + cov.fill(0.0); + cov(0, 0) = cov(1, 1) = cov(0, 1) = cov(1, 0) = 4.0e-8; + cov(2, 2) = 4.1e-8; + cov(0, 2) = cov(2, 0) = cov(1, 2) = cov(2, 1) = 3.9e-8; + return cov; +} + +//! Combination of two estimates, transcribed from Eq. 36 and Eq. 40 of +//! Urbatsch's LA-12658-MS. Those equations are written in terms of the matrix +//! S, so the sample covariance is converted with S = (n - 1) * Sigma before +//! being substituted. This is an independent statement of the same result and +//! is what pins down the standard deviation. +array two_estimate_reference(double e_i, double e_j, double sigma_ii, + double sigma_jj, double sigma_ij, int64_t n) +{ + double s_ii = (n - 1) * sigma_ii; + double s_jj = (n - 1) * sigma_jj; + double s_ij = (n - 1) * sigma_ij; + + double f = e_i - e_j; + double g = s_ii + s_jj - 2.0 * s_ij; + + double mean = e_i - (s_ii - s_ij) / g * f; + double variance = + (s_ii * s_jj - s_ij * s_ij) * (g + n * f * f) / (n * (n - 2) * g * g); + return {mean, std::sqrt(variance)}; +} + +} // namespace + +TEST_CASE("Test combine_estimates with three distinct estimates") +{ + auto cov = distinct_cov(); + int64_t n = 100; + + array result; + combine_estimates({0.980, 0.982, 0.981}, cov, n, result); + REQUIRE(result[1] > 0.0); + + // The weights sum to one, so estimates that all agree must combine to + // exactly that value + array agreed; + combine_estimates({0.975, 0.975, 0.975}, cov, n, agreed); + REQUIRE_THAT(agreed[0], WithinRel(0.975, 1e-12)); + + // and shifting every estimate must shift the combination by the same amount + array shifted; + combine_estimates({0.990, 0.992, 0.991}, cov, n, shifted); + REQUIRE_THAT(shifted[0] - result[0], WithinAbs(0.01, 1e-12)); +} + +TEST_CASE("Test combine_estimates with two coincident estimates") +{ + // The first two estimates match, so the three-estimate expression is + // singular and the two-estimate expression must be used instead + array estimates {0.950000, 0.950000, 0.950001}; + auto cov = coincident_cov(); + int64_t n = 100; + + array result; + combine_estimates(estimates, cov, n, result); + + auto reference = two_estimate_reference( + estimates[0], estimates[2], cov(0, 0), cov(2, 2), cov(0, 2), n); + REQUIRE_THAT(result[0], WithinRel(reference[0], 1e-12)); + REQUIRE_THAT(result[1], WithinRel(reference[1], 1e-12)); + + // With estimates this close the combination should be no worse than the + // best of them + REQUIRE(result[1] <= std::sqrt(cov(2, 2) / n)); +} + +TEST_CASE("Test combine_estimates standard deviation scales as 1/sqrt(n)") +{ + array estimates {0.950000, 0.950000, 0.950001}; + auto cov = coincident_cov(); + + // Holding the per-realization covariance fixed, the standard deviation of + // the mean falls as 1/sqrt(n). This is what fails if the two-estimate + // expression is evaluated with the sample covariance where the derivation + // calls for S, since the neglected factor carries its own dependence on n. + array low, high; + combine_estimates(estimates, cov, 100, low); + combine_estimates(estimates, cov, 400, high); + REQUIRE_THAT(low[1] / high[1], WithinRel(2.0, 0.02)); +} + +TEST_CASE("Test combine_estimates precondition") +{ + // A k by k sample covariance from n realizations has rank at most n - 1, so + // it is singular unless n exceeds k. Combining three estimates therefore + // needs four realizations, which is also where the n-3 factor in the + // standard deviation stops being positive. + REQUIRE(openmc::MIN_REALIZATIONS_TO_COMBINE == 4); + + auto cov = distinct_cov(); + array result; + combine_estimates( + {0.980, 0.982, 0.981}, cov, openmc::MIN_REALIZATIONS_TO_COMBINE, result); + REQUIRE(result[1] > 0.0); +} diff --git a/tests/regression_tests/mg_survival_biasing/results_true.dat b/tests/regression_tests/mg_survival_biasing/results_true.dat index 4b26978ada6..4f14ad167a1 100644 --- a/tests/regression_tests/mg_survival_biasing/results_true.dat +++ b/tests/regression_tests/mg_survival_biasing/results_true.dat @@ -1,2 +1,2 @@ k-combined: -9.889968E-01 9.144186E-03 +9.889968E-01 1.823442E-02