From c5e3007b57514c67ca8cba3a52a031e887bea288 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 3 Oct 2012 18:39:33 -0400 Subject: [PATCH] Enabled writing state points with option. State point format changed 4). --- src/DEPENDENCIES | 1 + src/constants.F90 | 2 +- src/hdf5_interface.F90 | 31 +++++---- src/main.F90 | 1 - src/state_point.F90 | 147 ++++++++++++++++++++++++++++------------ src/tally.F90 | 84 +---------------------- src/utils/statepoint.py | 36 ++++++---- 7 files changed, 146 insertions(+), 156 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index f31bc25f42..4188fa22a1 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -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 diff --git a/src/constants.F90 b/src/constants.F90 index ccede3e202..f69df99239 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 = 3 + integer, parameter :: REVISION_STATEPOINT = 4 ! ============================================================================ ! ADJUSTABLE PARAMETERS diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index b6406a15fa..4942a025c9 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -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 diff --git a/src/main.F90 b/src/main.F90 index ad0da8a182..dc72e5f8c5 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -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 diff --git a/src/state_point.F90 b/src/state_point.F90 index 7dca5c8967..e99358231c 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -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) diff --git a/src/tally.F90 b/src/tally.F90 index 0ed9b5b04f..083948af55 100755 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -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() diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 5ef88bc6d1..a1c12a880f 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -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]):