Extract the combined k-effective estimator and fix its two-estimate branch

include/openmc/math_functions.h

@@ -8,8 +8,10 @@
 #include <complex>
 #include <cstdlib>
 
+#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<double, 3>& estimates,
+  const tensor::StaticTensor2D<double, 3, 3>& cov, int64_t n,
+  array<double, 2>& combined);
+
 } // namespace openmc
 #endif // OPENMC_MATH_FUNCTIONS_H

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<double>::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<double>::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<double, 3> kv {};
-  tensor::Tensor<double> cov = tensor::zeros<double>({3, 3});
+  tensor::StaticTensor2D<double, 3, 3> 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<double, 2> 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<double, 3> 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;
 }
 

src/math_functions.cpp

@@ -1,10 +1,14 @@
 #include "openmc/math_functions.h"
 
+#include <cmath>  // for abs, sqrt
 #include <limits> // for numeric_limits
+#include <string> // 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<double, 3>& estimates,
+  const tensor::StaticTensor2D<double, 3, 3>& cov, int64_t n,
+  array<double, 2>& 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<double, 3> 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

tests/cpp_unit_tests/CMakeLists.txt

@@ -1,4 +1,5 @@
 set(TEST_NAMES
+  test_combine_estimates
   test_distribution
   test_file_utils
   test_tally