diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index e8b8f8b90..7551f4293 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -224,12 +224,12 @@ if run_mode == 'k-eigenvalue': Tallying moment orders for Legendre and spherical harmonic tally expansions (*e.g.*, 'P2', 'Y1,2', etc.). -**/tallies/tally /results** (Compound type) +**/tallies/tally /results** (*double[][][2]*) - Accumulated sum and sum-of-squares for each bin of the i-th tally. This is a - two-dimensional array, the first dimension of which represents combinations - of filter bins and the second dimensions of which represents scoring - bins. Each element of the array has fields 'sum' and 'sum_sq'. + Accumulated sum and sum-of-squares for each bin of the i-th tally. The first + dimension represents combinations of filter bins, the second dimensions + represents scoring bins, and the third dimension has two entries for the sum + and the sum-of-squares. **/source_present** (*int*) diff --git a/openmc/particle_restart.py b/openmc/particle_restart.py index a0410298b..cb1326586 100644 --- a/openmc/particle_restart.py +++ b/openmc/particle_restart.py @@ -1,3 +1,5 @@ +import h5py + class Particle(object): """Information used to restart a specific particle that caused a simulation to fail. @@ -33,12 +35,6 @@ class Particle(object): """ def __init__(self, filename): - import h5py - if h5py.__version__ == '2.6.0': - raise ImportError("h5py 2.6.0 has a known bug which makes it " - "incompatible with OpenMC's HDF5 files. " - "Please switch to a different version.") - self._f = h5py.File(filename, 'r') # Ensure filetype and revision are correct diff --git a/openmc/statepoint.py b/openmc/statepoint.py index aaa5c0b84..e05b18ebe 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -5,6 +5,7 @@ import warnings import glob import numpy as np +import h5py import openmc import openmc.checkvalue as cv @@ -104,12 +105,6 @@ class StatePoint(object): """ def __init__(self, filename, autolink=True): - import h5py - if h5py.__version__ == '2.6.0': - raise ImportError("h5py 2.6.0 has a known bug which makes it " - "incompatible with OpenMC's HDF5 files. " - "Please switch to a different version.") - self._f = h5py.File(filename, 'r') # Ensure filetype and revision are correct @@ -209,13 +204,13 @@ class StatePoint(object): def global_tallies(self): if self._global_tallies is None: data = self._f['global_tallies'].value - gt = np.zeros_like(data, dtype=[ + gt = np.zeros(data.shape[0], dtype=[ ('name', 'a14'), ('sum', 'f8'), ('sum_sq', 'f8'), ('mean', 'f8'), ('std_dev', 'f8')]) gt['name'] = ['k-collision', 'k-absorption', 'k-tracklength', 'leakage'] - gt['sum'] = data['sum'] - gt['sum_sq'] = data['sum_sq'] + gt['sum'] = data[:,1] + gt['sum_sq'] = data[:,2] # Calculate mean and sample standard deviation of mean n = self.n_realizations diff --git a/openmc/summary.py b/openmc/summary.py index fd1689eb2..37509ec37 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -2,6 +2,7 @@ from collections import Iterable import re import numpy as np +import h5py import openmc from openmc.region import Region @@ -23,15 +24,6 @@ class Summary(object): """ def __init__(self, filename): - # A user may not have h5py, but they can still use the rest of the - # Python API so we'll only try to import h5py if the user actually inits - # a Summary object. - import h5py - if h5py.__version__ == '2.6.0': - raise ImportError("h5py 2.6.0 has a known bug which makes it " - "incompatible with OpenMC's HDF5 files. " - "Please switch to a different version.") - openmc.reset_auto_ids() if not filename.endswith(('.h5', '.hdf5')): diff --git a/openmc/tallies.py b/openmc/tallies.py index 604713145..28da08bad 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -13,6 +13,7 @@ from xml.etree import ElementTree as ET from six import string_types import numpy as np +import h5py import openmc import openmc.checkvalue as cv @@ -269,20 +270,14 @@ class Tally(object): return None if not self._results_read: - import h5py - if h5py.__version__ == '2.6.0': - raise ImportError("h5py 2.6.0 has a known bug which makes it " - "incompatible with OpenMC's HDF5 files. " - "Please switch to a different version.") - # Open the HDF5 statepoint file f = h5py.File(self._sp_filename, 'r') # Extract Tally data from the file data = f['tallies/tally {0}/results'.format( self.id)].value - sum = data['sum'] - sum_sq = data['sum_sq'] + sum = data[:,:,0] + sum_sq = data[:,:,1] # Reshape the results arrays sum = np.reshape(sum, self.shape) @@ -1726,8 +1721,6 @@ class Tally(object): # HDF5 binary file if format == 'hdf5': - import h5py - filename = directory + '/' + filename + '.h5' if append: diff --git a/openmc/volume.py b/openmc/volume.py index 2a2a92f61..b1391c687 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -5,6 +5,7 @@ from warnings import warn import numpy as np import pandas as pd +import h5py import openmc import openmc.checkvalue as cv @@ -187,8 +188,6 @@ class VolumeCalculation(object): Results of the stochastic volume calculation """ - import h5py - with h5py.File(filename, 'r') as f: domain_type = f.attrs['domain_type'].decode() samples = f.attrs['samples'] diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index e1bf1f9c3..e5ed5df70 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -164,7 +164,7 @@ contains * t%stride) + 1 ! Get flux - flux = t % results(1,score_index) % sum + flux = t % results(RESULT_SUM,1,score_index) cmfd % flux(h,i,j,k) = flux ! Detect zero flux, abort if located @@ -175,10 +175,10 @@ contains end if ! Get total rr and convert to total xs - cmfd % totalxs(h,i,j,k) = t % results(2,score_index) % sum / flux + cmfd % totalxs(h,i,j,k) = t % results(RESULT_SUM,2,score_index) / flux ! Get p1 scatter rr and convert to p1 scatter xs - cmfd % p1scattxs(h,i,j,k) = t % results(3,score_index) % sum / flux + cmfd % p1scattxs(h,i,j,k) = t % results(RESULT_SUM,3,score_index) / flux ! Calculate diffusion coefficient cmfd % diffcof(h,i,j,k) = ONE/(3.0_8*(cmfd % totalxs(h,i,j,k) - & @@ -211,19 +211,18 @@ contains * t%stride) + 1 ! Get scattering - cmfd % scattxs(h,g,i,j,k) = t % results(1,score_index) % sum /& + cmfd % scattxs(h,g,i,j,k) = t % results(RESULT_SUM,1,score_index) /& cmfd % flux(h,i,j,k) ! Get nu-fission - cmfd % nfissxs(h,g,i,j,k) = t % results(2,score_index) % sum /& + cmfd % nfissxs(h,g,i,j,k) = t % results(RESULT_SUM,2,score_index) /& cmfd % flux(h,i,j,k) ! Bank source cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + & - t % results(2,score_index) % sum + t % results(RESULT_SUM,2,score_index) cmfd % keff_bal = cmfd % keff_bal + & - t % results(2,score_index) % sum / & - dble(t % n_realizations) + t % results(RESULT_SUM,2,score_index) / t % n_realizations end do INGROUP @@ -243,67 +242,67 @@ contains matching_bins(i_filter_surf) = OUT_LEFT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(1,h,i,j,k) = t % results(RESULT_SUM,1,score_index) matching_bins(i_filter_surf) = IN_LEFT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(2,h,i,j,k) = t % results(RESULT_SUM,1,score_index) ! Right surface matching_bins(i_filter_surf) = IN_RIGHT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(3,h,i,j,k) = t % results(RESULT_SUM,1,score_index) matching_bins(i_filter_surf) = OUT_RIGHT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(4,h,i,j,k) = t % results(RESULT_SUM,1,score_index) ! Back surface matching_bins(i_filter_surf) = OUT_BACK score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(5,h,i,j,k) = t % results(RESULT_SUM,1,score_index) matching_bins(i_filter_surf) = IN_BACK score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(6,h,i,j,k) = t % results(RESULT_SUM,1,score_index) ! Front surface matching_bins(i_filter_surf) = IN_FRONT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(7,h,i,j,k) = t % results(RESULT_SUM,1,score_index) matching_bins(i_filter_surf) = OUT_FRONT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(8,h,i,j,k) = t % results(RESULT_SUM,1,score_index) ! Bottom surface matching_bins(i_filter_surf) = OUT_BOTTOM score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(9,h,i,j,k) = t % results(RESULT_SUM,1,score_index) matching_bins(i_filter_surf) = IN_BOTTOM score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(10,h,i,j,k) = t % results(RESULT_SUM,1,score_index) ! Top surface matching_bins(i_filter_surf) = IN_TOP score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(11,h,i,j,k) = t % results(RESULT_SUM,1,score_index) matching_bins(i_filter_surf) = OUT_TOP score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(12,h,i,j,k) = t % results(RESULT_SUM,1,score_index) end if TALLY diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index d2631254b..a5d92002b 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -365,22 +365,18 @@ contains subroutine cmfd_tally_reset() - use global, only: n_cmfd_tallies, cmfd_tallies + use global, only: cmfd_tallies use output, only: write_message - use tally, only: reset_result integer :: i ! loop counter ! Print message call write_message("CMFD tallies reset", 7) - ! Begin loop around CMFD tallies - do i = 1, n_cmfd_tallies - - ! Reset that tally + ! Reset CMFD tallies + do i = 1, size(cmfd_tallies) cmfd_tallies(i) % n_realizations = 0 - call reset_result(cmfd_tallies(i) % results) - + cmfd_tallies(i) % results(:,:,:) = ZERO end do end subroutine cmfd_tally_reset diff --git a/src/constants.F90 b/src/constants.F90 index ac157d7ad..f2a9fb6f1 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -269,6 +269,12 @@ module constants ! ============================================================================ ! TALLY-RELATED CONSTANTS + ! Tally result entries + integer, parameter :: & + RESULT_VALUE = 1, & + RESULT_SUM = 2, & + RESULT_SUM_SQ = 3 + ! Tally type integer, parameter :: & TALLY_VOLUME = 1, & diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 713bbc351..05c8190a6 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -374,7 +374,7 @@ contains subroutine calculate_generation_keff() ! Get keff for this generation by subtracting off the starting value - keff_generation = global_tallies(K_TRACKLENGTH) % value - keff_generation + keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - keff_generation #ifdef MPI ! Combine values across all processors @@ -466,14 +466,14 @@ contains 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 - & + kv(1) = global_tallies(RESULT_SUM, K_COLLISION) / n + kv(2) = global_tallies(RESULT_SUM, K_ABSORPTION) / n + kv(3) = global_tallies(RESULT_SUM, K_TRACKLENGTH) / n + cov(1,1) = (global_tallies(RESULT_SUM_SQ, K_COLLISION) - & n * kv(1) * kv(1)) / (n - 1) - cov(2,2) = (global_tallies(K_ABSORPTION) % sum_sq - & + cov(2,2) = (global_tallies(RESULT_SUM_SQ, K_ABSORPTION) - & n * kv(2) * kv(2)) / (n - 1) - cov(3,3) = (global_tallies(K_TRACKLENGTH) % sum_sq - & + cov(3,3) = (global_tallies(RESULT_SUM_SQ, K_TRACKLENGTH) - & n * kv(3) * kv(3)) / (n - 1) ! Calculate covariances based on sums with Bessel's correction diff --git a/src/finalize.F90 b/src/finalize.F90 index 86100d195..83a5ac5fa 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -9,7 +9,7 @@ module finalize use message_passing #endif - use hdf5_interface, only: hdf5_bank_t, hdf5_tallyresult_t + use hdf5_interface, only: hdf5_bank_t use hdf5, only: h5tclose_f, h5close_f implicit none @@ -53,7 +53,6 @@ contains call free_memory() ! Release compound datatypes - call h5tclose_f(hdf5_tallyresult_t, hdf5_err) call h5tclose_f(hdf5_bank_t, hdf5_err) ! Close FORTRAN interface. @@ -62,7 +61,6 @@ contains #ifdef MPI ! Free all MPI types call MPI_TYPE_FREE(MPI_BANK, mpi_err) - call MPI_TYPE_FREE(MPI_TALLYRESULT, mpi_err) ! If MPI is in use and enabled, terminate it call MPI_FINALIZE(mpi_err) diff --git a/src/global.F90 b/src/global.F90 index 8111fd514..eb4b33d7d 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -1,5 +1,11 @@ module global + use, intrinsic :: ISO_C_BINDING + +#ifdef MPIF08 + use mpi_f08 +#endif + use bank_header, only: Bank use cmfd_header use constants @@ -14,15 +20,11 @@ module global use set_header, only: SetInt use surface_header, only: SurfaceContainer use source_header, only: SourceDistribution - use tally_header, only: TallyObject, TallyResult + use tally_header, only: TallyObject use trigger_header, only: KTrigger use timer_header, only: Timer use volume_header, only: VolumeCalculation -#ifdef MPIF08 - use mpi_f08 -#endif - implicit none ! ============================================================================ @@ -164,7 +166,7 @@ module global ! 3) track-length estimate of k-eff ! 4) leakage fraction - type(TallyResult), allocatable, target :: global_tallies(:) + real(C_DOUBLE), allocatable, target :: global_tallies(:,:) ! It is possible to protect accumulate operations on global tallies by using ! an atomic update. However, when multiple threads accumulate to the same @@ -272,10 +274,8 @@ module global integer :: mpi_err ! MPI error code #ifdef MPIF08 type(MPI_Datatype) :: MPI_BANK - type(MPI_Datatype) :: MPI_TALLYRESULT #else integer :: MPI_BANK ! MPI datatype for fission bank - integer :: MPI_TALLYRESULT ! MPI datatype for TallyResult #endif #ifdef _OPENMP diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 2a0711954..5e14605ae 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -16,7 +16,6 @@ module hdf5_interface use h5lt use error, only: fatal_error - use tally_header, only: TallyResult #ifdef PHDF5 use message_passing, only: MPI_COMM_WORLD, MPI_INFO_NULL #endif @@ -24,7 +23,6 @@ module hdf5_interface implicit none private - integer(HID_T), public :: hdf5_tallyresult_t ! Compound type for TallyResult integer(HID_T), public :: hdf5_bank_t ! Compound type for Bank integer(HID_T), public :: hdf5_integer8_t ! type for integer(8) @@ -42,8 +40,6 @@ module hdf5_interface module procedure write_long module procedure write_string module procedure write_string_1D - module procedure write_tally_result_1D - module procedure write_tally_result_2D end interface write_dataset interface read_dataset @@ -60,8 +56,6 @@ module hdf5_interface module procedure read_long module procedure read_string module procedure read_string_1D - module procedure read_tally_result_1D - module procedure read_tally_result_2D module procedure read_complex_2D end interface read_dataset @@ -2063,130 +2057,6 @@ contains call h5ltset_attribute_string_f(group_id, var, attr_type, attr_str, hdf5_err) end subroutine write_attribute_string -!=============================================================================== -! WRITE_TALLY_RESULT writes an OpenMC TallyResult type -!=============================================================================== - - subroutine write_tally_result_1D(group_id, name, buffer) - integer(HID_T), intent(in) :: group_id - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(in), target :: buffer(:) ! data to write - - integer(HSIZE_T) :: dims(1) - - dims(:) = shape(buffer) - call write_tally_result_1D_explicit(group_id, dims, name, buffer) - end subroutine write_tally_result_1D - - subroutine write_tally_result_1D_explicit(group_id, dims, name, buffer) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(1) - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(in), target :: buffer(dims(1)) - - integer :: hdf5_err - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr - - call h5screate_simple_f(1, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), hdf5_tallyresult_t, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - call h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_tally_result_1D_explicit - - subroutine write_tally_result_2D(group_id, name, buffer) - integer(HID_T), intent(in) :: group_id - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(in), target :: buffer(:,:) ! data to write - - integer(HSIZE_T) :: dims(2) - - dims(:) = shape(buffer) - call write_tally_result_2D_explicit(group_id, dims, name, buffer) - end subroutine write_tally_result_2D - - subroutine write_tally_result_2D_explicit(group_id, dims, name, buffer) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(2) - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(in), target :: buffer(dims(1),dims(2)) - - integer :: hdf5_err - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr - - call h5screate_simple_f(2, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), hdf5_tallyresult_t, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - call h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_tally_result_2D_explicit - -!=============================================================================== -! READ_TALLY_RESULT reads OpenMC TallyResult data -!=============================================================================== - - subroutine read_tally_result_1D(group_id, name, buffer) - integer(HID_T), intent(in) :: group_id - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(inout), target :: buffer(:) ! read data here - - integer(HSIZE_T) :: dims(1) - - dims(:) = shape(buffer) - call read_tally_result_1D_explicit(group_id, dims, name, buffer) - end subroutine read_tally_result_1D - - subroutine read_tally_result_1D_explicit(group_id, dims, name, buffer) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(1) - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(inout), target :: buffer(dims(1)) - - integer :: hdf5_err - integer(HID_T) :: dset ! data set handle - type(c_ptr) :: f_ptr - - call h5dopen_f(group_id, trim(name), dset, hdf5_err) - f_ptr = c_loc(buffer) - call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - call h5dclose_f(dset, hdf5_err) - end subroutine read_tally_result_1D_explicit - - subroutine read_tally_result_2D(group_id, name, buffer) - integer(HID_T), intent(in) :: group_id - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(inout), target :: buffer(:,:) - - integer(HSIZE_T) :: dims(2) - - dims(:) = shape(buffer) - call read_tally_result_2D_explicit(group_id, dims, name, buffer) - end subroutine read_tally_result_2D - - subroutine read_tally_result_2D_explicit(group_id, dims, name, buffer) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(2) - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(inout), target :: buffer(dims(1),dims(2)) - - integer :: hdf5_err - integer(HID_T) :: dset ! data set handle - type(c_ptr) :: f_ptr - - call h5dopen_f(group_id, trim(name), dset, hdf5_err) - f_ptr = c_loc(buffer) - call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - call h5dclose_f(dset, hdf5_err) - end subroutine read_tally_result_2D_explicit - subroutine read_attribute_double(buffer, obj_id, name) real(8), intent(inout), target :: buffer integer(HID_T), intent(in) :: obj_id diff --git a/src/initialize.F90 b/src/initialize.F90 index 4fe6f8769..692319acb 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -12,7 +12,7 @@ module initialize &BASE_UNIVERSE use global use hdf5_interface, only: file_open, read_dataset, file_close, hdf5_bank_t,& - hdf5_tallyresult_t, hdf5_integer8_t + hdf5_integer8_t use input_xml, only: read_input_xml, cells_in_univ_dict, read_plots_xml use material_header, only: Material use mgxs_data, only: read_mgxs, create_macro_xs @@ -22,7 +22,7 @@ module initialize use state_point, only: load_state_point use string, only: to_str, starts_with, ends_with, str_to_int use summary, only: write_summary - use tally_header, only: TallyObject, TallyResult + use tally_header, only: TallyObject use tally_initialize,only: configure_tallies use tally_filter use tally, only: init_tally_routines @@ -169,21 +169,11 @@ contains integer :: bank_blocks(5) ! Count for each datatype #ifdef MPIF08 type(MPI_Datatype) :: bank_types(5) - type(MPI_Datatype) :: result_types(1) - type(MPI_Datatype) :: temp_type #else integer :: bank_types(5) ! Datatypes - integer :: result_types(1) ! Datatypes - integer :: temp_type ! temporary derived type #endif integer(MPI_ADDRESS_KIND) :: bank_disp(5) ! Displacements - integer :: result_blocks(1) ! Count for each datatype - integer(MPI_ADDRESS_KIND) :: result_disp(1) ! Displacements - integer(MPI_ADDRESS_KIND) :: result_base_disp ! Base displacement - integer(MPI_ADDRESS_KIND) :: lower_bound ! Lower bound for TallyResult - integer(MPI_ADDRESS_KIND) :: extent ! Extent for TallyResult type(Bank) :: b - type(TallyResult) :: tr ! Indicate that MPI is turned on mpi_enabled = .true. @@ -222,34 +212,6 @@ contains bank_types, MPI_BANK, mpi_err) call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) - ! ========================================================================== - ! CREATE MPI_TALLYRESULT TYPE - - ! Determine displacements for MPI_BANK type - call MPI_GET_ADDRESS(tr%value, result_base_disp, mpi_err) - call MPI_GET_ADDRESS(tr%sum, result_disp(1), mpi_err) - - ! Adjust displacements - result_disp = result_disp - result_base_disp - - ! Define temporary type for TallyResult - result_blocks = (/ 2 /) - result_types = (/ MPI_REAL8 /) - call MPI_TYPE_CREATE_STRUCT(1, result_blocks, result_disp, result_types, & - temp_type, mpi_err) - - ! Adjust lower-bound and extent of type for tally score - lower_bound = 0 - extent = result_disp(1) + 16 - call MPI_TYPE_CREATE_RESIZED(temp_type, lower_bound, extent, & - MPI_TALLYRESULT, mpi_err) - - ! Commit derived type for tally scores - call MPI_TYPE_COMMIT(MPI_TALLYRESULT, mpi_err) - - ! Free temporary MPI type - call MPI_TYPE_FREE(temp_type, mpi_err) - end subroutine initialize_mpi #endif @@ -259,7 +221,6 @@ contains subroutine hdf5_initialize() - type(TallyResult), target :: tmp(2) ! temporary TallyResult type(Bank), target :: tmpb(2) ! temporary Bank integer :: hdf5_err integer(HID_T) :: coordinates_t ! HDF5 type for 3 reals @@ -268,14 +229,6 @@ contains ! Initialize FORTRAN interface. call h5open_f(hdf5_err) - ! Create the compound datatype for TallyResult - call h5tcreate_f(H5T_COMPOUND_F, h5offsetof(c_loc(tmp(1)), & - c_loc(tmp(2))), hdf5_tallyresult_t, hdf5_err) - call h5tinsert_f(hdf5_tallyresult_t, "sum", h5offsetof(c_loc(tmp(1)), & - c_loc(tmp(1)%sum)), H5T_NATIVE_DOUBLE, hdf5_err) - call h5tinsert_f(hdf5_tallyresult_t, "sum_sq", h5offsetof(c_loc(tmp(1)), & - c_loc(tmp(1)%sum_sq)), H5T_NATIVE_DOUBLE, hdf5_err) - ! Create compound type for xyz and uvw call h5tarray_create_f(H5T_NATIVE_DOUBLE, 1, dims, coordinates_t, hdf5_err) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index ff530c618..53cb11f98 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2039,7 +2039,6 @@ contains type(Library), allocatable :: libraries(:) type(VectorReal), allocatable :: nuc_temps(:) ! List of T to read for each nuclide type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b) - character(MAX_LINE_LEN) :: temp_str real(8), allocatable :: material_temps(:) logical :: file_exists character(MAX_FILE_LEN) :: env_variable diff --git a/src/math.F90 b/src/math.F90 index 3ab8195aa..11ace6976 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -1,8 +1,9 @@ module math + use, intrinsic :: ISO_C_BINDING + use constants use random_lcg, only: prn - use ISO_C_BINDING implicit none diff --git a/src/output.F90 b/src/output.F90 index 8127725a6..6a8ba23b5 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -611,7 +611,7 @@ contains t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1) ! Adjust sum_sq - global_tallies(:) % sum_sq = t_value * global_tallies(:) % sum_sq + global_tallies(RESULT_SUM_SQ,:) = t_value * global_tallies(RESULT_SUM_SQ,:) ! Adjust combined estimator if (n_realizations > 3) then @@ -623,26 +623,26 @@ contains ! write global tallies if (n_realizations > 1) then if (run_mode == MODE_EIGENVALUE) then - write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) & - % sum, global_tallies(K_COLLISION) % sum_sq - write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) & - % 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) "k-effective (Collision)", global_tallies(RESULT_SUM, & + K_COLLISION), global_tallies(RESULT_SUM_SQ, K_COLLISION) + write(ou,102) "k-effective (Track-length)", global_tallies(RESULT_SUM, & + K_TRACKLENGTH), global_tallies(RESULT_SUM_SQ, K_TRACKLENGTH) + write(ou,102) "k-effective (Absorption)", global_tallies(RESULT_SUM, & + K_ABSORPTION), global_tallies(RESULT_SUM_SQ, K_ABSORPTION) if (n_realizations > 3) write(ou,102) "Combined k-effective", k_combined end if - write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, & - global_tallies(LEAKAGE) % sum_sq + write(ou,102) "Leakage Fraction", global_tallies(RESULT_SUM, LEAKAGE), & + global_tallies(RESULT_SUM_SQ, LEAKAGE) else if (master) call warning("Could not compute uncertainties -- only one & &active batch simulated!") if (run_mode == MODE_EIGENVALUE) then - write(ou,103) "k-effective (Collision)", global_tallies(K_COLLISION) % sum - write(ou,103) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum - write(ou,103) "k-effective (Absorption)", global_tallies(K_ABSORPTION) % sum + write(ou,103) "k-effective (Collision)", global_tallies(RESULT_SUM, K_COLLISION) + write(ou,103) "k-effective (Track-length)", global_tallies(RESULT_SUM, K_TRACKLENGTH) + write(ou,103) "k-effective (Absorption)", global_tallies(RESULT_SUM, K_ABSORPTION) end if - write(ou,103) "Leakage Fraction", global_tallies(LEAKAGE) % sum + write(ou,103) "Leakage Fraction", global_tallies(RESULT_SUM, LEAKAGE) end if write(ou,*) @@ -765,7 +765,7 @@ contains end if ! Multiply uncertainty by t-value - t % results % sum_sq = t_value * t % results % sum_sq + t % results(RESULT_SUM_SQ,:,:) = t_value * t % results(RESULT_SUM_SQ,:,:) end if ! Write header block @@ -876,8 +876,8 @@ contains score_names(abs(t % score_bins(k))) write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,score_index,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,score_index,filter_index))) case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN) score_index = score_index - 1 do n_order = 0, t % moment_order(k) @@ -886,9 +886,8 @@ contains score_names(abs(t % score_bins(k))) write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) & - % sum_sq)) + to_str(t % results(RESULT_SUM,score_index,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,score_index,filter_index))) end do k = k + t % moment_order(k) case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, & @@ -902,9 +901,9 @@ contains // score_names(abs(t % score_bins(k))) write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index)& - % sum_sq)) + to_str(t % results(RESULT_SUM,score_index,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,score_index,& + filter_index))) end do end do k = k + (t % moment_order(k) + 1)**2 - 1 @@ -916,8 +915,8 @@ contains end if write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,score_index,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,score_index,filter_index))) end select end do indent = indent - 2 @@ -1020,16 +1019,16 @@ contains * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current on Left", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) matching_bins(i_filter_surf) = IN_LEFT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current on Left", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) ! Right Surface matching_bins(i_filter_surf) = OUT_RIGHT @@ -1037,16 +1036,16 @@ contains * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current on Right", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current on Right", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) if (n_dim >= 2) then @@ -1056,16 +1055,16 @@ contains * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current on Back", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) matching_bins(i_filter_surf) = IN_BACK filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current on Back", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) ! Front Surface matching_bins(i_filter_surf) = OUT_FRONT @@ -1073,16 +1072,16 @@ contains * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Net Current on Front", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Net Current on Front", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) end if if (n_dim == 3) then @@ -1092,16 +1091,16 @@ contains * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current on Bottom", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) matching_bins(i_filter_surf) = IN_BOTTOM filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current on Bottom", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) ! Top Surface matching_bins(i_filter_surf) = OUT_TOP @@ -1109,16 +1108,16 @@ contains * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current on Top", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) matching_bins(i_filter_surf) = IN_TOP filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current on Top", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) end if end do end do diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 776c7905a..082b337a4 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -119,16 +119,16 @@ contains ! Score implicit absorption estimate of keff !$omp atomic - global_tallies(K_ABSORPTION) % value = & - global_tallies(K_ABSORPTION) % value + p % absorb_wgt * & + global_tallies(RESULT_VALUE, K_ABSORPTION) = & + global_tallies(RESULT_VALUE, K_ABSORPTION) + p % absorb_wgt * & material_xs % nu_fission / material_xs % absorption else ! See if disappearance reaction happens if (material_xs % absorption > prn() * material_xs % total) then ! Score absorption estimate of keff !$omp atomic - global_tallies(K_ABSORPTION) % value = & - global_tallies(K_ABSORPTION) % value + p % wgt * & + global_tallies(RESULT_VALUE, K_ABSORPTION) = & + global_tallies(RESULT_VALUE, K_ABSORPTION) + p % wgt * & material_xs % nu_fission / material_xs % absorption p % alive = .false. diff --git a/src/simulation.F90 b/src/simulation.F90 index f59ce371c..ef055a549 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -20,8 +20,7 @@ module simulation use source, only: initialize_source, sample_external_source use state_point, only: write_state_point, write_source_point use string, only: to_str - use tally, only: synchronize_tallies, setup_active_usertallies, & - reset_result + use tally, only: synchronize_tallies, setup_active_usertallies use trigger, only: check_triggers use tracking, only: transport use volume_calc, only: run_volume_calculations @@ -220,7 +219,7 @@ contains if (ufs) call count_source_for_ufs() ! Store current value of tracklength k - keff_generation = global_tallies(K_TRACKLENGTH) % value + keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) end if end subroutine initialize_generation @@ -237,24 +236,24 @@ contains !$omp parallel !$omp critical if (run_mode == MODE_EIGENVALUE) then - global_tallies(K_COLLISION) % value = & - global_tallies(K_COLLISION) % value + global_tally_collision - global_tallies(K_ABSORPTION) % value = & - global_tallies(K_ABSORPTION) % value + global_tally_absorption - global_tallies(K_TRACKLENGTH) % value = & - global_tallies(K_TRACKLENGTH) % value + global_tally_tracklength + global_tallies(RESULT_VALUE, K_COLLISION) = & + global_tallies(RESULT_VALUE, K_COLLISION) + global_tally_collision + global_tallies(RESULT_VALUE, K_ABSORPTION) = & + global_tallies(RESULT_VALUE, K_ABSORPTION) + global_tally_absorption + global_tallies(RESULT_VALUE, K_TRACKLENGTH) = & + global_tallies(RESULT_VALUE, K_TRACKLENGTH) + global_tally_tracklength end if - global_tallies(LEAKAGE) % value = & - global_tallies(LEAKAGE) % value + global_tally_leakage + global_tallies(RESULT_VALUE, LEAKAGE) = & + global_tallies(RESULT_VALUE, LEAKAGE) + global_tally_leakage !$omp end critical ! reset private tallies if (run_mode == MODE_EIGENVALUE) then - global_tally_collision = 0 - global_tally_absorption = 0 - global_tally_tracklength = 0 + global_tally_collision = ZERO + global_tally_absorption = ZERO + global_tally_tracklength = ZERO end if - global_tally_leakage = 0 + global_tally_leakage = ZERO !$omp end parallel if (run_mode == MODE_EIGENVALUE) then @@ -302,7 +301,7 @@ contains ! Reset global tally results if (.not. active_batches) then - call reset_result(global_tallies) + global_tallies(:,:) = ZERO n_realizations = 0 end if diff --git a/src/state_point.F90 b/src/state_point.F90 index 6b7a467d4..f9616f228 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -353,7 +353,7 @@ contains ! Write sum and sum_sq for each bin tally_group = open_group(tallies_group, "tally " & // to_str(tally % id)) - call write_dataset(tally_group, "results", tally % results) + call tally % write_results_hdf5(tally_group) call close_group(tally_group) end do TALLY_RESULTS @@ -481,7 +481,7 @@ contains integer :: n_bins ! total number of bins integer(HID_T) :: tallies_group, tally_group real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of results - real(8), target :: global_temp(2,N_GLOBAL_TALLIES) + real(8), target :: global_temp(3,N_GLOBAL_TALLIES) #ifdef MPI real(8) :: dummy ! temporary receive buffer for non-root reduces #endif @@ -489,7 +489,7 @@ contains type(ElemKeyValueII), pointer :: current type(ElemKeyValueII), pointer :: next type(TallyObject), pointer :: tally - type(TallyResult), allocatable :: tallyresult_temp(:,:) + type(TallyObject) :: dummy_tally ! ========================================================================== ! COLLECT AND WRITE GLOBAL TALLIES @@ -505,9 +505,8 @@ contains 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 + n_bins = 3 * N_GLOBAL_TALLIES + global_temp(:,:) = global_tallies(:,:) if (master) then ! The MPI_IN_PLACE specifier allows the master to copy values into a @@ -519,20 +518,11 @@ contains ! Transfer values to value on master if (current_batch == n_max_batches .or. satisfy_triggers) then - global_tallies(:)%sum = global_temp(1,:) - global_tallies(:)%sum_sq = global_temp(2,:) + global_tallies(:,:) = global_temp(:,:) end if - ! Put reduced value in temporary tally result - allocate(tallyresult_temp(N_GLOBAL_TALLIES, 1)) - tallyresult_temp(:,1)%sum = global_temp(1,:) - tallyresult_temp(:,1)%sum_sq = global_temp(2,:) - ! Write out global tallies sum and sum_sq - call write_dataset(file_id, "global_tallies", tallyresult_temp) - - ! Deallocate temporary tally result - deallocate(tallyresult_temp) + call write_dataset(file_id, "global_tallies", global_temp) else ! Receive buffer not significant at other processors #ifdef MPI @@ -568,15 +558,15 @@ contains tally => tallies(i) ! Determine size of tally results array - m = size(tally%results, 1) - n = size(tally%results, 2) + m = size(tally%results, 2) + n = size(tally%results, 3) n_bins = m*n*2 ! Allocate array for storing sums and sums of squares, but ! contiguously in memory for each allocate(tally_temp(2,m,n)) - tally_temp(1,:,:) = tally%results(:,:)%sum - tally_temp(2,:,:) = tally%results(:,:)%sum_sq + tally_temp(1,:,:) = tally%results(RESULT_SUM,:,:) + tally_temp(2,:,:) = tally%results(RESULT_SUM_SQ,:,:) if (master) then tally_group = open_group(tallies_group, "tally " // & @@ -592,20 +582,20 @@ contains ! At the end of the simulation, store the results back in the ! regular TallyResults array if (current_batch == n_max_batches .or. satisfy_triggers) then - tally%results(:,:)%sum = tally_temp(1,:,:) - tally%results(:,:)%sum_sq = tally_temp(2,:,:) + tally%results(RESULT_SUM,:,:) = tally_temp(1,:,:) + tally%results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:) end if ! Put in temporary tally result - allocate(tallyresult_temp(m,n)) - tallyresult_temp(:,:)%sum = tally_temp(1,:,:) - tallyresult_temp(:,:)%sum_sq = tally_temp(2,:,:) + allocate(dummy_tally % results(3,m,n)) + dummy_tally % results(RESULT_SUM,:,:) = tally_temp(1,:,:) + dummy_tally % results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:) ! Write reduced tally results to file - call write_dataset(tally_group, "results", tally%results) + call dummy_tally % write_results_hdf5(tally_group) ! Deallocate temporary tally result - deallocate(tallyresult_temp) + deallocate(dummy_tally % results) else ! Receive buffer not significant at other processors #ifdef MPI @@ -771,7 +761,7 @@ contains call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.) ! Read global tally data - call read_dataset(file_id, "global_tallies", global_tallies) + call read_dataset(global_tallies, file_id, "global_tallies") ! Check if tally results are present tallies_group = open_group(file_id, "tallies") @@ -787,7 +777,7 @@ contains ! Read sum, sum_sq, and N for each bin tally_group = open_group(tallies_group, "tally " // & trim(to_str(tally % id))) - call read_dataset(tally_group, "results", tally % results) + call tally % read_results_hdf5(tally_group) call read_dataset(tally % n_realizations, tally_group, & "n_realizations") call close_group(tally_group) diff --git a/src/tally.F90 b/src/tally.F90 index 69e4e8ecd..d4968d9d6 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1,5 +1,7 @@ module tally + use, intrinsic :: ISO_C_BINDING + #ifdef MPI use message_passing #endif @@ -18,7 +20,6 @@ module tally use output, only: header use particle_header, only: LocalCoord, Particle use string, only: to_str - use tally_header, only: TallyResult use tally_filter implicit none @@ -1965,8 +1966,8 @@ contains score = score * calc_pn(t % moment_order(i), p % mu) endif !$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score + t % results(RESULT_VALUE, score_index, filter_index) = & + t % results(RESULT_VALUE, score_index, filter_index) + score case(SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN) @@ -1982,10 +1983,9 @@ contains ! multiply score by the angular flux moments and store !$omp critical (score_general_scatt_yn) - t % results(score_index: score_index + num_nm - 1, filter_index) & - % value = t & - % results(score_index: score_index + num_nm - 1, filter_index)& - % value & + t % results(RESULT_VALUE, score_index: score_index + num_nm - 1, & + filter_index) = t % results(RESULT_VALUE, & + score_index: score_index + num_nm - 1, filter_index) & + score * calc_pn(n, p % mu) * calc_rn(n, p % last_uvw) !$omp end critical (score_general_scatt_yn) end do @@ -2011,10 +2011,9 @@ contains ! multiply score by the angular flux moments and store !$omp critical (score_general_flux_tot_yn) - t % results(score_index: score_index + num_nm - 1, filter_index) & - % value = t & - % results(score_index: score_index + num_nm - 1, filter_index)& - % value & + t % results(RESULT_VALUE, score_index: score_index + num_nm - 1, & + filter_index) = t % results(RESULT_VALUE, & + score_index: score_index + num_nm - 1, filter_index) & + score * calc_rn(n, uvw) !$omp end critical (score_general_flux_tot_yn) end do @@ -2031,8 +2030,8 @@ contains ! get the score and tally it !$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value & + t % results(RESULT_VALUE, score_index, filter_index) = & + t % results(RESULT_VALUE, score_index, filter_index) & + score * calc_pn(n, p % mu) end do i = i + t % moment_order(i) @@ -2040,8 +2039,8 @@ contains case default !$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score + t % results(RESULT_VALUE, score_index, filter_index) = & + t % results(RESULT_VALUE, score_index, filter_index) + score end select @@ -2458,8 +2457,8 @@ contains ! Add score to tally !$omp atomic - t % results(i_score, i_filter) % value = & - t % results(i_score, i_filter) % value + score * filter_weight + t % results(RESULT_VALUE, i_score, i_filter) = & + t % results(RESULT_VALUE, i_score, i_filter) + score * filter_weight ! Case for tallying delayed emissions else if (score_bin == SCORE_DELAYED_NU_FISSION .and. g /= 0) then @@ -2498,8 +2497,8 @@ contains ! Add score to tally !$omp atomic - t % results(i_score, i_filter) % value = & - t % results(i_score, i_filter) % value + score * filter_weight + t % results(RESULT_VALUE, i_score, i_filter) = & + t % results(RESULT_VALUE, i_score, i_filter) + score * filter_weight end if end if end do @@ -2536,8 +2535,8 @@ contains filter_weight = product(filter_weights(:size(t % filters))) !$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score * filter_weight + t % results(RESULT_VALUE, score_index, filter_index) = & + t % results(RESULT_VALUE, score_index, filter_index) + score * filter_weight ! reset original delayed group bin matching_bins(t % find_filter(FILTER_DELAYEDGROUP)) = bin_original @@ -2997,8 +2996,8 @@ contains filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt + t % results(RESULT_VALUE, 1, filter_index) = & + t % results(RESULT_VALUE, 1, filter_index) + p % wgt end if ! Inward current on d1 min surface @@ -3033,8 +3032,8 @@ contains filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt + t % results(RESULT_VALUE, 1, filter_index) = & + t % results(RESULT_VALUE, 1, filter_index) + p % wgt ijk0(d1) = ijk0(d1) - 1 end if @@ -3053,8 +3052,8 @@ contains filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt + t % results(RESULT_VALUE, 1, filter_index) = & + t % results(RESULT_VALUE, 1, filter_index) + p % wgt end if ! Inward current on d1 max surface @@ -3089,8 +3088,8 @@ contains filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt + t % results(RESULT_VALUE, 1, filter_index) = & + t % results(RESULT_VALUE, 1, filter_index) + p % wgt ijk0(d1) = ijk0(d1) + 1 end if @@ -3116,9 +3115,10 @@ 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 + real(C_DOUBLE) :: k_col ! Copy of batch collision estimate of keff + real(C_DOUBLE) :: k_abs ! Copy of batch absorption estimate of keff + real(C_DOUBLE) :: k_tra ! Copy of batch tracklength estimate of keff + real(C_DOUBLE) :: val #ifdef MPI ! Combine tally results onto master process @@ -3142,9 +3142,9 @@ contains if (run_mode == MODE_EIGENVALUE) then 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 = global_tallies(RESULT_VALUE, K_COLLISION) / total_weight + k_abs = global_tallies(RESULT_VALUE, K_ABSORPTION) / total_weight + k_tra = global_tallies(RESULT_VALUE, K_TRACKLENGTH) / 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 @@ -3152,7 +3152,14 @@ contains end if ! Accumulate results for global tallies - call accumulate_result(global_tallies) + do i = 1, size(global_tallies, 2) + val = global_tallies(RESULT_VALUE, i)/total_weight + global_tallies(RESULT_VALUE, i) = ZERO + + global_tallies(RESULT_SUM, i) = global_tallies(RESULT_SUM, i) + val + global_tallies(RESULT_SUM_SQ, i) = & + global_tallies(RESULT_SUM_SQ, i) + val*val + end do end if end subroutine synchronize_tallies @@ -3182,7 +3189,7 @@ contains allocate(tally_temp(m,n)) - tally_temp = t % results(:,:) % value + tally_temp = t % results(RESULT_VALUE,:,:) if (master) then ! The MPI_IN_PLACE specifier allows the master to copy values into @@ -3191,35 +3198,35 @@ contains MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Transfer values to value on master - t % results(:,:) % value = tally_temp + t % results(RESULT_VALUE,:,:) = tally_temp 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) ! Reset value on other processors - t % results(:,:) % value = 0 + t % results(RESULT_VALUE,:,:) = ZERO end if deallocate(tally_temp) end do ! Copy global tallies into array to be reduced - global_temp = global_tallies(:) % value + global_temp = global_tallies(RESULT_VALUE, :) if (master) then call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, & MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Transfer values back to global_tallies on master - global_tallies(:) % value = global_temp + global_tallies(RESULT_VALUE, :) = global_temp else ! Receive buffer not significant at other processors call MPI_REDUCE(global_temp, dummy, N_GLOBAL_TALLIES, & MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Reset value on other processors - global_tallies(:) % value = ZERO + global_tallies(RESULT_VALUE, :) = ZERO end if ! We also need to determine the total starting weight of particles from the @@ -3244,6 +3251,9 @@ contains type(TallyObject), intent(inout) :: t + integer :: i, j + real(C_DOUBLE) :: val + ! Increment number of realizations if (reduce_tallies) then t % n_realizations = t % n_realizations + 1 @@ -3251,92 +3261,59 @@ contains t % n_realizations = t % n_realizations + n_procs end if - ! Accumulate each TallyResult - call accumulate_result(t % results) + ! Accumulate each result + do j = 1, size(t % results, 3) + do i = 1, size(t % results, 2) + val = t % results(RESULT_VALUE, i, j)/total_weight + t % results(RESULT_VALUE, i, j) = ZERO + + t % results(RESULT_SUM, i, j) = & + t % results(RESULT_SUM, i, j) + val + t % results(RESULT_SUM_SQ, i, j) = & + t % results(RESULT_SUM_SQ, i, j) + val*val + end do + end do end subroutine accumulate_tally !=============================================================================== ! TALLY_STATISTICS computes the mean and standard deviation of the mean of each -! tally and stores them in the val and val_sq attributes of the TallyResults -! respectively +! tally and stores them in the RESULT_SUM and RESULT_SUM_SQ positions !=============================================================================== subroutine tally_statistics() - integer :: i ! index in tallies array - type(TallyObject), pointer :: t - - ! Calculate statistics for user-defined tallies - do i = 1, n_tallies - t => tallies(i) - - call statistics_result(t % results, t % n_realizations) - end do - - ! Calculate statistics for global tallies - call statistics_result(global_tallies, n_realizations) - - end subroutine tally_statistics - -!=============================================================================== -! ACCUMULATE_RESULT accumulates results from many histories (or many generations) -! into a single realization of a random variable. -!=============================================================================== - - elemental subroutine accumulate_result(this) - - type(TallyResult), intent(inout) :: this - - real(8) :: val - - ! Add the sum and square of the sum of contributions from a tally result to - ! the variables sum and sum_sq. This will later allow us to calculate a - ! variance on the tallies. - - val = this % value/total_weight - this % sum = this % sum + val - this % sum_sq = this % sum_sq + val*val - - ! Reset the single batch estimate - this % value = ZERO - - end subroutine accumulate_result - -!=============================================================================== -! STATISTICS_RESULT determines the sample mean and the standard deviation of the -! mean for a TallyResult. -!=============================================================================== - - elemental subroutine statistics_result(this, n) - - type(TallyResult), intent(inout) :: this - integer, intent(in) :: n + integer :: j, k ! score/filter indices + integer :: n ! number of realizations ! Calculate sample mean and standard deviation of the mean -- note that we ! have used Bessel's correction so that the estimator of the variance of the ! sample mean is unbiased. - this % sum = this % sum/n - this % sum_sq = sqrt((this % sum_sq/n - this % sum * & - this % sum) / (n - 1)) + do i = 1, n_tallies + n = tallies(i) % n_realizations - end subroutine statistics_result + associate (r => tallies(i) % results) + do k = 1, size(r, 3) + do j = 1, size(r, 2) + r(RESULT_SUM, j, k) = r(RESULT_SUM, j, k) / n + r(RESULT_SUM_SQ, j, k) = sqrt((r(RESULT_SUM_SQ, j, k)/n - & + r(RESULT_SUM, j, k) * r(RESULT_SUM, j, k))/(n - 1)) + end do + end do + end associate + end do -!=============================================================================== -! RESET_RESULT zeroes out the value and accumulated sum and sum-squared for a -! single TallyResult. -!=============================================================================== - - elemental subroutine reset_result(this) - - type(TallyResult), intent(inout) :: this - - this % value = ZERO - this % sum = ZERO - this % sum_sq = ZERO - - end subroutine reset_result + ! Calculate statistics for global tallies + n = n_realizations + associate (r => global_tallies) + do j = 1, size(r, 2) + r(RESULT_SUM, j) = r(RESULT_SUM, j) / n + r(RESULT_SUM_SQ, j) = sqrt((r(RESULT_SUM_SQ, j)/n - & + r(RESULT_SUM, j) * r(RESULT_SUM, j))/(n - 1)) + end do + end associate + end subroutine tally_statistics !=============================================================================== ! SETUP_ACTIVE_USERTALLIES diff --git a/src/tally_header.F90 b/src/tally_header.F90 index fe385458b..eeda03e76 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -1,23 +1,15 @@ module tally_header + use, intrinsic :: ISO_C_BINDING + + use hdf5 + use constants, only: NONE, N_FILTER_TYPES use tally_filter_header, only: TallyFilterContainer use trigger_header, only: TriggerObject - use, intrinsic :: ISO_C_BINDING - implicit none -!=============================================================================== -! TALLYRESULT provides accumulation of results in a particular tally bin -!=============================================================================== - - type, bind(C) :: TallyResult - real(C_DOUBLE) :: value = 0. - real(C_DOUBLE) :: sum = 0. - real(C_DOUBLE) :: sum_sq = 0. - end type TallyResult - !=============================================================================== ! TALLYOBJECT describes a user-specified tally. The region of phase space to ! tally in is given by the TallyFilters and the results are stored in a @@ -68,7 +60,7 @@ module tally_header integer :: total_filter_bins integer :: total_score_bins - type(TallyResult), allocatable :: results(:,:) + real(C_DOUBLE), allocatable :: results(:,:,:) ! reset property - allows a tally to be reset after every batch logical :: reset = .false. @@ -79,6 +71,79 @@ module tally_header ! Tally precision triggers integer :: n_triggers = 0 ! # of triggers type(TriggerObject), allocatable :: triggers(:) ! Array of triggers + + contains + procedure :: write_results_hdf5 + procedure :: read_results_hdf5 end type TallyObject +contains + + subroutine write_results_hdf5(this, group_id) + class(TallyObject), intent(in) :: this + integer(HID_T), intent(in) :: group_id + + integer :: hdf5_err + integer(HID_T) :: dset, dspace + integer(HID_T) :: memspace + integer(HSIZE_T) :: dims(3) + integer(HSIZE_T) :: dims_slab(3) + integer(HSIZE_T) :: offset(3) = [1,0,0] + + ! Create file dataspace + dims_slab(:) = shape(this % results) + dims_slab(1) = 2 + call h5screate_simple_f(3, dims_slab, dspace, hdf5_err) + + ! Create memory dataspace that contains only SUM and SUM_SQ values + dims(:) = shape(this % results) + call h5screate_simple_f(3, dims, memspace, hdf5_err) + call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offset, dims_slab, & + hdf5_err) + + ! Create and write to dataset + call h5dcreate_f(group_id, "results", H5T_NATIVE_DOUBLE, dspace, dset, & + hdf5_err) + call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, this % results, dims_slab, & + hdf5_err, mem_space_id=memspace) + + ! Close identifiers + call h5dclose_f(dset, hdf5_err) + call h5sclose_f(memspace, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + end subroutine write_results_hdf5 + + subroutine read_results_hdf5(this, group_id) + class(TallyObject), intent(inout) :: this + integer(HID_T), intent(in) :: group_id + + integer :: hdf5_err + integer(HID_T) :: dset, dspace + integer(HID_T) :: memspace + integer(HSIZE_T) :: dims(3) + integer(HSIZE_T) :: dims_slab(3) + integer(HSIZE_T) :: offset(3) = [1,0,0] + + ! Create file dataspace + dims_slab(:) = shape(this % results) + dims_slab(1) = 2 + call h5screate_simple_f(3, dims_slab, dspace, hdf5_err) + + ! Create memory dataspace that contains only SUM and SUM_SQ values + dims(:) = shape(this % results) + call h5screate_simple_f(3, dims, memspace, hdf5_err) + call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offset, dims_slab, & + hdf5_err) + + ! Create and write to dataset + call h5dopen_f(group_id, "results", dset, hdf5_err) + call h5dread_f(dset, H5T_NATIVE_DOUBLE, this % results, dims_slab, & + hdf5_err, mem_space_id=memspace) + + ! Close identifiers + call h5dclose_f(dset, hdf5_err) + call h5sclose_f(memspace, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + end subroutine read_results_hdf5 + end module tally_header diff --git a/src/tally_initialize.F90 b/src/tally_initialize.F90 index 1ae403bd6..0e7a50a9c 100644 --- a/src/tally_initialize.F90 +++ b/src/tally_initialize.F90 @@ -20,7 +20,8 @@ contains subroutine configure_tallies() ! Allocate global tallies - allocate(global_tallies(N_GLOBAL_TALLIES)) + allocate(global_tallies(3, N_GLOBAL_TALLIES)) + global_tallies(:,:) = ZERO call setup_tally_arrays() @@ -62,7 +63,8 @@ contains t % total_score_bins = t % n_score_bins * t % n_nuclide_bins ! Allocate results array - allocate(t % results(t % total_score_bins, t % total_filter_bins)) + allocate(t % results(3, t % total_score_bins, t % total_filter_bins)) + t % results(:,:,:) = ZERO end do TALLY_LOOP diff --git a/src/trigger.F90 b/src/trigger.F90 index 39786809b..f41cc46cf 100644 --- a/src/trigger.F90 +++ b/src/trigger.F90 @@ -432,17 +432,19 @@ contains real(8), intent(inout) :: rel_err ! tally relative error integer, intent(in) :: score_index ! tally results score index integer, intent(in) :: filter_index ! tally results filter index - integer :: n ! number of realizations - real(8) :: mean ! tally mean - type(TallyResult) :: tally_result ! pointer to TallyResult - type(TallyObject), pointer :: t ! tally pointer + type(TallyObject), intent(in) :: t ! tally + + integer :: n ! number of realizations + real(8) :: mean ! tally mean + real(8) :: tally_sum, tally_sum_sq ! results for a single tally bin n = t % n_realizations - tally_result = t % results(score_index, filter_index) + tally_sum = t % results(RESULT_SUM, score_index, filter_index) + tally_sum_sq = t % results(RESULT_SUM_SQ, score_index, filter_index) ! Compute the tally mean and standard deviation - mean = tally_result % sum / n - std_dev = sqrt((tally_result % sum_sq / n - mean * mean) / (n - 1)) + mean = tally_sum / n + std_dev = sqrt((tally_sum_sq / n - mean * mean) / (n - 1)) ! Compute the relative error if the mean is non-zero if (mean == ZERO) then