diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index ed790f4c8b..48ec5f8220 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -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 diff --git a/src/OBJECTS b/src/OBJECTS index 7b369d3dba..66a05ef4f2 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -24,6 +24,7 @@ interpolation.o \ input_xml.o \ main.o \ material_header.o \ +math.o \ mesh_header.o \ mesh.o \ output.o \ diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 5a47411c91..b2729df997 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -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 diff --git a/src/math.F90 b/src/math.F90 new file mode 100644 index 0000000000..e60528e66a --- /dev/null +++ b/src/math.F90 @@ -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 diff --git a/src/output.F90 b/src/output.F90 index 063787bf12..b672b494fd 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/src/state_point.F90 b/src/state_point.F90 index 1cfd93b34d..debf55ec9b 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -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) diff --git a/src/tally.F90 b/src/tally.F90 index d261b37948..c4836e85f5 100644 --- a/src/tally.F90 +++ b/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