crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Enabled writing state points with <no_reduce> option. State point format changed 4).

Browse source
This commit is contained in:
Paul Romano 2012-10-03 18:39:33 -04:00
parent 93d80ecd35
commit c5e3007b57
7 changed files with 146 additions and 156 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/DEPENDENCIES
View file

@ -251,6 +251,7 @@ state_point.o: output.o
state_point.o: source.o
state_point.o: string.o
state_point.o: tally_header.o
state_point.o: tally.o
string.o: constants.o
string.o: error.o

2
src/constants.F90
View file

@ -12,7 +12,7 @@ module constants
! Revision numbers for binary files
integer, parameter :: REVISION_SOURCE = 1
integer, parameter :: REVISION_STATEPOINT = 3
integer, parameter :: REVISION_STATEPOINT = 4
! ============================================================================
! ADJUSTABLE PARAMETERS

31
src/hdf5_interface.F90
View file

@ -886,20 +886,6 @@ contains
dims, entropy, hdf5_err)
end if
! Write out global tallies sum and sum_sq
call hdf5_make_integer(hdf5_state_point, "n_global_tallies", &
N_GLOBAL_TALLIES)
! Write global tallies
dims(1) = N_GLOBAL_TALLIES
call h5screate_simple_f(1, dims, dspace, hdf5_err)
call h5dcreate_f(hdf5_state_point, "global_tallies", hdf5_tallyscore_t, &
dspace, dset, hdf5_err)
f_ptr = c_loc(global_tallies(1))
CALL h5dwrite_f(dset, hdf5_tallyscore_t, f_ptr, hdf5_err)
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
! Create group for tallies
call h5gcreate_f(hdf5_state_point, "tallies", tallies_group, hdf5_err)
@ -1050,10 +1036,27 @@ contains
call h5gclose_f(temp_group, hdf5_err)
end do TALLY_METADATA
! Write out global tallies sum and sum_sq
call hdf5_make_integer(hdf5_state_point, "n_global_tallies", &
N_GLOBAL_TALLIES)
! Write global tallies
dims(1) = N_GLOBAL_TALLIES
call h5screate_simple_f(1, dims, dspace, hdf5_err)
call h5dcreate_f(hdf5_state_point, "global_tallies", hdf5_tallyscore_t, &
dspace, dset, hdf5_err)
f_ptr = c_loc(global_tallies(1))
CALL h5dwrite_f(dset, hdf5_tallyscore_t, f_ptr, hdf5_err)
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
if (tallies_on) then
! Indicate that tallies are on
call hdf5_make_integer(tallies_group, "tallies_present", 1)
! Write number of realizations
call hdf5_make_integer(tallies_group, "n_realizations", n_realizations)
! Write tally sum and sum_sq
TALLY_SCORES: do i = 1, n_tallies
! Get pointer to tally

1
src/main.F90
View file

@ -24,7 +24,6 @@ program main
case (MODE_TALLIES)
! For tallies-only mode, we just skip straight to finalize_run to write out
! the tally results
n_realizations = restart_batch - n_inactive
end select
! finalize run

147
src/state_point.F90
View file

@ -7,6 +7,7 @@ module state_point
use source, only: write_source_binary
use string, only: to_str
use tally_header, only: TallyObject
use tally, only: setup_active_tallies
#ifdef MPI
use mpi
@ -66,12 +67,22 @@ contains
call write_tally_scores_nr(fh)
elseif (master) then
! Write global tallies
call MPI_FILE_WRITE(fh, N_GLOBAL_TALLIES, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
call MPI_FILE_WRITE(fh, global_tallies, N_GLOBAL_TALLIES, &
MPI_TALLYSCORE, MPI_STATUS_IGNORE, mpi_err)
if (tallies_on) then
! Indicate that tallies are on
temp = 1
call MPI_FILE_WRITE(fh, temp, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
! Write number of realizations
call MPI_FILE_WRITE(fh, n_realizations, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
! Write all tally scores
TALLY_SCORES: do i = 1, n_tallies
t => tallies(i)
@ -157,13 +168,6 @@ contains
write(UNIT_STATE) entropy(1:current_batch)
end if
! Write out global tallies sum and sum_sq
write(UNIT_STATE) N_GLOBAL_TALLIES
GLOBAL_TALLIES_LOOP: do i = 1, N_GLOBAL_TALLIES
write(UNIT_STATE) global_tallies(i) % sum
write(UNIT_STATE) global_tallies(i) % sum_sq
end do GLOBAL_TALLIES_LOOP
! Write number of meshes
write(UNIT_STATE) n_meshes
@ -236,10 +240,20 @@ contains
write(UNIT_STATE) t % score_bins
end do TALLY_METADATA
! Write out global tallies sum and sum_sq
write(UNIT_STATE) N_GLOBAL_TALLIES
GLOBAL_TALLIES_LOOP: do i = 1, N_GLOBAL_TALLIES
write(UNIT_STATE) global_tallies(i) % sum
write(UNIT_STATE) global_tallies(i) % sum_sq
end do GLOBAL_TALLIES_LOOP
if (tallies_on) then
! Indicate that tallies are on
write(UNIT_STATE) 1
! Number of realizations
write(UNIT_STATE) n_realizations
TALLY_SCORES: do i = 1, n_tallies
! Get pointer to tally
t => tallies(i)
@ -338,12 +352,6 @@ contains
MPI_STATUS_IGNORE, mpi_err)
end if
! Write out global tallies sum and sum_sq
call MPI_FILE_WRITE(fh, N_GLOBAL_TALLIES, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
call MPI_FILE_WRITE(fh, global_tallies, N_GLOBAL_TALLIES, &
MPI_TALLYSCORE, MPI_STATUS_IGNORE, mpi_err)
! Write number of meshes
call MPI_FILE_WRITE(fh, n_meshes, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
@ -455,21 +463,62 @@ contains
subroutine write_tally_scores_nr(fh)
integer, intent(in) :: fh ! file handle
integer :: i ! loop index
integer :: n ! number of filter bins
integer :: m ! number of score bins
integer :: temp ! temporary variable
integer :: n_bins ! total number of bins
real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of scores
real(8) :: global_temp(2,N_GLOBAL_TALLIES)
real(8) :: dummy ! temporary receive buffer for non-root reduces
type(TallyObject), pointer :: t => null()
! ==========================================================================
! COLLECT AND WRITE GLOBAL TALLIES
! Write number of global tallies
if (master) then
call MPI_FILE_WRITE(fh, N_GLOBAL_TALLIES, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
end if
! Copy global tallies into temporary array for reducing
n_bins = 2 * N_GLOBAL_TALLIES
global_temp(1,:) = global_tallies(:) % sum
global_temp(2,:) = global_tallies(:) % sum_sq
if (master) then
! The MPI_IN_PLACE specifier allows the master to copy values into a
! receive buffer without having a temporary variable
call MPI_REDUCE(MPI_IN_PLACE, global_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
! Write out global tallies sum and sum_sq
call MPI_FILE_WRITE(fh, global_temp, n_bins, MPI_REAL8, &
MPI_STATUS_IGNORE, mpi_err)
! Transfer values to value on master
if (current_batch == n_batches) then
global_tallies(:) % sum = global_temp(1,:)
global_tallies(:) % sum_sq = global_temp(2,:)
end if
else
! Receive buffer not significant at other processors
call MPI_REDUCE(global_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
end if
if (tallies_on) then
! Indicate that tallies are on
if (master) then
temp = 1
call MPI_FILE_WRITE(fh, temp, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
! Write number of realizations
call MPI_FILE_WRITE(fh, n_realizations, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
end if
! Write all tally scores
@ -611,17 +660,6 @@ contains
end if
if (master) then
! Read number of global tallies and make sure it matches
call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err)
if (temp(1) /= N_GLOBAL_TALLIES) then
message = "Number of global tallies does not match in state point."
call fatal_error()
end if
! Read global tally data
call MPI_FILE_READ(fh, global_tallies, N_GLOBAL_TALLIES, &
MPI_TALLYSCORE, MPI_STATUS_IGNORE, mpi_err)
! Read number of meshes
call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err)
if (temp(1) /= n_meshes) then
@ -708,6 +746,17 @@ contains
deallocate(int_array)
end do TALLY_METADATA
! Read number of global tallies and make sure it matches
call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err)
if (temp(1) /= N_GLOBAL_TALLIES) then
message = "Number of global tallies does not match in state point."
call fatal_error()
end if
! Read global tally data
call MPI_FILE_READ(fh, global_tallies, N_GLOBAL_TALLIES, &
MPI_TALLYSCORE, MPI_STATUS_IGNORE, mpi_err)
! Check if tally results are present
call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err)
@ -716,6 +765,10 @@ contains
! Read sum and sum squared
if (temp(1) == 1) then
! Read number of realizations
call MPI_FILE_READ(fh, n_realizations, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
TALLY_SCORES: do i = 1, n_tallies
n = size(tallies(i) % scores, 1) * size(tallies(i) % scores, 2)
call MPI_FILE_READ(fh, tallies(i) % scores, n, MPI_TALLYSCORE, &
@ -795,19 +848,6 @@ contains
end if
if (master) then
! Read number of global tallies and make sure it matches
read(UNIT_STATE) temp(1)
if (temp(1) /= N_GLOBAL_TALLIES) then
message = "Number of global tallies does not match in state point."
call fatal_error()
end if
! Read global tally data
do i = 1, N_GLOBAL_TALLIES
read(UNIT_STATE) global_tallies(i) % sum
read(UNIT_STATE) global_tallies(i) % sum_sq
end do
! Read number of meshes
read(UNIT_STATE) temp(1)
if (temp(1) /= n_meshes) then
@ -892,9 +932,25 @@ contains
deallocate(int_array)
end do TALLY_METADATA
! Read number of global tallies and make sure it matches
read(UNIT_STATE) temp(1)
if (temp(1) /= N_GLOBAL_TALLIES) then
message = "Number of global tallies does not match in state point."
call fatal_error()
end if
! Read global tally data
do i = 1, N_GLOBAL_TALLIES
read(UNIT_STATE) global_tallies(i) % sum
read(UNIT_STATE) global_tallies(i) % sum_sq
end do
! Read sum and sum squared
read(UNIT_STATE) temp(1)
if (temp(1) == 1) then
! Read number of realizations
read(UNIT_STATE) n_realizations
TALLY_SCORES: do i = 1, n_tallies
do k = 1, size(tallies(i) % scores, 2)
do j = 1, size(tallies(i) % scores, 1)
@ -924,6 +980,7 @@ contains
subroutine replay_batch_history
integer :: n = 0 ! number of realizations
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
@ -936,29 +993,31 @@ contains
! For criticality calculations, turn on tallies if we've reached active
! batches
if (current_batch == n_inactive) tallies_on = .true.
if (current_batch == n_inactive) then
tallies_on = .true.
call setup_active_tallies()
end if
! Add to number of realizations
if (current_batch > n_inactive) then
n_realizations = n_realizations + 1
n = n + 1
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 = temp(1) / n
if (n_realizations > 1) then
if (n > 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 - 1)
t_value = t_percentile(ONE - alpha/TWO, n - 1)
else
t_value = ONE
end if
keff_std = t_value * sqrt((temp(2)/n_realizations - keff*keff) &
/ (n_realizations - 1))
keff_std = t_value * sqrt((temp(2)/n - keff*keff)/(n - 1))
end if
else
keff = k_batch(current_batch)

84
src/tally.F90
View file

@ -1668,11 +1668,6 @@ contains
call accumulate_score(global_tallies)
end if
#ifdef MPI
! if ((.not. reduce_tallies) .and. current_batch == n_batches) &
! call reduce_tally_sums()
#endif
if (associated(curr_ptr)) nullify(curr_ptr)
end subroutine synchronize_tallies
@ -1760,83 +1755,6 @@ contains
end subroutine reduce_tally_values
#endif
!===============================================================================
! REDUCE_TALLY_SUMS gathers data from the sum and sum_sq member variables of
! TallyScores rather than the data from values. This is used if the user
! specified that tallies should not be reduced and scores are accumulated on
! every processor *before* being reduced.
! ===============================================================================
#ifdef MPI
subroutine reduce_tally_sums()
integer :: n ! number of filter bins
integer :: m ! number of score bins
integer :: n_bins ! total number of bins
real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of scores
real(8) :: global_temp(2,N_GLOBAL_TALLIES)
real(8) :: dummy ! temporary receive buffer for non-root reduces
type(TallyObject), pointer :: t => null()
type(ListInt), pointer :: curr_ptr => null()
curr_ptr => active_tallies
do while(associated(curr_ptr))
t => tallies(curr_ptr % data)
m = t % n_score_bins * t % n_nuclide_bins
n = t % n_total_bins
n_bins = m*n*2
allocate(tally_temp(2,m,n))
tally_temp(1,:,:) = t % scores(:,:) % sum
tally_temp(2,:,:) = t % scores(:,:) % sum_sq
if (master) then
! The MPI_IN_PLACE specifier allows the master to copy values into a
! receive buffer without having a temporary variable
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
! Transfer values to value on master
t % scores(:,:) % sum = tally_temp(1,:,:)
t % scores(:,:) % sum_sq = tally_temp(2,:,:)
else
! Receive buffer not significant at other processors
call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
end if
deallocate(tally_temp)
curr_ptr => curr_ptr % next
end do
n_bins = 2 * N_GLOBAL_TALLIES
global_temp(1,:) = global_tallies(:) % sum
global_temp(2,:) = global_tallies(:) % sum_sq
if (master) then
! The MPI_IN_PLACE specifier allows the master to copy values into a
! receive buffer without having a temporary variable
call MPI_REDUCE(MPI_IN_PLACE, global_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
! Transfer values to value on master
global_tallies(:) % sum = global_temp(1,:)
global_tallies(:) % sum_sq = global_temp(2,:)
else
! Receive buffer not significant at other processors
call MPI_REDUCE(global_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
end if
if (associated(curr_ptr)) nullify(curr_ptr)
end subroutine reduce_tally_sums
#endif
!===============================================================================
! WRITE_TALLIES creates an output file and writes out the mean values of all
! tallies and their standard deviations
@ -2371,7 +2289,7 @@ contains
end subroutine reset_score
!===============================================================================
! SETUP_ACTIVE_USERTALLIES
! SETUP_ACTIVE_TALLIES
!===============================================================================
subroutine setup_active_tallies()

36
src/utils/statepoint.py
View file

@ -117,11 +117,6 @@ class StatePoint(BinaryFile):
self.k_batch = self._get_double(self.current_batch)
self.entropy = self._get_double(self.current_batch)
# Read global tallies
n_global_tallies = self._get_int()[0]
self.global_tallies = np.array(self._get_double(2*n_global_tallies))
self.global_tallies.shape = (n_global_tallies, 2)
# Read number of meshes
n_meshes = self._get_int()[0]
@ -191,11 +186,18 @@ class StatePoint(BinaryFile):
if not self._metadata:
self._read_metadata()
# Read global tallies
n_global_tallies = self._get_int()[0]
self.global_tallies = np.array(self._get_double(2*n_global_tallies))
self.global_tallies.shape = (n_global_tallies, 2)
# Flag indicating if tallies are present
tallies_present = self._get_int()[0]
# Read tally results
if tallies_present:
self.n_realizations = self._get_int()[0]
for t in self.tallies:
n = t.n_score_bins * t.n_filter_bins
t.values = np.array(self._get_double(2*n))
@ -223,10 +225,7 @@ class StatePoint(BinaryFile):
"""Calculates confidence intervals for each tally bin."""
# Determine number of realizations
if self.run_mode == 2:
n = self.current_batch - self.n_inactive
else:
n = self.current_batch
n = self.n_realizations
# Determine significance level and percentile for two-sided CI
alpha = 1 - confidence
@ -243,11 +242,22 @@ class StatePoint(BinaryFile):
"""
# Determine number of realizations
if self.run_mode == 2:
n = self.current_batch - self.n_inactive
else:
n = self.current_batch
n = self.n_realizations
# Global tallies
for i in range(len(self.global_tallies)):
# Get sum and sum of squares
s, s2 = self.global_tallies[i]
# Calculate sample mean and replace value
s /= n
self.global_tallies[i,0] = s
# Calculate standard deviation
if s != 0.0:
self.global_tallies[i,1] = t_value*sqrt((s2/n - s*s)/(n-1))
# Regular tallies
for t in self.tallies:
for i in range(t.values.shape[0]):
for j in range(t.values.shape[1]):