mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Convert openmc_get_keff to C++
This commit is contained in:
Paul Romano
parent
05ecbb153a
commit
1c7ebf198e
3 changed files with 171 additions and 164 deletions
|
|
@ -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();
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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 <algorithm> // for min
|
||||
#include <cmath> // for sqrt
|
||||
#include <array>
|
||||
#include <cmath> // for sqrt, abs, pow
|
||||
#include <string>
|
||||
|
||||
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<double, 3> kv {};
|
||||
xt::xtensor<double, 2> cov = xt::zeros<double>({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<double, 3> 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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue