Merge branch 'combined'

This commit is contained in:
Paul Romano 2013-02-21 16:26:10 -05:00
commit 919afc3cae
9 changed files with 430 additions and 15 deletions

View file

@ -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*

View file

@ -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
@ -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
! ============================================================================

View file

@ -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,102 @@ 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
! 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
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

View file

@ -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_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_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength
real(8) :: k_combined(2) ! combined best estimate of k-effective
! Shannon entropy
logical :: entropy_on = .false.

View file

@ -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

View file

@ -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
@ -1429,6 +1432,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,*)

View file

@ -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

View file

@ -1742,6 +1742,9 @@ contains
subroutine synchronize_tallies()
integer :: i
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
! Combine tally results onto master process
@ -1768,6 +1771,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_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_col_abs = k_col_abs + k_col * k_abs
k_col_tra = k_col_tra + k_col * k_tra
k_abs_tra = k_abs_tra + k_abs * k_tra
end if
end if
! Accumulate results for global tallies

View file

@ -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]