From 67a4a592a8cf453a210fece001702116ac59c698 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 8 Sep 2026 02:53:18 +0300 Subject: [PATCH 1/5] refactor out the function to combine estimators --- include/openmc/math_functions.h | 33 +++ src/eigenvalue.cpp | 205 +++++------------- src/math_functions.cpp | 132 +++++++++++ tests/cpp_unit_tests/CMakeLists.txt | 1 + .../cpp_unit_tests/test_combine_estimates.cpp | 136 ++++++++++++ 5 files changed, 360 insertions(+), 147 deletions(-) create mode 100644 tests/cpp_unit_tests/test_combine_estimates.cpp diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index d3bccca7b5d..b75718a1111 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,36 @@ 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. +//! +//! \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 +//! \return Whether a combination was formed. False when there are too few +//! realizations, or when the covariance is degenerate enough that the result +//! is not finite; in either case the caller must supply its own estimate. +//============================================================================== + +bool 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..edb64dee152 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -427,160 +427,71 @@ void calculate_average_keff() int openmc_get_keff(double* k_combined) { - k_combined[0] = 0.0; - k_combined[1] = 0.0; - - // 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) { - k_combined[0] = simulation::keff; - k_combined[1] = simulation::keff_std; - if (simulation::n_realizations <= 1) { - k_combined[1] = 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}); - 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; - cov(0, 0) = - (gt(GlobalTally::K_COLLISION, TallyResult::SUM_SQ) - n * kv[0] * kv[0]) / - (n - 1); - cov(1, 1) = - (gt(GlobalTally::K_ABSORPTION, TallyResult::SUM_SQ) - n * kv[1] * kv[1]) / - (n - 1); - cov(2, 2) = - (gt(GlobalTally::K_TRACKLENGTH, TallyResult::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; - 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; - - } else { - // No two estimators match, so set boolean to use all three estimators. - use_three = true; + // The covariance needs at least two realizations to be formed at all; a + // combination needs more still, which combine_estimates() decides. Random + // ray produces a single estimate of k rather than three independent ones, + // so there is nothing to combine there either. + bool combined = false; + if (n > 1 && settings::solver_type != SolverType::RANDOM_RAY) { + // Copy estimates of k-effective and its variance (not variance of the + // mean) + const auto& gt = simulation::global_tallies; + + array kv {}; + 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; + cov(0, 0) = + (gt(GlobalTally::K_COLLISION, TallyResult::SUM_SQ) - n * kv[0] * kv[0]) / + (n - 1); + cov(1, 1) = (gt(GlobalTally::K_ABSORPTION, TallyResult::SUM_SQ) - + n * kv[1] * kv[1]) / + (n - 1); + cov(2, 2) = (gt(GlobalTally::K_TRACKLENGTH, TallyResult::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); + + // 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 + array result; + combined = combine_estimates(kv, cov, n, result); + if (combined) { + k_combined[0] = result[0]; + k_combined[1] = result[1]; + } } - 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; - } + if (!combined) { + // Report the average over generations. This function has to return a value + // whenever it is called -- including from a statepoint written during the + // inactive batches, when no realization has been accumulated -- and the + // generation average is the only estimate of k that is 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 fallback at all, but the only estimate there is. + k_combined[0] = simulation::keff; + k_combined[1] = simulation::keff_std; - // Complete calculations of S sums - for (auto& S_i : S) { - S_i *= (n - 1); + // 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 + if (n <= 1) { + k_combined[1] = std::numeric_limits::infinity(); } - 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..efe15723338 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -1,9 +1,11 @@ #include "openmc/math_functions.h" +#include // for abs, sqrt, isfinite #include // for numeric_limits #include "openmc/external/Faddeeva.hh" +#include "openmc/array.h" #include "openmc/constants.h" #include "openmc/random_lcg.h" @@ -1009,4 +1011,134 @@ 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); } +bool combine_estimates(const array& estimates, + const tensor::StaticTensor2D& cov, int64_t n, + array& combined) +{ + combined[0] = 0.0; + combined[1] = 0.0; + + // The three-estimate expression has an n-3 term in a denominator, and the + // two-estimate expression an n-2 term + if (n <= 3) + return false; + + // 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]); + } + + return std::isfinite(combined[0]) && std::isfinite(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..8fa87b0632b --- /dev/null +++ b/tests/cpp_unit_tests/test_combine_estimates.cpp @@ -0,0 +1,136 @@ +#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; + REQUIRE(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; + REQUIRE(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; + REQUIRE(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; + REQUIRE(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; + REQUIRE(combine_estimates(estimates, cov, 100, low)); + REQUIRE(combine_estimates(estimates, cov, 400, high)); + REQUIRE_THAT(low[1] / high[1], WithinRel(2.0, 0.02)); +} + +TEST_CASE("Test combine_estimates rejects too few realizations") +{ + auto cov = distinct_cov(); + array result; + + // The three-estimate expression has an n-3 term in a denominator and the + // two-estimate expression an n-2 term, so a combination is only defined + // above three realizations. Callers supply their own estimate below that. + for (int64_t n : {int64_t(0), int64_t(1), int64_t(2), int64_t(3)}) { + REQUIRE_FALSE(combine_estimates({0.980, 0.982, 0.981}, cov, n, result)); + } + REQUIRE(combine_estimates({0.980, 0.982, 0.981}, cov, 4, result)); +} From 92481daf9500bcbd15cd4a38e69b22f9419f713d Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 8 Sep 2026 03:01:41 +0300 Subject: [PATCH 2/5] clang format --- src/eigenvalue.cpp | 6 +++--- tests/cpp_unit_tests/test_combine_estimates.cpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index edb64dee152..e73dcccf19b 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -448,9 +448,9 @@ int openmc_get_keff(double* k_combined) cov(0, 0) = (gt(GlobalTally::K_COLLISION, TallyResult::SUM_SQ) - n * kv[0] * kv[0]) / (n - 1); - cov(1, 1) = (gt(GlobalTally::K_ABSORPTION, TallyResult::SUM_SQ) - - n * kv[1] * kv[1]) / - (n - 1); + cov(1, 1) = + (gt(GlobalTally::K_ABSORPTION, TallyResult::SUM_SQ) - n * kv[1] * kv[1]) / + (n - 1); cov(2, 2) = (gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM_SQ) - n * kv[2] * kv[2]) / (n - 1); diff --git a/tests/cpp_unit_tests/test_combine_estimates.cpp b/tests/cpp_unit_tests/test_combine_estimates.cpp index 8fa87b0632b..638bdcf7e6f 100644 --- a/tests/cpp_unit_tests/test_combine_estimates.cpp +++ b/tests/cpp_unit_tests/test_combine_estimates.cpp @@ -46,8 +46,8 @@ openmc::tensor::StaticTensor2D coincident_cov() //! 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) +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; From d72acaa851705e30392f714d68c83bfab652b2b7 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 8 Sep 2026 03:08:14 +0300 Subject: [PATCH 3/5] updated mg_survival_biasing results --- tests/regression_tests/mg_survival_biasing/results_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 239d78ca04c87993d0621b41fcc35409042cf5e4 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 8 Sep 2026 04:30:06 +0300 Subject: [PATCH 4/5] fix a regression problem --- include/openmc/math_functions.h | 5 ++--- src/math_functions.cpp | 4 ++-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index b75718a1111..72ccef2acc4 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -283,9 +283,8 @@ bool isclose(double a, double b, double rel_tol, double abs_tol); //! 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 -//! \return Whether a combination was formed. False when there are too few -//! realizations, or when the covariance is degenerate enough that the result -//! is not finite; in either case the caller must supply its own estimate. +//! \return Whether there were enough realizations to form a combination at all. +//! When false, the caller must supply its own estimate. //============================================================================== bool combine_estimates(const array& estimates, diff --git a/src/math_functions.cpp b/src/math_functions.cpp index efe15723338..797a84ebd6e 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -1,6 +1,6 @@ #include "openmc/math_functions.h" -#include // for abs, sqrt, isfinite +#include // for abs, sqrt #include // for numeric_limits #include "openmc/external/Faddeeva.hh" @@ -1138,7 +1138,7 @@ bool combine_estimates(const array& estimates, combined[1] = std::sqrt(combined[1]); } - return std::isfinite(combined[0]) && std::isfinite(combined[1]); + return true; } } // namespace openmc From c53f9baffe8cf47bf7444b9936d4f9f38607ec96 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 8 Sep 2026 04:47:11 +0300 Subject: [PATCH 5/5] simplify design a bit --- include/openmc/math_functions.h | 15 ++- src/eigenvalue.cpp | 105 ++++++++---------- src/math_functions.cpp | 15 +-- .../cpp_unit_tests/test_combine_estimates.cpp | 31 +++--- 4 files changed, 84 insertions(+), 82 deletions(-) diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index 72ccef2acc4..400f30f701b 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -278,16 +278,25 @@ bool isclose(double a, double b, double rel_tol, double abs_tol); //! 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 -//! \return Whether there were enough realizations to form a combination at all. -//! When false, the caller must supply its own estimate. //============================================================================== -bool combine_estimates(const array& estimates, +//! 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); diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index e73dcccf19b..8b5e60c8469 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -429,69 +429,60 @@ int openmc_get_keff(double* k_combined) { int64_t n = simulation::n_realizations; - // The covariance needs at least two realizations to be formed at all; a - // combination needs more still, which combine_estimates() decides. Random - // ray produces a single estimate of k rather than three independent ones, - // so there is nothing to combine there either. - bool combined = false; - if (n > 1 && settings::solver_type != SolverType::RANDOM_RAY) { - // Copy estimates of k-effective and its variance (not variance of the - // mean) - const auto& gt = simulation::global_tallies; - - array kv {}; - 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; - cov(0, 0) = - (gt(GlobalTally::K_COLLISION, TallyResult::SUM_SQ) - n * kv[0] * kv[0]) / - (n - 1); - cov(1, 1) = - (gt(GlobalTally::K_ABSORPTION, TallyResult::SUM_SQ) - n * kv[1] * kv[1]) / - (n - 1); - cov(2, 2) = (gt(GlobalTally::K_TRACKLENGTH, TallyResult::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); - - // 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 - array result; - combined = combine_estimates(kv, cov, n, result); - if (combined) { - k_combined[0] = result[0]; - k_combined[1] = result[1]; - } - } - - if (!combined) { - // Report the average over generations. This function has to return a value - // whenever it is called -- including from a statepoint written during the - // inactive batches, when no realization has been accumulated -- and the - // generation average is the only estimate of k that is 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 fallback at all, but the only estimate there is. + // 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; // 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 - if (n <= 1) { - k_combined[1] = std::numeric_limits::infinity(); - } + k_combined[1] = + n > 1 ? simulation::keff_std : std::numeric_limits::infinity(); + return 0; } + // Copy estimates of k-effective and its variance (not variance of the mean) + const auto& gt = simulation::global_tallies; + + array kv {}; + 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; + cov(0, 0) = + (gt(GlobalTally::K_COLLISION, TallyResult::SUM_SQ) - n * kv[0] * kv[0]) / + (n - 1); + cov(1, 1) = + (gt(GlobalTally::K_ABSORPTION, TallyResult::SUM_SQ) - n * kv[1] * kv[1]) / + (n - 1); + cov(2, 2) = + (gt(GlobalTally::K_TRACKLENGTH, TallyResult::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); + + // 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]; + return 0; } diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 797a84ebd6e..92562c4e70d 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -2,11 +2,13 @@ #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 { @@ -1011,17 +1013,18 @@ 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); } -bool combine_estimates(const array& estimates, +void combine_estimates(const array& estimates, const tensor::StaticTensor2D& cov, int64_t n, array& combined) { combined[0] = 0.0; combined[1] = 0.0; - // The three-estimate expression has an n-3 term in a denominator, and the - // two-estimate expression an n-2 term - if (n <= 3) - return false; + 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 @@ -1137,8 +1140,6 @@ bool combine_estimates(const array& estimates, ((n - 1) * g + n * f * f) / (n * (n - 2) * g * g); combined[1] = std::sqrt(combined[1]); } - - return true; } } // namespace openmc diff --git a/tests/cpp_unit_tests/test_combine_estimates.cpp b/tests/cpp_unit_tests/test_combine_estimates.cpp index 638bdcf7e6f..bada38f1b28 100644 --- a/tests/cpp_unit_tests/test_combine_estimates.cpp +++ b/tests/cpp_unit_tests/test_combine_estimates.cpp @@ -70,18 +70,18 @@ TEST_CASE("Test combine_estimates with three distinct estimates") int64_t n = 100; array result; - REQUIRE(combine_estimates({0.980, 0.982, 0.981}, cov, n, 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; - REQUIRE(combine_estimates({0.975, 0.975, 0.975}, cov, n, 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; - REQUIRE(combine_estimates({0.990, 0.992, 0.991}, cov, n, shifted)); + combine_estimates({0.990, 0.992, 0.991}, cov, n, shifted); REQUIRE_THAT(shifted[0] - result[0], WithinAbs(0.01, 1e-12)); } @@ -94,7 +94,7 @@ TEST_CASE("Test combine_estimates with two coincident estimates") int64_t n = 100; array result; - REQUIRE(combine_estimates(estimates, cov, n, 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); @@ -116,21 +116,22 @@ TEST_CASE("Test combine_estimates standard deviation scales as 1/sqrt(n)") // 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; - REQUIRE(combine_estimates(estimates, cov, 100, low)); - REQUIRE(combine_estimates(estimates, cov, 400, 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 rejects too few realizations") +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; - - // The three-estimate expression has an n-3 term in a denominator and the - // two-estimate expression an n-2 term, so a combination is only defined - // above three realizations. Callers supply their own estimate below that. - for (int64_t n : {int64_t(0), int64_t(1), int64_t(2), int64_t(3)}) { - REQUIRE_FALSE(combine_estimates({0.980, 0.982, 0.981}, cov, n, result)); - } - REQUIRE(combine_estimates({0.980, 0.982, 0.981}, cov, 4, result)); + combine_estimates( + {0.980, 0.982, 0.981}, cov, openmc::MIN_REALIZATIONS_TO_COMBINE, result); + REQUIRE(result[1] > 0.0); }