mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Confidence intervals for keff. Moved quantile functions to new math module.
This commit is contained in:
Paul Romano
parent
5d983d4be8
commit
c1e2a1cfcd
7 changed files with 168 additions and 107 deletions
|
|
@ -145,6 +145,7 @@ input_xml.o: xml-fortran/templates/tallies_t.o
|
|||
|
||||
intercycle.o: error.o
|
||||
intercycle.o: global.o
|
||||
intercycle.o: math.o
|
||||
intercycle.o: mesh.o
|
||||
intercycle.o: mesh_header.o
|
||||
intercycle.o: output.o
|
||||
|
|
@ -170,6 +171,8 @@ main.o: global.o
|
|||
main.o: initialize.o
|
||||
main.o: plot.o
|
||||
|
||||
math.o: constants.o
|
||||
|
||||
mesh.o: constants.o
|
||||
mesh.o: global.o
|
||||
mesh.o: mesh_header.o
|
||||
|
|
@ -182,6 +185,7 @@ output.o: datatypes.o
|
|||
output.o: endf.o
|
||||
output.o: geometry_header.o
|
||||
output.o: global.o
|
||||
output.o: math.o
|
||||
output.o: mesh_header.o
|
||||
output.o: particle_header.o
|
||||
output.o: plot_header.o
|
||||
|
|
@ -239,6 +243,7 @@ source.o: string.o
|
|||
|
||||
state_point.o: error.o
|
||||
state_point.o: global.o
|
||||
state_point.o: math.o
|
||||
state_point.o: output.o
|
||||
state_point.o: string.o
|
||||
state_point.o: tally_header.o
|
||||
|
|
@ -250,6 +255,7 @@ string.o: global.o
|
|||
tally.o: constants.o
|
||||
tally.o: error.o
|
||||
tally.o: global.o
|
||||
tally.o: math.o
|
||||
tally.o: mesh.o
|
||||
tally.o: mesh_header.o
|
||||
tally.o: output.o
|
||||
|
|
|
|||
|
|
@ -24,6 +24,7 @@ interpolation.o \
|
|||
input_xml.o \
|
||||
main.o \
|
||||
material_header.o \
|
||||
math.o \
|
||||
mesh_header.o \
|
||||
mesh.o \
|
||||
output.o \
|
||||
|
|
|
|||
|
|
@ -4,6 +4,7 @@ module intercycle
|
|||
|
||||
use error, only: fatal_error, warning
|
||||
use global
|
||||
use math, only: t_percentile
|
||||
use mesh, only: count_bank_sites
|
||||
use mesh_header, only: StructuredMesh
|
||||
use output, only: write_message, print_batch_keff
|
||||
|
|
@ -361,6 +362,8 @@ contains
|
|||
subroutine calculate_keff()
|
||||
|
||||
real(8) :: temp(2) ! used to reduce sum and sum_sq
|
||||
real(8) :: alpha ! significance level for CI
|
||||
real(8) :: t_value ! t-value for confidence intervals
|
||||
|
||||
message = "Calculate batch keff..."
|
||||
call write_message(8)
|
||||
|
|
@ -401,8 +404,16 @@ contains
|
|||
keff = global_tallies(K_ANALOG) % sum / n_realizations
|
||||
|
||||
if (n_realizations > 1) then
|
||||
if (confidence_intervals) then
|
||||
! Calculate t-value for confidence intervals
|
||||
alpha = ONE - CONFIDENCE_LEVEL
|
||||
t_value = t_percentile(ONE - alpha/TWO, n_realizations)
|
||||
else
|
||||
t_value = ONE
|
||||
end if
|
||||
|
||||
! Standard deviation of the sample mean of k
|
||||
keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq / &
|
||||
keff_std = t_value * sqrt((global_tallies(K_ANALOG) % sum_sq / &
|
||||
n_realizations - keff * keff) / (n_realizations - 1))
|
||||
end if
|
||||
else
|
||||
|
|
@ -427,8 +438,16 @@ contains
|
|||
keff = temp(1) / n_realizations
|
||||
|
||||
if (n_realizations > 1) then
|
||||
if (confidence_intervals) then
|
||||
! Calculate t-value for confidence intervals
|
||||
alpha = ONE - CONFIDENCE_LEVEL
|
||||
t_value = t_percentile(ONE - alpha/TWO, n_realizations)
|
||||
else
|
||||
t_value = ONE
|
||||
end if
|
||||
|
||||
! Standard deviation of the sample mean of k
|
||||
keff_std = sqrt((temp(2)/n_realizations - keff*keff) / &
|
||||
keff_std = t_value * sqrt((temp(2)/n_realizations - keff*keff) / &
|
||||
(n_realizations - 1))
|
||||
end if
|
||||
end if
|
||||
|
|
|
|||
112
src/math.F90
Normal file
112
src/math.F90
Normal file
|
|
@ -0,0 +1,112 @@
|
|||
module math
|
||||
|
||||
use constants, only: PI, ONE, TWO
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! NORMAL_PERCENTILE calculates the percentile of the standard normal
|
||||
! distribution with a specified probability level
|
||||
!===============================================================================
|
||||
|
||||
function normal_percentile(p) result(z)
|
||||
|
||||
real(8), intent(in) :: p ! probability level
|
||||
real(8) :: z ! corresponding z-value
|
||||
|
||||
real(8) :: q
|
||||
real(8) :: r
|
||||
real(8), parameter :: p_low = 0.02425_8
|
||||
real(8), parameter :: a(6) = (/ &
|
||||
-3.969683028665376e1_8, 2.209460984245205e2_8, -2.759285104469687e2_8, &
|
||||
1.383577518672690e2_8, -3.066479806614716e1_8, 2.506628277459239e0_8 /)
|
||||
real(8), parameter :: b(5) = (/ &
|
||||
-5.447609879822406e1_8, 1.615858368580409e2_8, -1.556989798598866e2_8, &
|
||||
6.680131188771972e1_8, -1.328068155288572e1_8 /)
|
||||
real(8), parameter :: c(6) = (/ &
|
||||
-7.784894002430293e-3_8, -3.223964580411365e-1_8, -2.400758277161838_8, &
|
||||
-2.549732539343734_8, 4.374664141464968_8, 2.938163982698783_8 /)
|
||||
real(8), parameter :: d(4) = (/ &
|
||||
7.784695709041462e-3_8, 3.224671290700398e-1_8, &
|
||||
2.445134137142996_8, 3.754408661907416_8 /)
|
||||
|
||||
if (p < p_low) then
|
||||
! Rational approximation for lower region.
|
||||
|
||||
q = sqrt(-TWO*log(p))
|
||||
z = (((((c(1)*q + c(2))*q + c(3))*q + c(4))*q + c(5))*q + c(6)) / &
|
||||
((((d(1)*q + d(2))*q + d(3))*q + d(4))*q + 1.)
|
||||
|
||||
elseif (p <= 1. - p_low) then
|
||||
! Rational approximation for central region
|
||||
|
||||
q = p - 0.5
|
||||
r = q*q
|
||||
z = (((((a(1)*r + a(2))*r + a(3))*r + a(4))*r + a(5))*r + a(6))*q / &
|
||||
(((((b(1)*r + b(2))*r + b(3))*r + b(4))*r + b(5))*r + 1.)
|
||||
|
||||
else
|
||||
! Rational approximation for upper region
|
||||
|
||||
q = sqrt(-2*log(1. - p))
|
||||
z = -(((((c(1)*q + c(2))*q + c(3))*q + c(4))*q + c(5))*q + c(6)) / &
|
||||
((((d(1)*q + d(2))*q + d(3))*q + d(4))*q + 1.)
|
||||
endif
|
||||
|
||||
! Refinement based on Newton's method
|
||||
#ifndef NO_F2008
|
||||
z = z - (0.5 * erfc(-z/sqrt(TWO)) - p) * sqrt(TWO*PI) * exp(0.5*z*z)
|
||||
#endif
|
||||
|
||||
end function normal_percentile
|
||||
|
||||
!===============================================================================
|
||||
! T_PERCENTILE calculates the percentile of the Student's t distribution with a
|
||||
! specified probability level and number of degrees of freedom
|
||||
!===============================================================================
|
||||
|
||||
function t_percentile(p, df) result(t)
|
||||
|
||||
real(8), intent(in) :: p ! probability level
|
||||
integer, intent(in) :: df ! degrees of freedom
|
||||
real(8) :: t ! corresponding t-value
|
||||
|
||||
real(8) :: n ! degrees of freedom as a real(8)
|
||||
real(8) :: k ! n - 2
|
||||
real(8) :: z ! percentile of normal distribution
|
||||
real(8) :: z2 ! z * z
|
||||
|
||||
if (df == 1) then
|
||||
! For one degree of freedom, the t-distribution becomes a Cauchy
|
||||
! distribution whose cdf we can invert directly
|
||||
|
||||
t = tan(PI*(p - 0.5))
|
||||
|
||||
elseif (df == 2) then
|
||||
! For two degrees of freedom, the cdf is given by 1/2 + x/(2*sqrt(x^2 +
|
||||
! 2)). This can be directly inverted to yield the solution below
|
||||
|
||||
t = TWO*sqrt(TWO)*(p - 0.5)/sqrt(ONE - 4.*(p - 0.5)**2)
|
||||
|
||||
else
|
||||
|
||||
! This approximation is from E. Olusegun George and Meenakshi Sivaram, "A
|
||||
! modification of the Fisher-Cornish approximation for the student t
|
||||
! percentiles," Communication in Statistics - Simulation and Computation,
|
||||
! 16 (4), pp. 1123-1132 (1987).
|
||||
|
||||
n = real(df,8)
|
||||
k = 1./(n - 2.)
|
||||
z = normal_percentile(p)
|
||||
z2 = z * z
|
||||
t = sqrt(n*k) * (z + (z2 - 3.)*z*k/4. + ((5.*z2 - 56.)*z2 + &
|
||||
75.)*z*k*k/96. + (((z2 - 27.)*3.*z2 + 417.)*z2 - 315.) &
|
||||
*z*k*k*k/384.)
|
||||
|
||||
end if
|
||||
|
||||
end function t_percentile
|
||||
|
||||
end module math
|
||||
|
|
@ -8,6 +8,7 @@ module output
|
|||
use endf, only: reaction_name
|
||||
use geometry_header, only: Cell, Universe, Surface, BASE_UNIVERSE
|
||||
use global
|
||||
use math, only: t_percentile
|
||||
use mesh_header, only: StructuredMesh
|
||||
use particle_header, only: LocalCoord
|
||||
use plot_header
|
||||
|
|
@ -1194,6 +1195,8 @@ contains
|
|||
|
||||
integer(8) :: total_particles ! total # of particles simulated
|
||||
real(8) :: speed ! # of neutrons/second
|
||||
real(8) :: alpha ! significance level for CI
|
||||
real(8) :: t_value ! t-value for confidence intervals
|
||||
character(15) :: string
|
||||
|
||||
! display header block
|
||||
|
|
@ -1231,6 +1234,15 @@ contains
|
|||
! display header block for results
|
||||
call header("Results")
|
||||
|
||||
if (confidence_intervals) then
|
||||
! Calculate t-value for confidence intervals
|
||||
alpha = ONE - CONFIDENCE_LEVEL
|
||||
t_value = t_percentile(ONE - alpha/TWO, n_realizations)
|
||||
|
||||
! Adjust sum_sq
|
||||
global_tallies(:) % sum_sq = t_value * global_tallies(:) % sum_sq
|
||||
end if
|
||||
|
||||
! write global tallies
|
||||
write(ou,102) "k-effective (Analog)", global_tallies(K_ANALOG) % sum, &
|
||||
global_tallies(K_ANALOG) % sum_sq
|
||||
|
|
|
|||
|
|
@ -2,6 +2,7 @@ module state_point
|
|||
|
||||
use error, only: warning, fatal_error
|
||||
use global
|
||||
use math, only: t_percentile
|
||||
use output, only: write_message, print_batch_keff
|
||||
use string, only: to_str
|
||||
use tally_header, only: TallyObject
|
||||
|
|
@ -318,7 +319,9 @@ contains
|
|||
|
||||
subroutine replay_batch_history
|
||||
|
||||
real(8), save :: temp(2) = ZERO
|
||||
real(8), save :: temp(2) = ZERO ! temporary values for keff
|
||||
real(8) :: alpha ! significance level for CI
|
||||
real(8) :: t_value ! t-value for confidence intervals
|
||||
|
||||
! Write message at beginning
|
||||
if (current_batch == 1) then
|
||||
|
|
@ -337,8 +340,18 @@ contains
|
|||
temp(1) = temp(1) + k_batch(current_batch)
|
||||
temp(2) = temp(2) + k_batch(current_batch)*k_batch(current_batch)
|
||||
|
||||
! calculate mean keff
|
||||
keff = temp(1) / n_realizations
|
||||
keff_std = sqrt((temp(2)/n_realizations - keff*keff) &
|
||||
|
||||
if (confidence_intervals) then
|
||||
! Calculate t-value for confidence intervals
|
||||
alpha = ONE - CONFIDENCE_LEVEL
|
||||
t_value = t_percentile(ONE - alpha/TWO, n_realizations)
|
||||
else
|
||||
t_value = ONE
|
||||
end if
|
||||
|
||||
keff_std = t_value * sqrt((temp(2)/n_realizations - keff*keff) &
|
||||
/ (n_realizations - 1))
|
||||
else
|
||||
keff = k_batch(current_batch)
|
||||
|
|
|
|||
104
src/tally.F90
104
src/tally.F90
|
|
@ -3,6 +3,7 @@ module tally
|
|||
use constants
|
||||
use error, only: fatal_error
|
||||
use global
|
||||
use math, only: t_percentile
|
||||
use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices, &
|
||||
mesh_indices_to_bin, mesh_intersects
|
||||
use mesh_header, only: StructuredMesh
|
||||
|
|
@ -2320,107 +2321,4 @@ contains
|
|||
|
||||
end subroutine reset_score
|
||||
|
||||
!===============================================================================
|
||||
! NORMAL_PERCENTILE calculates the percentile of the standard normal
|
||||
! distribution with a specified probability level
|
||||
!===============================================================================
|
||||
|
||||
function normal_percentile(p) result(z)
|
||||
|
||||
real(8), intent(in) :: p ! probability level
|
||||
real(8) :: z ! corresponding z-value
|
||||
|
||||
real(8) :: q
|
||||
real(8) :: r
|
||||
real(8), parameter :: p_low = 0.02425_8
|
||||
real(8), parameter :: a(6) = (/ &
|
||||
-3.969683028665376e1_8, 2.209460984245205e2_8, -2.759285104469687e2_8, &
|
||||
1.383577518672690e2_8, -3.066479806614716e1_8, 2.506628277459239e0_8 /)
|
||||
real(8), parameter :: b(5) = (/ &
|
||||
-5.447609879822406e1_8, 1.615858368580409e2_8, -1.556989798598866e2_8, &
|
||||
6.680131188771972e1_8, -1.328068155288572e1_8 /)
|
||||
real(8), parameter :: c(6) = (/ &
|
||||
-7.784894002430293e-3_8, -3.223964580411365e-1_8, -2.400758277161838_8, &
|
||||
-2.549732539343734_8, 4.374664141464968_8, 2.938163982698783_8 /)
|
||||
real(8), parameter :: d(4) = (/ &
|
||||
7.784695709041462e-3_8, 3.224671290700398e-1_8, &
|
||||
2.445134137142996_8, 3.754408661907416_8 /)
|
||||
|
||||
if (p < p_low) then
|
||||
! Rational approximation for lower region.
|
||||
|
||||
q = sqrt(-TWO*log(p))
|
||||
z = (((((c(1)*q + c(2))*q + c(3))*q + c(4))*q + c(5))*q + c(6)) / &
|
||||
((((d(1)*q + d(2))*q + d(3))*q + d(4))*q + 1.)
|
||||
|
||||
elseif (p <= 1. - p_low) then
|
||||
! Rational approximation for central region
|
||||
|
||||
q = p - 0.5
|
||||
r = q*q
|
||||
z = (((((a(1)*r + a(2))*r + a(3))*r + a(4))*r + a(5))*r + a(6))*q / &
|
||||
(((((b(1)*r + b(2))*r + b(3))*r + b(4))*r + b(5))*r + 1.)
|
||||
|
||||
else
|
||||
! Rational approximation for upper region
|
||||
|
||||
q = sqrt(-2*log(1. - p))
|
||||
z = -(((((c(1)*q + c(2))*q + c(3))*q + c(4))*q + c(5))*q + c(6)) / &
|
||||
((((d(1)*q + d(2))*q + d(3))*q + d(4))*q + 1.)
|
||||
endif
|
||||
|
||||
! Refinement based on Newton's method
|
||||
#ifndef NO_F2008
|
||||
z = z - (0.5 * erfc(-z/sqrt(TWO)) - p) * sqrt(TWO*PI) * exp(0.5*z*z)
|
||||
#endif
|
||||
|
||||
end function normal_percentile
|
||||
|
||||
!===============================================================================
|
||||
! T_PERCENTILE calculates the percentile of the Student's t distribution with a
|
||||
! specified probability level and number of degrees of freedom
|
||||
!===============================================================================
|
||||
|
||||
function t_percentile(p, df) result(t)
|
||||
|
||||
real(8), intent(in) :: p ! probability level
|
||||
integer, intent(in) :: df ! degrees of freedom
|
||||
real(8) :: t ! corresponding t-value
|
||||
|
||||
real(8) :: n ! degrees of freedom as a real(8)
|
||||
real(8) :: k ! n - 2
|
||||
real(8) :: z ! percentile of normal distribution
|
||||
real(8) :: z2 ! z * z
|
||||
|
||||
if (df == 1) then
|
||||
! For one degree of freedom, the t-distribution becomes a Cauchy
|
||||
! distribution whose cdf we can invert directly
|
||||
|
||||
t = tan(PI*(p - 0.5))
|
||||
|
||||
elseif (df == 2) then
|
||||
! For two degrees of freedom, the cdf is given by 1/2 + x/(2*sqrt(x^2 +
|
||||
! 2)). This can be directly inverted to yield the solution below
|
||||
|
||||
t = TWO*sqrt(TWO)*(p - 0.5)/sqrt(ONE - 4.*(p - 0.5)**2)
|
||||
|
||||
else
|
||||
|
||||
! This approximation is from E. Olusegun George and Meenakshi Sivaram, "A
|
||||
! modification of the Fisher-Cornish approximation for the student t
|
||||
! percentiles," Communication in Statistics - Simulation and Computation,
|
||||
! 16 (4), pp. 1123-1132 (1987).
|
||||
|
||||
n = real(df,8)
|
||||
k = 1./(n - 2.)
|
||||
z = normal_percentile(p)
|
||||
z2 = z * z
|
||||
t = sqrt(n*k) * (z + (z2 - 3.)*z*k/4. + ((5.*z2 - 56.)*z2 + &
|
||||
75.)*z*k*k/96. + (((z2 - 27.)*3.*z2 + 417.)*z2 - 315.) &
|
||||
*z*k*k*k/384.)
|
||||
|
||||
end if
|
||||
|
||||
end function t_percentile
|
||||
|
||||
end module tally
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue