From e6608f4dd42e28d0cbb531cffa2d21929512e6b2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Feb 2013 15:23:00 -0500 Subject: [PATCH 1/8] Implemented combined estimator for k-effective. --- src/finalize.F90 | 105 +++++++++++++++++++++++++++++++++++++++++++++++ src/global.F90 | 8 +++- src/output.F90 | 1 + src/tally.F90 | 13 ++++++ 4 files changed, 125 insertions(+), 2 deletions(-) diff --git a/src/finalize.F90 b/src/finalize.F90 index 1bfdfe71f4..f65f6986c7 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -32,6 +32,8 @@ contains if (run_mode /= MODE_PLOTTING) then ! Calculate statistics for tallies and write to tallies.out if (master) call tally_statistics() + if (master .and. run_mode == MODE_EIGENVALUE) & + call calculate_combined_keff() if (output_tallies) then if (master) call write_tallies() end if @@ -66,4 +68,107 @@ contains end subroutine finalize_run +!=============================================================================== +! CALCULATE_COMBINED_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). +!=============================================================================== + + subroutine calculate_combined_keff() + + integer :: i + integer :: n + real(8) :: k_abs, k_col, k_tra + real(8) :: s_abs, s_col, s_tra + real(8) :: s_abs_col + real(8) :: s_col_tra + real(8) :: s_tra_abs + real(8) :: s_ij, s_ik, s_jk + real(8) :: s_kk, s_jj + real(8) :: s_1i + real(8) :: k_i, k_j + real(8) :: f, g + real(8) :: S(3) + + ! Initialize variables + n = n_realizations + g = ZERO + S = ZERO + + ! Copy estimates of k-effective and its variance (not variance of the mean) + k_abs = global_tallies(K_ABSORPTION) % sum + k_col = global_tallies(K_COLLISION) % sum + k_tra = global_tallies(K_TRACKLENGTH) % sum + s_abs = global_tallies(K_ABSORPTION) % sum_sq**2 * n + s_col = global_tallies(K_COLLISION) % sum_sq**2 * n + s_tra = global_tallies(K_TRACKLENGTH) % sum_sq**2 * n + + ! Calculate covariances based on sums with Bessel's correction + s_abs_col = (k_abs_col - n * k_abs * k_col)/(n - 1) + s_col_tra = (k_col_tra - n * k_col * k_tra)/(n - 1) + s_tra_abs = (k_tra_abs - n * k_tra * k_abs)/(n - 1) + + do i = 1, 3 + if (i == 1) then + ! i = collision, j = absorption, k = tracklength + k_i = k_col + k_j = k_abs + s_1i = s_col + s_jj = s_abs + s_kk = s_tra + s_ij = s_abs_col + s_ik = s_col_tra + s_jk = s_tra_abs + elseif (i == 2) then + ! i = absortion, j = tracklength, k = collision + k_i = k_abs + k_j = k_tra + s_1i = s_abs_col + s_jj = s_tra + s_kk = s_col + s_ij = s_tra_abs + s_ik = s_abs_col + s_jk = s_col_tra + elseif (i == 3) then + ! i = tracklength, j = collision, k = absorption + k_i = k_tra + k_j = k_col + s_1i = s_col_tra + s_jj = s_col + s_kk = s_abs + s_ij = s_col_tra + s_ik = s_tra_abs + s_jk = s_abs_col + end if + + ! Calculate weighting + f = s_jj*(s_kk - s_ik) - s_kk*s_ij + s_jk*(s_ij + s_ik - s_jk) + + ! Add to S sums for variance of combined estimate + S(1) = S(1) + f*s_1i + S(2) = S(2) + (s_jj + s_kk - TWO*s_jk)*k_i*k_i + S(3) = S(3) + (s_kk + s_ij - s_jk - s_ik)*k_i*k_j + + ! Add to sum for combined k-effective + k_combined(1) = k_combined(1) + f*k_i + g = g + f + end do + + ! Complete calculations of S sums + S = (n - 1)*S + S(1) = (n - 1)**2 * S(1) + + ! Calculate combined estimate of k-effective + k_combined(1) = k_combined(1) / g + + ! Calculate standard deviation of combined estimate + g = (n - 1)**2 * g + k_combined(2) = sqrt(S(1)/(g*n*(n-3)) * (ONE + n*((S(2) - TWO*S(3))/g))) + + end subroutine calculate_combined_keff + end module finalize diff --git a/src/global.F90 b/src/global.F90 index ec537117fe..314187b45b 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -169,8 +169,12 @@ module global ! Temporary k-effective values real(8), allocatable :: k_batch(:) ! batch estimates of k - real(8) :: keff = ONE ! average k over active batches - real(8) :: keff_std ! standard deviation of average k + real(8) :: keff = ONE ! average k over active batches + real(8) :: keff_std ! standard deviation of average k + real(8) :: k_abs_col = ZERO ! sum over batches of k_absorption * k_collision + real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength + real(8) :: k_tra_abs = ZERO ! sum over batches of k_trackelngth * k_absorption + real(8) :: k_combined(2) ! combined best estimate of k-effective ! Shannon entropy logical :: entropy_on = .false. diff --git a/src/output.F90 b/src/output.F90 index ad0bc2cdaa..b20d923609 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1429,6 +1429,7 @@ contains % sum, global_tallies(K_TRACKLENGTH) % sum_sq write(ou,102) "k-effective (Absorption)", global_tallies(K_ABSORPTION) & % sum, global_tallies(K_ABSORPTION) % sum_sq + write(ou,102) "Combined k-effective", k_combined write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, & global_tallies(LEAKAGE) % sum_sq write(ou,*) diff --git a/src/tally.F90 b/src/tally.F90 index 528b0dff10..68b6e44bd3 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1735,6 +1735,9 @@ contains subroutine synchronize_tallies() integer :: i + real(8) :: k_abs ! Copy of batch absorption estimate of keff + real(8) :: k_col ! Copy of batch collision estimate of keff + real(8) :: k_tra ! Copy of batch tracklength estimate of keff #ifdef MPI ! Combine tally results onto master process @@ -1761,6 +1764,16 @@ contains ! accumulate_result k_batch(current_batch) = global_tallies(K_TRACKLENGTH) % value + + if (active_batches) then + ! Accumulate products of different estimators of k + k_abs = global_tallies(K_ABSORPTION) % value / total_weight + k_col = global_tallies(K_COLLISION) % value / total_weight + k_tra = global_tallies(K_TRACKLENGTH) % value / total_weight + k_abs_col = k_abs_col + k_abs * k_col + k_col_tra = k_col_tra + k_col * k_tra + k_tra_abs = k_tra_abs + k_tra * k_abs + end if end if ! Accumulate results for global tallies From ea1b5661093759c121d8cd1e3feeb69569d33302 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Feb 2013 15:31:55 -0500 Subject: [PATCH 2/8] Consistently order collision, absorption, tracklength. --- src/constants.F90 | 4 ++-- src/finalize.F90 | 30 +++++++++++++++--------------- src/global.F90 | 4 ++-- src/tally.F90 | 8 ++++---- 4 files changed, 23 insertions(+), 23 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index d1b913a21d..373a1f9c3d 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -299,8 +299,8 @@ module constants integer, parameter :: N_GLOBAL_TALLIES = 4 integer, parameter :: & K_COLLISION = 1, & - K_TRACKLENGTH = 2, & - K_ABSORPTION = 3, & + K_ABSORPTION = 2, & + K_TRACKLENGTH = 3, & LEAKAGE = 4 ! ============================================================================ diff --git a/src/finalize.F90 b/src/finalize.F90 index f65f6986c7..4bb850bcf6 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -82,11 +82,11 @@ contains integer :: i integer :: n - real(8) :: k_abs, k_col, k_tra - real(8) :: s_abs, s_col, s_tra - real(8) :: s_abs_col + real(8) :: k_col, k_abs, k_tra + real(8) :: s_col, s_abs, s_tra + real(8) :: s_col_abs real(8) :: s_col_tra - real(8) :: s_tra_abs + real(8) :: s_abs_tra real(8) :: s_ij, s_ik, s_jk real(8) :: s_kk, s_jj real(8) :: s_1i @@ -100,17 +100,17 @@ contains S = ZERO ! Copy estimates of k-effective and its variance (not variance of the mean) - k_abs = global_tallies(K_ABSORPTION) % sum k_col = global_tallies(K_COLLISION) % sum + k_abs = global_tallies(K_ABSORPTION) % sum k_tra = global_tallies(K_TRACKLENGTH) % sum - s_abs = global_tallies(K_ABSORPTION) % sum_sq**2 * n s_col = global_tallies(K_COLLISION) % sum_sq**2 * n + s_abs = global_tallies(K_ABSORPTION) % sum_sq**2 * n s_tra = global_tallies(K_TRACKLENGTH) % sum_sq**2 * n ! Calculate covariances based on sums with Bessel's correction - s_abs_col = (k_abs_col - n * k_abs * k_col)/(n - 1) + s_col_abs = (k_col_abs - n * k_col * k_abs)/(n - 1) s_col_tra = (k_col_tra - n * k_col * k_tra)/(n - 1) - s_tra_abs = (k_tra_abs - n * k_tra * k_abs)/(n - 1) + s_abs_tra = (k_abs_tra - n * k_abs * k_tra)/(n - 1) do i = 1, 3 if (i == 1) then @@ -120,18 +120,18 @@ contains s_1i = s_col s_jj = s_abs s_kk = s_tra - s_ij = s_abs_col + s_ij = s_col_abs s_ik = s_col_tra - s_jk = s_tra_abs + s_jk = s_abs_tra elseif (i == 2) then ! i = absortion, j = tracklength, k = collision k_i = k_abs k_j = k_tra - s_1i = s_abs_col + s_1i = s_col_abs s_jj = s_tra s_kk = s_col - s_ij = s_tra_abs - s_ik = s_abs_col + s_ij = s_abs_tra + s_ik = s_col_abs s_jk = s_col_tra elseif (i == 3) then ! i = tracklength, j = collision, k = absorption @@ -141,8 +141,8 @@ contains s_jj = s_col s_kk = s_abs s_ij = s_col_tra - s_ik = s_tra_abs - s_jk = s_abs_col + s_ik = s_abs_tra + s_jk = s_col_abs end if ! Calculate weighting diff --git a/src/global.F90 b/src/global.F90 index 314187b45b..31cd3535e3 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -171,9 +171,9 @@ module global real(8), allocatable :: k_batch(:) ! batch estimates of k real(8) :: keff = ONE ! average k over active batches real(8) :: keff_std ! standard deviation of average k - real(8) :: k_abs_col = ZERO ! sum over batches of k_absorption * k_collision + real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength - real(8) :: k_tra_abs = ZERO ! sum over batches of k_trackelngth * k_absorption + real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength real(8) :: k_combined(2) ! combined best estimate of k-effective ! Shannon entropy diff --git a/src/tally.F90 b/src/tally.F90 index 68b6e44bd3..ff46934f58 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1735,8 +1735,8 @@ contains subroutine synchronize_tallies() integer :: i - real(8) :: k_abs ! Copy of batch absorption estimate of keff real(8) :: k_col ! Copy of batch collision estimate of keff + real(8) :: k_abs ! Copy of batch absorption estimate of keff real(8) :: k_tra ! Copy of batch tracklength estimate of keff #ifdef MPI @@ -1767,12 +1767,12 @@ contains if (active_batches) then ! Accumulate products of different estimators of k - k_abs = global_tallies(K_ABSORPTION) % value / total_weight k_col = global_tallies(K_COLLISION) % value / total_weight + k_abs = global_tallies(K_ABSORPTION) % value / total_weight k_tra = global_tallies(K_TRACKLENGTH) % value / total_weight - k_abs_col = k_abs_col + k_abs * k_col + k_col_abs = k_col_abs + k_col * k_abs k_col_tra = k_col_tra + k_col * k_tra - k_tra_abs = k_tra_abs + k_tra * k_abs + k_abs_tra = k_abs_tra + k_abs * k_tra end if end if From 5b2e09391a2bedcfa39952324d24111f6ce8edc2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Feb 2013 17:19:13 -0500 Subject: [PATCH 3/8] Simplify permutations in calculate_combined_keff. --- src/finalize.F90 | 95 ++++++++++++++++++++---------------------------- 1 file changed, 40 insertions(+), 55 deletions(-) diff --git a/src/finalize.F90 b/src/finalize.F90 index 4bb850bcf6..1e9bedf422 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -80,19 +80,14 @@ contains subroutine calculate_combined_keff() - integer :: i - integer :: n - real(8) :: k_col, k_abs, k_tra - real(8) :: s_col, s_abs, s_tra - real(8) :: s_col_abs - real(8) :: s_col_tra - real(8) :: s_abs_tra - real(8) :: s_ij, s_ik, s_jk - real(8) :: s_kk, s_jj - real(8) :: s_1i - real(8) :: k_i, k_j - real(8) :: f, g - real(8) :: S(3) + integer :: l ! loop index + integer :: i, j, k ! indices referring to collision, absorption, or track + integer :: 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 ! Initialize variables n = n_realizations @@ -100,61 +95,51 @@ contains S = ZERO ! Copy estimates of k-effective and its variance (not variance of the mean) - k_col = global_tallies(K_COLLISION) % sum - k_abs = global_tallies(K_ABSORPTION) % sum - k_tra = global_tallies(K_TRACKLENGTH) % sum - s_col = global_tallies(K_COLLISION) % sum_sq**2 * n - s_abs = global_tallies(K_ABSORPTION) % sum_sq**2 * n - s_tra = global_tallies(K_TRACKLENGTH) % sum_sq**2 * n + kv(1) = global_tallies(K_COLLISION) % sum + kv(2) = global_tallies(K_ABSORPTION) % sum + kv(3) = global_tallies(K_TRACKLENGTH) % sum + cov(1,1) = global_tallies(K_COLLISION) % sum_sq**2 * n + cov(2,2) = global_tallies(K_ABSORPTION) % sum_sq**2 * n + cov(3,3) = global_tallies(K_TRACKLENGTH) % sum_sq**2 * n ! Calculate covariances based on sums with Bessel's correction - s_col_abs = (k_col_abs - n * k_col * k_abs)/(n - 1) - s_col_tra = (k_col_tra - n * k_col * k_tra)/(n - 1) - s_abs_tra = (k_abs_tra - n * k_abs * k_tra)/(n - 1) + cov(1,2) = (k_col_abs - n * kv(1) * kv(2))/(n - 1) + cov(1,3) = (k_col_tra - n * kv(1) * kv(3))/(n - 1) + cov(2,3) = (k_abs_tra - n * kv(2) * kv(3))/(n - 1) + cov(2,1) = cov(1,2) + cov(3,1) = cov(1,3) + cov(3,2) = cov(2,3) - do i = 1, 3 - if (i == 1) then + do l = 1, 3 + ! Permutations of estimates + if (l == 1) then ! i = collision, j = absorption, k = tracklength - k_i = k_col - k_j = k_abs - s_1i = s_col - s_jj = s_abs - s_kk = s_tra - s_ij = s_col_abs - s_ik = s_col_tra - s_jk = s_abs_tra - elseif (i == 2) then + i = 1 + j = 2 + k = 3 + elseif (l == 2) then ! i = absortion, j = tracklength, k = collision - k_i = k_abs - k_j = k_tra - s_1i = s_col_abs - s_jj = s_tra - s_kk = s_col - s_ij = s_abs_tra - s_ik = s_col_abs - s_jk = s_col_tra - elseif (i == 3) then + i = 2 + j = 3 + k = 1 + elseif (l == 3) then ! i = tracklength, j = collision, k = absorption - k_i = k_tra - k_j = k_col - s_1i = s_col_tra - s_jj = s_col - s_kk = s_abs - s_ij = s_col_tra - s_ik = s_abs_tra - s_jk = s_col_abs + i = 3 + j = 1 + k = 2 end if ! Calculate weighting - f = s_jj*(s_kk - s_ik) - s_kk*s_ij + s_jk*(s_ij + s_ik - s_jk) + 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*s_1i - S(2) = S(2) + (s_jj + s_kk - TWO*s_jk)*k_i*k_i - S(3) = S(3) + (s_kk + s_ij - s_jk - s_ik)*k_i*k_j + 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*k_i + k_combined(1) = k_combined(1) + f * kv(l) g = g + f end do From 4a55dae476a4ef8384651ffc142132af43826c2a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Feb 2013 17:22:23 -0500 Subject: [PATCH 4/8] Add confidence interval for combined keff estimator. --- src/output.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/output.F90 b/src/output.F90 index b20d923609..28b067c780 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1420,6 +1420,9 @@ contains ! Adjust sum_sq global_tallies(:) % sum_sq = t_value * global_tallies(:) % sum_sq + + ! Adjust combined estimator + k_combined(2) = t_value * k_combined(2) end if ! write global tallies From b4870354833f2df562e4d98b357cef9bc7aa4ad3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Feb 2013 18:21:41 -0500 Subject: [PATCH 5/8] Moved calculation of combined keff before tally_statistics. --- src/finalize.F90 | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/finalize.F90 b/src/finalize.F90 index 1e9bedf422..30609b971c 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -31,9 +31,9 @@ contains if (run_mode /= MODE_PLOTTING) then ! Calculate statistics for tallies and write to tallies.out - if (master) call tally_statistics() if (master .and. run_mode == MODE_EIGENVALUE) & call calculate_combined_keff() + if (master) call tally_statistics() if (output_tallies) then if (master) call write_tallies() end if @@ -93,14 +93,18 @@ contains n = n_realizations g = ZERO S = ZERO + k_combined = ZERO ! Copy estimates of k-effective and its variance (not variance of the mean) - kv(1) = global_tallies(K_COLLISION) % sum - kv(2) = global_tallies(K_ABSORPTION) % sum - kv(3) = global_tallies(K_TRACKLENGTH) % sum - cov(1,1) = global_tallies(K_COLLISION) % sum_sq**2 * n - cov(2,2) = global_tallies(K_ABSORPTION) % sum_sq**2 * n - cov(3,3) = global_tallies(K_TRACKLENGTH) % sum_sq**2 * n + kv(1) = global_tallies(K_COLLISION) % sum / n + kv(2) = global_tallies(K_ABSORPTION) % sum / n + kv(3) = global_tallies(K_TRACKLENGTH) % sum / n + cov(1,1) = (global_tallies(K_COLLISION) % sum_sq - & + n * kv(1) * kv(1)) / (n - 1) + cov(2,2) = (global_tallies(K_ABSORPTION) % sum_sq - & + n * kv(2) * kv(2)) / (n - 1) + cov(3,3) = (global_tallies(K_TRACKLENGTH) % sum_sq - & + n * kv(3) * kv(3)) / (n - 1) ! Calculate covariances based on sums with Bessel's correction cov(1,2) = (k_col_abs - n * kv(1) * kv(2))/(n - 1) From cf2af635dff59e4f6404163d87000cc43f18b7ff Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Feb 2013 18:31:25 -0500 Subject: [PATCH 6/8] Move calculate_combined_keff to eigenvalue module. --- src/eigenvalue.F90 | 101 +++++++++++++++++++++++++++++++++++++++++++++ src/finalize.F90 | 94 ----------------------------------------- 2 files changed, 101 insertions(+), 94 deletions(-) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index bc41a9128e..dd509b8396 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -184,6 +184,9 @@ contains ! Write out state point if it's been specified for this batch do i = 1, n_state_points if (current_batch == statepoint_batch(i)) then + ! Calculate combined estimate of k-effective + if (master) call calculate_combined_keff() + ! Create state point file #ifdef HDF5 call hdf5_write_state_point() @@ -194,6 +197,12 @@ contains end if end do + if (master .and. current_batch == n_batches) then + ! Make sure combined estimate of k-effective is calculated at the last + ! batch in case no state point is written + call calculate_combined_keff() + end if + end subroutine finalize_batch !=============================================================================== @@ -636,6 +645,98 @@ contains end subroutine calculate_keff +!=============================================================================== +! CALCULATE_COMBINED_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). +!=============================================================================== + + subroutine calculate_combined_keff() + + integer :: l ! loop index + integer :: i, j, k ! indices referring to collision, absorption, or track + integer :: 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 + + ! Initialize variables + n = n_realizations + g = ZERO + S = ZERO + k_combined = ZERO + + ! Copy estimates of k-effective and its variance (not variance of the mean) + kv(1) = global_tallies(K_COLLISION) % sum / n + kv(2) = global_tallies(K_ABSORPTION) % sum / n + kv(3) = global_tallies(K_TRACKLENGTH) % sum / n + cov(1,1) = (global_tallies(K_COLLISION) % sum_sq - & + n * kv(1) * kv(1)) / (n - 1) + cov(2,2) = (global_tallies(K_ABSORPTION) % sum_sq - & + n * kv(2) * kv(2)) / (n - 1) + cov(3,3) = (global_tallies(K_TRACKLENGTH) % sum_sq - & + n * kv(3) * kv(3)) / (n - 1) + + ! Calculate covariances based on sums with Bessel's correction + cov(1,2) = (k_col_abs - n * kv(1) * kv(2))/(n - 1) + cov(1,3) = (k_col_tra - n * kv(1) * kv(3))/(n - 1) + cov(2,3) = (k_abs_tra - n * kv(2) * kv(3))/(n - 1) + cov(2,1) = cov(1,2) + cov(3,1) = cov(1,3) + cov(3,2) = cov(2,3) + + 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 - 1)*S + S(1) = (n - 1)**2 * S(1) + + ! Calculate combined estimate of k-effective + k_combined(1) = k_combined(1) / g + + ! Calculate standard deviation of combined estimate + g = (n - 1)**2 * g + k_combined(2) = sqrt(S(1)/(g*n*(n-3)) * (ONE + n*((S(2) - TWO*S(3))/g))) + + end subroutine calculate_combined_keff + !=============================================================================== ! COUNT_SOURCE_FOR_UFS determines the source fraction in each UFS mesh cell and ! reweights the source bank so that the sum of the weights is equal to diff --git a/src/finalize.F90 b/src/finalize.F90 index 30609b971c..1bfdfe71f4 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -31,8 +31,6 @@ contains if (run_mode /= MODE_PLOTTING) then ! Calculate statistics for tallies and write to tallies.out - if (master .and. run_mode == MODE_EIGENVALUE) & - call calculate_combined_keff() if (master) call tally_statistics() if (output_tallies) then if (master) call write_tallies() @@ -68,96 +66,4 @@ contains end subroutine finalize_run -!=============================================================================== -! CALCULATE_COMBINED_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). -!=============================================================================== - - subroutine calculate_combined_keff() - - integer :: l ! loop index - integer :: i, j, k ! indices referring to collision, absorption, or track - integer :: 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 - - ! Initialize variables - n = n_realizations - g = ZERO - S = ZERO - k_combined = ZERO - - ! Copy estimates of k-effective and its variance (not variance of the mean) - kv(1) = global_tallies(K_COLLISION) % sum / n - kv(2) = global_tallies(K_ABSORPTION) % sum / n - kv(3) = global_tallies(K_TRACKLENGTH) % sum / n - cov(1,1) = (global_tallies(K_COLLISION) % sum_sq - & - n * kv(1) * kv(1)) / (n - 1) - cov(2,2) = (global_tallies(K_ABSORPTION) % sum_sq - & - n * kv(2) * kv(2)) / (n - 1) - cov(3,3) = (global_tallies(K_TRACKLENGTH) % sum_sq - & - n * kv(3) * kv(3)) / (n - 1) - - ! Calculate covariances based on sums with Bessel's correction - cov(1,2) = (k_col_abs - n * kv(1) * kv(2))/(n - 1) - cov(1,3) = (k_col_tra - n * kv(1) * kv(3))/(n - 1) - cov(2,3) = (k_abs_tra - n * kv(2) * kv(3))/(n - 1) - cov(2,1) = cov(1,2) - cov(3,1) = cov(1,3) - cov(3,2) = cov(2,3) - - 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 - 1)*S - S(1) = (n - 1)**2 * S(1) - - ! Calculate combined estimate of k-effective - k_combined(1) = k_combined(1) / g - - ! Calculate standard deviation of combined estimate - g = (n - 1)**2 * g - k_combined(2) = sqrt(S(1)/(g*n*(n-3)) * (ONE + n*((S(2) - TWO*S(3))/g))) - - end subroutine calculate_combined_keff - end module finalize From 839276ef07bfe3d3d781e2b8dd113066c28e5536 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Feb 2013 18:35:36 -0500 Subject: [PATCH 7/8] Make sure n_realizations is at least 4 for combined keff calculation. --- src/eigenvalue.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index dd509b8396..d97b5c3511 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -666,6 +666,10 @@ contains real(8) :: g ! sum of weighting factors real(8) :: S(3) ! sums used for variance calculation + ! 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 + ! Initialize variables n = n_realizations g = ZERO From 5d5425b0500e33d399b15702e7ec1611980d49b5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Feb 2013 19:08:28 -0500 Subject: [PATCH 8/8] Added combined estimate of k-effective to state point. --- docs/source/devguide/statepoint.rst | 268 ++++++++++++++++++++++++++-- src/constants.F90 | 2 +- src/hdf5_interface.F90 | 9 + src/state_point.F90 | 28 +++ src/utils/statepoint.py | 4 + 5 files changed, 300 insertions(+), 11 deletions(-) diff --git a/docs/source/devguide/statepoint.rst b/docs/source/devguide/statepoint.rst index 8b48e908d5..7fbfe14a15 100644 --- a/docs/source/devguide/statepoint.rst +++ b/docs/source/devguide/statepoint.rst @@ -4,6 +4,254 @@ State Point Binary File Specifications ====================================== +---------- +Revision 8 +---------- + +**integer(4) REVISION_STATEPOINT** + + Revision of the binary state point file. Any time a change is made in the + format of the state-point file, this integer is incremented. + +**integer(4) VERSION_MAJOR** + + Major version number for OpenMC + +**integer(4) VERSION_MINOR** + + Minor version number for OpenMC + +**integer(4) VERSION_RELEASE** + + Release version number for OpenMC + +**character(19) time_stamp** + + Date and time the state point was written. + +**character(255) path** + + Absolute path to directory containing input files. + +**integer(8) seed** + + Pseudo-random number generator seed. + +**integer(4) run_mode** + + run mode used. The modes are described in constants.F90. + +**integer(8) n_particles** + + Number of particles used per generation. + +**integer(4) n_batches** + + Total number of batches (active + inactive). + +**integer(4) current_batch** + + The number of batches already simulated. + +if (run_mode == MODE_EIGENVALUE) + + **integer(4) n_inactive** + + Number of inactive batches + + **integer(4) gen_per_batch** + + Number of generations per batch for criticality calculations + + *do i = 1, current_batch* + + **real(8) k_batch(i)** + + k-effective for the i-th batch + + *do i = 1, current_batch* + + **real(8) entropy(i)** + + Shannon entropy for the i-th batch + + **real(8) k_col_abs** + + Sum of product of collision/absorption estimates of k-effective + + **real(8) k_col_tra** + + Sum of product of collision/track-length estimates of k-effective + + **real(8) k_abs_tra** + + Sum of product of absorption/track-length estimates of k-effective + + **real(8) k_combined(2)** + + Mean and standard deviation of a combined estimate of k-effective + +**integer(4) n_meshes** + + Number of meshes in tallies.xml file + +*do i = 1, n_meshes* + + **integer(4) meshes(i) % id** + + Unique ID of mesh. + + **integer(4) meshes(i) % type** + + Type of mesh. + + **integer(4) meshes(i) % n_dimension** + + Number of dimensions for mesh (2 or 3). + + **integer(4) meshes(i) % dimension(:)** + + Number of mesh cells in each dimension. + + **real(8) meshes(i) % lower_left(:)** + + Coordinates of lower-left corner of mesh. + + **real(8) meshes(i) % upper_right(:)** + + Coordinates of upper-right corner of mesh. + + **real(8) meshes(i) % width(:)** + + Width of each mesh cell in each dimension. + +**integer(4) n_tallies** + +*do i = 1, n_tallies* + + **integer(4) tallies(i) % id** + + Unique ID of tally. + + **integer(4) tallies(i) % n_realizations** + + Number of realizations for the i-th tally. + + **integer(4) size(tallies(i) % scores, 1)** + + Total number of score bins for the i-th tally + + **integer(4) size(tallies(i) % scores, 2)** + + Total number of filter bins for the i-th tally + + **integer(4) tallies(i) % n_filters** + + *do j = 1, tallies(i) % n_filters* + + **integer(4) tallies(i) % filter(j) % type** + + Type of tally filter. + + **integer(4) tallies(i) % filter(j) % n_bins** + + Number of bins for filter. + + **integer(4)/real(8) tallies(i) % filter(j) % bins(:)** + + Value for each filter bin of this type. + + **integer(4) tallies(i) % n_nuclide_bins** + + Number of nuclide bins. If none are specified, this is just one. + + *do j = 1, tallies(i) % n_nuclide_bins* + + **integer(4) tallies(i) % nuclide_bins(j)** + + Values of specified nuclide bins + + **integer(4) tallies(i) % n_score_bins** + + Number of scoring bins. + + *do j = 1, tallies(i) % n_score_bins* + + **integer(4) tallies(i) % score_bins(j)** + + Values of specified scoring bins (e.g. SCORE_FLUX). + + *do j = 1, tallies(i) % n_score_bins* + + **integer(4) tallies(i) % scatt_order(j)** + + Scattering Order specified scoring bins. + + **integer(4) tallies(i) % n_score_bins** + + Number of scoring bins without accounting for those added by + the scatter-pn command. + +**integer(4) n_realizations** + + Number of realizations for global tallies. + +**integer(4) N_GLOBAL_TALLIES** + + Number of global tally scores + +*do i = 1, N_GLOBAL_TALLIES* + + **real(8) global_tallies(i) % sum** + + Accumulated sum for the i-th global tally + + **real(8) global_tallies(i) % sum_sq** + + Accumulated sum of squares for the i-th global tally + +**integer(4) tallies_on** + + Flag indicated if tallies are present in the file. + +if (tallies_on > 0) + + *do i = 1, n_tallies* + + *do k = 1, size(tallies(i) % scores, 2)* + + *do j = 1, size(tallies(i) % scores, 1)* + + **real(8) tallies(i) % scores(j,k) % sum** + + Accumulated sum for the j-th score and k-th filter of the + i-th tally + + **real(8) tallies(i) % scores(j,k) % sum_sq** + + Accumulated sum of squares for the j-th score and k-th + filter of the i-th tally + +if (run_mode == MODE_EIGENVALUE) + + *do i = 1, n_particles* + + **real(8) source_bank(i) % wgt** + + Weight of the i-th source particle + + **real(8) source_bank(i) % xyz(1:3)** + + Coordinates of the i-th source particle. + + **real(8) source_bank(i) % uvw(1:3)** + + Direction of the i-th source particle + + **real(8) source_bank(i) % E** + + Energy of the i-th source particle. + ---------- Revision 7 ---------- @@ -53,7 +301,7 @@ Revision 7 The number of batches already simulated. -if (run_mode == MODE_CRITICALITY) +if (run_mode == MODE_EIGENVALUE) **integer(4) n_inactive** @@ -216,7 +464,7 @@ if (tallies_on > 0) Accumulated sum of squares for the j-th score and k-th filter of the i-th tally -if (run_mode == MODE_CRITICALITY) +if (run_mode == MODE_EIGENVALUE) *do i = 1, n_particles* @@ -285,7 +533,7 @@ Revision 6 The number of batches already simulated. -if (run_mode == MODE_CRITICALITY) +if (run_mode == MODE_EIGENVALUE) **integer(4) n_inactive** @@ -437,7 +685,7 @@ if (tallies_on > 0) Accumulated sum of squares for the j-th score and k-th filter of the i-th tally -if (run_mode == MODE_CRITICALITY) +if (run_mode == MODE_EIGENVALUE) *do i = 1, n_particles* @@ -502,7 +750,7 @@ Revision 5 The number of batches already simulated. -if (run_mode == MODE_CRITICALITY) +if (run_mode == MODE_EIGENVALUE) **integer(4) n_inactive** @@ -646,7 +894,7 @@ if (tallies_on > 0) Accumulated sum of squares for the j-th score and k-th filter of the i-th tally -if (run_mode == MODE_CRITICALITY) +if (run_mode == MODE_EIGENVALUE) *do i = 1, n_particles* @@ -711,7 +959,7 @@ Revision 4 The number of batches already simulated. -if (run_mode == MODE_CRITICALITY) +if (run_mode == MODE_EIGENVALUE) **integer(4) n_inactive** @@ -851,7 +1099,7 @@ if (tallies_on > 0) Accumulated sum of squares for the j-th score and k-th filter of the i-th tally -if (run_mode == MODE_CRITICALITY) +if (run_mode == MODE_EIGENVALUE) *do i = 1, n_particles* @@ -916,7 +1164,7 @@ Revision 3 The number of batches already simulated. -if (run_mode == MODE_CRITICALITY) +if (run_mode == MODE_EIGENVALUE) **integer(4) n_inactive** @@ -1052,7 +1300,7 @@ if (tallies_on > 0) Accumulated sum of squares for the j-th score and k-th filter of the i-th tally -if (run_mode == MODE_CRITICALITY) +if (run_mode == MODE_EIGENVALUE) *do i = 1, n_particles* diff --git a/src/constants.F90 b/src/constants.F90 index 373a1f9c3d..95c10f5281 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -12,7 +12,7 @@ module constants ! Revision numbers for binary files integer, parameter :: REVISION_SOURCE = 1 - integer, parameter :: REVISION_STATEPOINT = 7 + integer, parameter :: REVISION_STATEPOINT = 8 ! ============================================================================ ! ADJUSTABLE PARAMETERS diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index bd093434d8..5ca2b0c682 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -887,6 +887,12 @@ contains dims(1) = current_batch*gen_per_batch call h5ltmake_dataset_double_f(hdf5_state_point, "entropy", 1, & dims, entropy, hdf5_err) + call hdf5_write_double(hdf5_state_point, "k_col_abs", k_col_abs) + call hdf5_write_double(hdf5_state_point, "k_col_tra", k_col_tra) + call hdf5_write_double(hdf5_state_point, "k_abs_tra", k_abs_tra) + dims(1) = 2 + call h5ltmake_dataset_double_f(hdf5_state_point, "k_combined", 1, & + dims, k_combined, hdf5_err) end if ! Create group for tallies @@ -1189,6 +1195,9 @@ contains k_batch(1:restart_batch), dims, hdf5_err) call h5ltread_dataset_double_f(hdf5_state_point, "entropy", & entropy(1:restart_batch*gen_per_batch), dims, hdf5_err) + call hdf5_read_double(hdf5_state_point, "k_col_abs", k_col_abs) + call hdf5_read_double(hdf5_state_point, "k_col_tra", k_col_tra) + call hdf5_read_double(hdf5_state_point, "k_abs_tra", k_abs_tra) end if ! Read number of realizations for global tallies diff --git a/src/state_point.F90 b/src/state_point.F90 index 7b3e4a5a94..9a92e64f09 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -187,6 +187,10 @@ contains write(UNIT_STATE) n_inactive, gen_per_batch write(UNIT_STATE) k_batch(1:current_batch) write(UNIT_STATE) entropy(1:current_batch*gen_per_batch) + write(UNIT_STATE) k_col_abs + write(UNIT_STATE) k_col_tra + write(UNIT_STATE) k_abs_tra + write(UNIT_STATE) k_combined end if ! Write number of meshes @@ -374,6 +378,14 @@ contains MPI_STATUS_IGNORE, mpi_err) call MPI_FILE_WRITE(fh, entropy, current_batch*gen_per_batch, MPI_REAL8, & MPI_STATUS_IGNORE, mpi_err) + call MPI_FILE_WRITE(fh, k_col_abs, 1, MPI_REAL8, & + MPI_STATUS_IGNORE, mpi_err) + call MPI_FILE_WRITE(fh, k_col_tra, 1, MPI_REAL8, & + MPI_STATUS_IGNORE, mpi_err) + call MPI_FILE_WRITE(fh, k_abs_tra, 1, MPI_REAL8, & + MPI_STATUS_IGNORE, mpi_err) + call MPI_FILE_WRITE(fh, k_combined, 2, MPI_REAL8, & + MPI_STATUS_IGNORE, mpi_err) end if ! Write number of meshes @@ -685,6 +697,16 @@ contains MPI_STATUS_IGNORE, mpi_err) call MPI_FILE_READ_ALL(fh, entropy, restart_batch*gen_per_batch, & MPI_REAL8, MPI_STATUS_IGNORE, mpi_err) + call MPI_FILE_READ_ALL(fh, k_col_abs, 1, MPI_REAL8, & + MPI_STATUS_IGNORE, mpi_err) + call MPI_FILE_READ_ALL(fh, k_col_tra, 1, MPI_REAL8, & + MPI_STATUS_IGNORE, mpi_err) + call MPI_FILE_READ_ALL(fh, k_abs_tra, 1, MPI_REAL8, & + MPI_STATUS_IGNORE, mpi_err) + allocate(real_array(2)) + call MPI_FILE_READ_ALL(fh, real_array, 2, MPI_REAL8, & + MPI_STATUS_IGNORE, mpi_err) + deallocate(real_array) end if if (master) then @@ -892,6 +914,12 @@ contains read(UNIT_STATE) n_inactive, gen_per_batch read(UNIT_STATE) k_batch(1:restart_batch) read(UNIT_STATE) entropy(1:restart_batch*gen_per_batch) + read(UNIT_STATE) k_col_abs + read(UNIT_STATE) k_col_tra + read(UNIT_STATE) k_abs_tra + allocate(real_array(2)) + read(UNIT_STATE) real_array + deallocate(real_array) end if if (master) then diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 7e6d11d30a..1fb814c5fc 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -195,6 +195,10 @@ class StatePoint(BinaryFile): self.n_inactive, self.gen_per_batch = self._get_int(2) self.k_batch = self._get_double(self.current_batch) self.entropy = self._get_double(self.current_batch) + self.k_col_abs = self._get_double()[0] + self.k_col_tra = self._get_double()[0] + self.k_abs_tra = self._get_double()[0] + self.k_combined = self._get_double(2) # Read number of meshes n_meshes = self._get_int()[0]