diff --git a/src/finalize.F90 b/src/finalize.F90 index b795cdb4f2..86100d1950 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -9,7 +9,8 @@ module finalize use message_passing #endif - use hdf5_interface, only: h5tclose_f, h5close_f, hdf5_err + use hdf5_interface, only: hdf5_bank_t, hdf5_tallyresult_t + use hdf5, only: h5tclose_f, h5close_f implicit none @@ -22,8 +23,10 @@ contains subroutine finalize_run() + integer :: hdf5_err + ! Start finalization timer - call time_finalize % start() + call time_finalize%start() if (run_mode /= MODE_PLOTTING .and. run_mode /= MODE_PARTICLE) then ! Calculate statistics for tallies and write to tallies.out @@ -37,8 +40,8 @@ contains end if ! Stop timers and show timing statistics - call time_finalize % stop() - call time_total % stop() + call time_finalize%stop() + call time_total%stop() if (master .and. (run_mode /= MODE_PLOTTING .and. & run_mode /= MODE_PARTICLE)) then call print_runtime() diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 041e948ff0..1684c352b7 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -64,10 +64,9 @@ module hdf5_interface public :: file_create public :: file_open public :: file_close + public :: create_group public :: open_group public :: close_group - public :: write_source_bank - public :: read_source_bank public :: write_attribute_string contains @@ -243,7 +242,9 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr @@ -287,9 +288,10 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr ! Set up collective vs. independent I/O @@ -344,7 +346,9 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr @@ -404,9 +408,10 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr ! Set up collective vs. independent I/O @@ -461,7 +466,9 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr @@ -521,9 +528,10 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr ! Set up collective vs. independent I/O @@ -578,7 +586,9 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr @@ -638,9 +648,10 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr ! Set up collective vs. independent I/O @@ -695,7 +706,9 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr @@ -755,9 +768,10 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr ! Set up collective vs. independent I/O @@ -795,7 +809,9 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr @@ -839,9 +855,10 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr ! Set up collective vs. independent I/O @@ -896,7 +913,9 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr @@ -956,9 +975,10 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr ! Set up collective vs. independent I/O @@ -1013,7 +1033,9 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr @@ -1073,9 +1095,10 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr ! Set up collective vs. independent I/O @@ -1130,7 +1153,9 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr @@ -1190,9 +1215,10 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr ! Set up collective vs. independent I/O @@ -1247,7 +1273,9 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr @@ -1307,9 +1335,10 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr ! Set up collective vs. independent I/O @@ -1347,7 +1376,9 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr @@ -1391,9 +1422,10 @@ contains integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle type(c_ptr) :: f_ptr ! Set up collective vs. independent I/O @@ -1432,12 +1464,13 @@ contains integer :: n integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle integer(HSIZE_T) :: dims1(1) integer(HSIZE_T) :: dims2(2) - type(c_ptr) :: f_ptr character(len=len_trim(buffer)), dimension(1) :: str_tmp ! Set up collective vs. independent I/O @@ -1492,12 +1525,13 @@ contains integer :: n integer :: hdf5_err integer :: data_xfer_mode +#ifdef PHDF5 integer(HID_T) :: plist ! property list +#endif integer(HID_T) :: dset ! data set handle integer(HID_T) :: dspace ! data or file space handle integer(HSIZE_T) :: dims1(1) integer(HSIZE_T) :: dims2(2) - type(c_ptr) :: f_ptr character(len=len_trim(buffer)), dimension(1) :: str_tmp ! Set up collective vs. independent I/O @@ -1673,164 +1707,6 @@ contains call h5dclose_f(dset, hdf5_err) end subroutine read_tally_result_2D_explicit -!=============================================================================== -! WRITE_SOURCE_BANK writes OpenMC source_bank data -!=============================================================================== - - subroutine write_source_bank(group_id) - use bank_header, only: Bank - use global, only: n_particles, work, source_bank - - integer(HID_T), intent(in) :: group_id - - integer :: hdf5_err - integer :: data_xfer_mode - integer(HID_T) :: plist ! property list - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - integer(HID_T) :: memspace ! memory space handle - integer(HSIZE_T) :: dims(1) - type(c_ptr) :: f_ptr -#ifdef PHDF5 - integer(HSIZE_T) :: offset(1) ! source data offset -#endif - -#ifdef PHDF5 - ! Set size of total dataspace for all procs and rank - dims(1) = n_particles - call h5screate_simple_f(1, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, "source_bank", hdf5_bank_t, dspace, dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - - ! Create another data space but for each proc individually - dims(1) = work - call h5screate_simple_f(rank, dims, memspace, hdf5_err) - - ! Get the individual local proc dataspace - call h5dget_space_f(dset, dspace, hdf5_err) - - ! Select hyperslab for this dataspace - offset(1) = work_index(rank) - call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims, hdf5_err) - - ! Set up the property list for parallel writing - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, H5FD_MPIO_COLLECTIVE_F, hdf5_err) - - ! Set up pointer to data - f_ptr = c_loc(source_bank) - - ! Write data to file in parallel - call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err, & - file_space_id=dspace, mem_space_id=memspace, & - xfer_prp=plist) - - ! Close all ids - call h5sclose_f(dspace, hdf5_err) - call h5sclose_f(memspace, hdf5_err) - call h5dclose_f(dset, hdf5_err) - call h5pclose_f(plist, hdf5_err) - -#else - - ! Set size - dims(1) = work - - ! Create dataspace - call h5screate_simple_f(1, dims, dspace, hdf5_err) - - ! Create dataset - call h5dcreate_f(group_id, "source_bank", hdf5_bank_t, & - dspace, dset, hdf5_err) - - ! Set up pointer to data - f_ptr = c_loc(source_bank) - - ! Write dataset to file - call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err) - - ! Close all ids - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - -#endif - - end subroutine write_source_bank - -!=============================================================================== -! READ_SOURCE_BANK reads OpenMC source_bank data -!=============================================================================== - - subroutine read_source_bank(group_id) - use bank_header, only: Bank - use global, only: work, source_bank - - integer(HID_T), intent(in) :: group_id - - integer :: hdf5_err - integer :: data_xfer_mode - integer(HID_T) :: plist ! property list - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data space handle - integer(HID_T) :: memspace ! memory space handle - integer(HSIZE_T) :: dims(1) - type(c_ptr) :: f_ptr -#ifdef PHDF5 - integer(HSIZE_T) :: offset(1) ! offset of data -#endif - -#ifdef PHDF5 - - ! Open the dataset - call h5dopen_f(group_id, "source_bank", dset, hdf5_err) - - ! Create another data space but for each proc individually - dims(1) = work - call h5screate_simple_f(1, dims, memspace, hdf5_err) - - ! Get the individual local proc dataspace - call h5dget_space_f(dset, dspace, hdf5_err) - - ! Select hyperslab for this dataspace - offset(1) = work_index(rank) - call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims, hdf5_err) - - ! Set up the property list for parallel writing - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, H5FD_MPIO_COLLECTIVE_F, hdf5_err) - - ! Set up pointer to data - f_ptr = c_loc(source_bank) - - ! Read data from file in parallel - call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err, & - file_space_id=dspace, mem_space_id=memspace, & - xfer_prp=plist) - - ! Close all ids - call h5sclose_f(dspace, hdf5_err) - call h5sclose_f(memspace, hdf5_err) - call h5dclose_f(dset, hdf5_err) - call h5pclose_f(plist, hdf5_err) - -#else - - ! Open dataset - call h5dopen_f(group_id, "source_bank", dset, hdf5_err) - - ! Set up pointer to data - f_ptr = c_loc(source_bank) - - ! Read dataset from file - call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err) - - ! Close all ids - call h5dclose_f(dset, hdf5_err) - -#endif - - end subroutine read_source_bank - function using_mpio_device(obj_id) result(mpio) integer(HID_T), intent(in) :: obj_id logical :: mpio diff --git a/src/hdf5_summary.F90 b/src/hdf5_summary.F90 index a9dab46f50..e8f566c434 100644 --- a/src/hdf5_summary.F90 +++ b/src/hdf5_summary.F90 @@ -6,16 +6,16 @@ module hdf5_summary use geometry_header, only: Cell, Surface, Universe, Lattice, RectLattice, & &HexLattice use global + use hdf5_interface use material_header, only: Material use mesh_header, only: StructuredMesh - use output_interface use output, only: time_stamp use string, only: to_str use tally_header, only: TallyObject - implicit none + use hdf5 - type(BinaryOutput) :: su + implicit none contains @@ -25,48 +25,48 @@ contains subroutine hdf5_write_summary() - character(MAX_FILE_LEN) :: filename = "summary.h5" + integer(HID_T) :: file_id ! Create a new file using default properties. - call su % file_create(filename) + file_id = file_create("summary.h5") ! Write header information - call hdf5_write_header() + call hdf5_write_header(file_id) ! Write eigenvalue information if (run_mode == MODE_EIGENVALUE) then ! Write number of particles - call su % write_data(n_particles, "n_particles") + call write_dataset(file_id, "n_particles", n_particles) ! Use H5LT interface to write n_batches, n_inactive, and n_active - call su % write_data(n_batches, "n_batches") - call su % write_data(n_inactive, "n_inactive") - call su % write_data(n_active, "n_active") - call su % write_data(gen_per_batch, "gen_per_batch") + call write_dataset(file_id, "n_batches", n_batches) + call write_dataset(file_id, "n_inactive", n_inactive) + call write_dataset(file_id, "n_active", n_active) + call write_dataset(file_id, "gen_per_batch", gen_per_batch) ! Add description of each variable - call su % write_attribute_string("n_particles", & + call write_attribute_string(file_id, "n_particles", & "description", "Number of particles per generation") - call su % write_attribute_string("n_batches", & + call write_attribute_string(file_id, "n_batches", & "description", "Total number of batches") - call su % write_attribute_string("n_inactive", & + call write_attribute_string(file_id, "n_inactive", & "description", "Number of inactive batches") - call su % write_attribute_string("n_active", & + call write_attribute_string(file_id, "n_active", & "description", "Number of active batches") - call su % write_attribute_string("gen_per_batch", & + call write_attribute_string(file_id, "gen_per_batch", & "description", "Number of generations per batch") end if - call hdf5_write_geometry() - call hdf5_write_materials() - call hdf5_write_nuclides() + call hdf5_write_geometry(file_id) + call hdf5_write_materials(file_id) + call hdf5_write_nuclides(file_id) if (n_tallies > 0) then - call hdf5_write_tallies() + call hdf5_write_tallies(file_id) end if ! Terminate access to the file. - call su % file_close() + call file_close(file_id) end subroutine hdf5_write_summary @@ -74,19 +74,20 @@ contains ! HDF5_WRITE_HEADER !=============================================================================== - subroutine hdf5_write_header() + subroutine hdf5_write_header(file_id) + integer(HID_T), intent(in) :: file_id ! Write version information - call su % write_data(VERSION_MAJOR, "version_major") - call su % write_data(VERSION_MINOR, "version_minor") - call su % write_data(VERSION_RELEASE, "version_release") + call write_dataset(file_id, "version_major", VERSION_MAJOR) + call write_dataset(file_id, "version_minor", VERSION_MINOR) + call write_dataset(file_id, "version_release", VERSION_RELEASE) ! Write current date and time - call su % write_data(time_stamp(), "date_and_time") + call write_dataset(file_id, "date_and_time", time_stamp()) ! Write MPI information - call su % write_data(n_procs, "n_procs") - call su % write_attribute_string("n_procs", "description", & + call write_dataset(file_id, "n_procs", n_procs) + call write_attribute_string(file_id, "n_procs", "description", & "Number of MPI processes") end subroutine hdf5_write_header @@ -95,782 +96,659 @@ contains ! HDF5_WRITE_GEOMETRY !=============================================================================== - subroutine hdf5_write_geometry() + subroutine hdf5_write_geometry(file_id) + integer(HID_T), intent(in) :: file_id integer :: i, j, k, m integer, allocatable :: lattice_universes(:,:,:) - type(Cell), pointer :: c => null() - type(Surface), pointer :: s => null() - type(Universe), pointer :: u => null() - class(Lattice), pointer :: lat => null() + integer(HID_T) :: geom_group + integer(HID_T) :: cells_group, cell_group + integer(HID_T) :: surfaces_group, surface_group + integer(HID_T) :: universes_group, univ_group + integer(HID_T) :: lattices_group, lattice_group + type(Cell), pointer :: c + type(Surface), pointer :: s + type(Universe), pointer :: u + class(Lattice), pointer :: lat ! Use H5LT interface to write number of geometry objects - call su % write_data(n_cells, "n_cells", group="geometry") - call su % write_data(n_surfaces, "n_surfaces", group="geometry") - call su % write_data(n_universes, "n_universes", group="geometry") - call su % write_data(n_lattices, "n_lattices", group="geometry") + geom_group = create_group(file_id, "geometry") + call write_dataset(geom_group, "n_cells", n_cells) + call write_dataset(geom_group, "n_surfaces", n_surfaces) + call write_dataset(geom_group, "n_universes", n_universes) + call write_dataset(geom_group, "n_lattices", n_lattices) ! ========================================================================== ! WRITE INFORMATION ON CELLS ! Create a cell group (nothing directly written in this group) then close - call su % open_group("geometry/cells") - call su % close_group() + cells_group = create_group(geom_group, "cells") ! Write information on each cell CELL_LOOP: do i = 1, n_cells c => cells(i) + cell_group = create_group(cells_group, "cell " // trim(to_str(c%id))) ! Write internal OpenMC index for this cell - call su % write_data(i, "index", & - group="geometry/cells/cell " // trim(to_str(c % id))) + call write_dataset(cell_group, "index", i) ! Write name for this cell - call su % write_data(c % name, "name", & - group="geometry/cells/cell " // trim(to_str(c % id))) + call write_dataset(cell_group, "name", c%name) ! Write universe for this cell - call su % write_data(universes(c % universe) % id, "universe", & - group="geometry/cells/cell " // trim(to_str(c % id))) + call write_dataset(cell_group, "universe", universes(c%universe)%id) ! Write information on what fills this cell - select case (c % type) + select case (c%type) case (CELL_NORMAL) - call su % write_data("normal", "fill_type", & - group="geometry/cells/cell " // trim(to_str(c % id))) - if (c % material == MATERIAL_VOID) then - call su % write_data(-1, "material", & - group="geometry/cells/cell " // trim(to_str(c % id))) + call write_dataset(cell_group, "fill_type", "normal") + if (c%material == MATERIAL_VOID) then + call write_dataset(cell_group, "material", -1) else - call su % write_data(materials(c % material) % id, "material", & - group="geometry/cells/cell " // trim(to_str(c % id))) + call write_dataset(cell_group, "material", materials(c%material)%id) end if case (CELL_FILL) - call su % write_data("universe", "fill_type", & - group="geometry/cells/cell " // trim(to_str(c % id))) - call su % write_data(universes(c % fill) % id, "fill", & - group="geometry/cells/cell " // trim(to_str(c % id))) - - call su % write_data(size(c % offset), "maps", & - group="geometry/cells/cell " // trim(to_str(c % id))) - if (size(c % offset) > 0) then - call su % write_data(c % offset, "offset", & - length=size(c % offset), & - group="geometry/cells/cell " // trim(to_str(c % id))) + call write_dataset(cell_group, "fill_type", "universe") + call write_dataset(cell_group, "fill", universes(c%fill)%id) + call write_dataset(cell_group, "maps", size(c%offset)) + if (size(c%offset) > 0) then + call write_dataset(cell_group, "offset", c%offset) end if - if (allocated(c % translation)) then - call su % write_data(1, "translated", & - group="geometry/cells/cell " // trim(to_str(c % id))) - call su % write_data(c % translation, "translation", length=3, & - group="geometry/cells/cell " // trim(to_str(c % id))) + if (allocated(c%translation)) then + call write_dataset(cell_group, "translated", 1) + call write_dataset(cell_group, "translation", c%translation) else - call su % write_data(0, "translated", & - group="geometry/cells/cell " // trim(to_str(c % id))) + call write_dataset(cell_group, "translated", 0) end if - if (allocated(c % rotation)) then - call su % write_data(1, "rotated", & - group="geometry/cells/cell " // trim(to_str(c % id))) - call su % write_data(c % rotation, "rotation", length=3, & - group="geometry/cells/cell " // trim(to_str(c % id))) + if (allocated(c%rotation)) then + call write_dataset(cell_group, "rotated", 1) + call write_dataset(cell_group, "rotation", c%rotation) else - call su % write_data(0, "rotated", & - group="geometry/cells/cell " // trim(to_str(c % id))) + call write_dataset(cell_group, "rotated", 0) end if case (CELL_LATTICE) - call su % write_data("lattice", "fill_type", & - group="geometry/cells/cell " // trim(to_str(c % id))) - call su % write_data(lattices(c % fill) % obj % id, "lattice", & - group="geometry/cells/cell " // trim(to_str(c % id))) + call write_dataset(cell_group, "fill_type", "lattice") + call write_dataset(cell_group, "lattice", lattices(c%fill)%obj%id) end select ! Write list of bounding surfaces - if (c % n_surfaces > 0) then - call su % write_data(c % surfaces, "surfaces", length= c % n_surfaces, & - group="geometry/cells/cell " // trim(to_str(c % id))) + if (c%n_surfaces > 0) then + call write_dataset(cell_group, "surfaces", c%surfaces) end if + call close_group(cell_group) end do CELL_LOOP + call close_group(cells_group) + ! ========================================================================== ! WRITE INFORMATION ON SURFACES - ! Create surfaces group (nothing directly written here) then close - call su % open_group("geometry/surfaces") - call su % close_group() + ! Create surfaces group + surfaces_group = create_group(file_id, "surfaces") ! Write information on each surface SURFACE_LOOP: do i = 1, n_surfaces s => surfaces(i) + surface_group = create_group(surfaces_group, "surface " // & + trim(to_str(s%id))) ! Write internal OpenMC index for this surface - call su % write_data(i, "index", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "index", i) ! Write name for this surface - call su % write_data(s % name, "name", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "name", s%name) ! Write surface type - select case (s % type) + select case (s%type) case (SURF_PX) - call su % write_data("X Plane", "type", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "type", "X Plane") case (SURF_PY) - call su % write_data("Y Plane", "type", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "type", "Y Plane") case (SURF_PZ) - call su % write_data("Z Plane", "type", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "type", "Z Plane") case (SURF_PLANE) - call su % write_data("Plane", "type", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "type", "Plane") case (SURF_CYL_X) - call su % write_data("X Cylinder", "type", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "type", "X Cylinder") case (SURF_CYL_Y) - call su % write_data("Y Cylinder", "type", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "type", "Y Cylinder") case (SURF_CYL_Z) - call su % write_data("Z Cylinder", "type", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "type", "Z Cylinder") case (SURF_SPHERE) - call su % write_data("Sphere", "type", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "type", "Sphere") case (SURF_CONE_X) - call su % write_data("X Cone", "type", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "type", "X Cone") case (SURF_CONE_Y) - call su % write_data("Y Cone", "type", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "type", "Y Cone") case (SURF_CONE_Z) - call su % write_data("Z Cone", "type", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "type", "Z Cone") end select ! Write coefficients for surface - call su % write_data(s % coeffs, "coefficients", length=size(s % coeffs), & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "coefficients", s%coeffs) ! Write positive neighbors - if (allocated(s % neighbor_pos)) then - call su % write_data(s % neighbor_pos, "neighbors_positive", & - length=size(s % neighbor_pos), & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + if (allocated(s%neighbor_pos)) then + call write_dataset(surface_group, "neighbors_positive", s%neighbor_pos) end if ! Write negative neighbors - if (allocated(s % neighbor_neg)) then - call su % write_data(s % neighbor_neg, "neighbors_negative", & - length=size(s % neighbor_neg), & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + if (allocated(s%neighbor_neg)) then + call write_dataset(surface_group, "neighbors_negative", s%neighbor_neg) end if ! Write boundary condition - select case (s % bc) + select case (s%bc) case (BC_TRANSMIT) - call su % write_data("transmission", "boundary_condition", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "boundary_condition", "transmission") case (BC_VACUUM) - call su % write_data("vacuum", "boundary_condition", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "boundary_condition", "vacuum") case (BC_REFLECT) - call su % write_data("reflective", "boundary_condition", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "boundary_condition", "reflective") case (BC_PERIODIC) - call su % write_data("periodic", "boundary_condition", & - group="geometry/surfaces/surface " // trim(to_str(s % id))) + call write_dataset(surface_group, "boundary_condition", "periodic") end select + call close_group(surface_group) end do SURFACE_LOOP + call close_group(surfaces_group) + ! ========================================================================== ! WRITE INFORMATION ON UNIVERSES ! Create universes group (nothing directly written here) then close - call su % open_group("geometry/universes") - call su % close_group() + universes_group = create_group(geom_group, "universes") ! Write information on each universe UNIVERSE_LOOP: do i = 1, n_universes u => universes(i) + univ_group = create_group(universes_group, "universe " // & + trim(to_str(u%id))) ! Write internal OpenMC index for this universe - call su % write_data(i, "index", & - group="geometry/universes/universe " // trim(to_str(u % id))) + call write_dataset(univ_group, "index", i) ! Write list of cells in this universe - if (u % n_cells > 0) then - call su % write_data(u % cells, "cells", length=u % n_cells, & - group="geometry/universes/universe " // trim(to_str(u % id))) - end if + if (u%n_cells > 0) call write_dataset(univ_group, "cells", u%cells) + call close_group(univ_group) end do UNIVERSE_LOOP + call close_group(universes_group) + ! ========================================================================== ! WRITE INFORMATION ON LATTICES ! Create lattices group (nothing directly written here) then close - call su % open_group("geometry/lattices") - call su % close_group() + lattices_group = create_group(geom_group, "lattices") ! Write information on each lattice LATTICE_LOOP: do i = 1, n_lattices - lat => lattices(i) % obj + lat => lattices(i)%obj + lattice_group = create_group(lattices_group, "lattice " // trim(to_str(lat%id))) ! Write internal OpenMC index for this lattice - call su % write_data(i, "index", & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) + call write_dataset(lattice_group, "index", i) ! Write name for this lattice - call su % write_data(lat % name, "name", & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) + call write_dataset(lattice_group, "name", lat%name) ! Write lattice type select type (lat) type is (RectLattice) ! Write lattice type. - call su % write_data("rectangular", "type", & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) + call write_dataset(lattice_group, "type", "rectangular") ! Write lattice dimensions, lower left corner, and pitch - call su % write_data(lat % n_cells, "dimension", length=3, & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) + call write_dataset(lattice_group, "dimension", lat%n_cells) + call write_dataset(lattice_group, "lower_left", lat%lower_left) + call write_dataset(lattice_group, "pitch", lat%pitch) - if (lat % is_3d) then - call su % write_data(lat % lower_left, "lower_left", length=3, & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - else - call su % write_data(lat % lower_left, "lower_left", length=2, & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - end if + call write_dataset(lattice_group, "outer", lat%outer) + call write_dataset(lattice_group, "offset_size", size(lat%offset)) + call write_dataset(lattice_group, "maps", size(lat%offset,1)) - if (lat % is_3d) then - call su % write_data(lat % pitch, "pitch", length=3, & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - else - call su % write_data(lat % pitch, "pitch", length=2, & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - end if - - call su % write_data(lat % outer, "outer", & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - call su % write_data(size(lat % offset), "offset_size", & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - call su % write_data(size(lat % offset,1), "maps", & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - - if (size(lat % offset) > 0) then - call su % write_data(lat % offset, "offsets", & - length=shape(lat % offset), & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) + if (size(lat%offset) > 0) then + call write_dataset(lattice_group, "offsets", lat%offset) end if ! Write lattice universes. - allocate(lattice_universes(lat % n_cells(1), lat % n_cells(2), & - &lat % n_cells(3))) - do j = 1, lat % n_cells(1) - do k = 1, lat % n_cells(2) - do m = 1, lat % n_cells(3) - lattice_universes(j,k,m) = universes(lat % universes(j,k,m)) % id + allocate(lattice_universes(lat%n_cells(1), lat%n_cells(2), & + &lat%n_cells(3))) + do j = 1, lat%n_cells(1) + do k = 1, lat%n_cells(2) + do m = 1, lat%n_cells(3) + lattice_universes(j,k,m) = universes(lat%universes(j,k,m))%id end do end do end do - call su % write_data(lattice_universes, "universes", & - length=lat % n_cells, & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) + call write_dataset(lattice_group, "universes", lattice_universes) deallocate(lattice_universes) type is (HexLattice) ! Write lattice type. - call su % write_data("hexagonal", "type", & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) + call write_dataset(lattice_group, "type", "hexagonal") ! Write number of lattice cells. - call su % write_data(lat % n_rings, "n_rings", & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - call su % write_data(lat % n_axial, "n_axial", & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) + call write_dataset(lattice_group, "n_rings", lat%n_rings) + call write_dataset(lattice_group, "n_axial", lat%n_axial) ! Write lattice center, pitch and outer universe. - if (lat % is_3d) then - call su % write_data(lat % center, "center", length=3, & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - else - call su % write_data(lat % center, "center", length=2, & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - end if + call write_dataset(lattice_group, "center", lat%center) + call write_dataset(lattice_group, "pitch", lat%pitch) - if (lat % is_3d) then - call su % write_data(lat % pitch, "pitch", length=2, & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - else - call su % write_data(lat % pitch, "pitch", length=1, & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - end if + call write_dataset(lattice_group, "outer", lat%outer) + call write_dataset(lattice_group, "offset_size", size(lat%offset)) + call write_dataset(lattice_group, "maps", size(lat%offset,1)) - call su % write_data(lat % outer, "outer", & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - call su % write_data(size(lat % offset), "offset_size", & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - call su % write_data(size(lat % offset,1), "maps", & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) - - if (size(lat % offset) > 0) then - call su % write_data(lat % offset, "offsets", & - length=shape(lat % offset), & - group="geometry/lattices/lattice " // trim(to_str(lat % id))) + if (size(lat%offset) > 0) then + call write_dataset(lattice_group, "offsets", lat%offset) end if ! Write lattice universes. - allocate(lattice_universes(2*lat % n_rings - 1, 2*lat % n_rings - 1, & - &lat % n_axial)) - do m = 1, lat % n_axial - do k = 1, 2*lat % n_rings - 1 - do j = 1, 2*lat % n_rings - 1 - if (j + k < lat % n_rings + 1) then + allocate(lattice_universes(2*lat%n_rings - 1, 2*lat%n_rings - 1, & + &lat%n_axial)) + do m = 1, lat%n_axial + do k = 1, 2*lat%n_rings - 1 + do j = 1, 2*lat%n_rings - 1 + if (j + k < lat%n_rings + 1) then ! This array position is never used; put a -1 to indicate this lattice_universes(j,k,m) = -1 cycle - else if (j + k > 3*lat % n_rings - 1) then + else if (j + k > 3*lat%n_rings - 1) then ! This array position is never used; put a -1 to indicate this lattice_universes(j,k,m) = -1 cycle end if - lattice_universes(j,k,m) = universes(lat % universes(j,k,m)) % id + lattice_universes(j,k,m) = universes(lat%universes(j,k,m))%id end do end do end do - call su % write_data(lattice_universes, "universes", & - &length=(/2*lat % n_rings-1, 2*lat % n_rings-1, lat % n_axial/), & - &group="geometry/lattices/lattice " // trim(to_str(lat % id))) + call write_dataset(lattice_group, "universes", lattice_universes) deallocate(lattice_universes) end select + + call close_group(lattice_group) end do LATTICE_LOOP + call close_group(lattices_group) + call close_group(geom_group) + end subroutine hdf5_write_geometry !=============================================================================== ! HDF5_WRITE_MATERIALS !=============================================================================== - subroutine hdf5_write_materials() + subroutine hdf5_write_materials(file_id) + integer(HID_T), intent(in) :: file_id integer :: i integer :: j integer, allocatable :: zaids(:) - type(Material), pointer :: m => null() + integer(HID_T) :: materials_group + integer(HID_T) :: material_group + integer(HID_T) :: sab_group + type(Material), pointer :: m - ! Use H5LT interface to write number of materials - call su % write_data(n_materials, "n_materials", group="materials") + materials_group = create_group(file_id, "materials") + + ! write number of materials + call write_dataset(file_id, "n_materials", n_materials) ! Write information on each material do i = 1, n_materials m => materials(i) + material_group = create_group(materials_group, "material " // & + trim(to_str(m%id))) ! Write internal OpenMC index for this material - call su % write_data(i, "index", & - group="materials/material " // trim(to_str(m % id))) + call write_dataset(material_group, "index", i) ! Write name for this material - call su % write_data(m % name, "name", & - group="materials/material " // trim(to_str(m % id))) + call write_dataset(material_group, "name", m%name) ! Write atom density with units - call su % write_data(m % density, "atom_density", & - group="materials/material " // trim(to_str(m % id))) - call su % write_attribute_string("atom_density", "units", "atom/b-cm", & - group="materials/material " // trim(to_str(m % id))) + call write_dataset(material_group, "atom_density", m%density) + call write_attribute_string(material_group, "atom_density", "units", & + "atom/b-cm") ! Copy ZAID for each nuclide to temporary array - allocate(zaids(m % n_nuclides)) - do j = 1, m % n_nuclides - zaids(j) = nuclides(m % nuclide(j)) % zaid + allocate(zaids(m%n_nuclides)) + do j = 1, m%n_nuclides + zaids(j) = nuclides(m%nuclide(j))%zaid end do ! Write temporary array to 'nuclides' - call su % write_data(zaids, "nuclides", length=m % n_nuclides, & - group="materials/material " // trim(to_str(m % id))) + call write_dataset(material_group, "nuclides", zaids) ! Deallocate temporary array deallocate(zaids) ! Write atom densities - call su % write_data(m % atom_density, "nuclide_densities", & - length=m % n_nuclides, & - group="materials/material " // trim(to_str(m % id))) + call write_dataset(material_group, "nuclide_densities", m%atom_density) ! Write S(a,b) information if present - call su % write_data(m % n_sab, "n_sab", & - group="materials/material " // trim(to_str(m % id))) + call write_dataset(material_group, "n_sab", m%n_sab) - if (m % n_sab > 0) then - call su % write_data(m % i_sab_nuclides, "i_sab_nuclides", & - length=m % n_sab, & - group="materials/material " // trim(to_str(m % id))) - call su % write_data(m % i_sab_tables, "i_sab_tables", & - length=m % n_sab, & - group="materials/material " // trim(to_str(m % id))) + if (m%n_sab > 0) then + call write_dataset(material_group, "i_sab_nuclides", m%i_sab_nuclides) + call write_dataset(material_group, "i_sab_tables", m%i_sab_tables) - do j = 1, m % n_sab - call su % write_data(m % sab_names(j), to_str(j), & - group="materials/material " // & - trim(to_str(m % id)) // "/sab_tables") + sab_group = create_group(material_group, "sab_tables") + do j = 1, m%n_sab + call write_dataset(sab_group, to_str(j), m%sab_names(j)) end do + call close_group(sab_group) end if + call close_group(material_group) end do + call close_group(materials_group) + end subroutine hdf5_write_materials !=============================================================================== ! HDF5_WRITE_TALLIES !=============================================================================== - subroutine hdf5_write_tallies() + subroutine hdf5_write_tallies(file_id) + integer(HID_T), intent(in) :: file_id integer :: i, j integer, allocatable :: temp_array(:) ! nuclide bin array - type(StructuredMesh), pointer :: m => null() - type(TallyObject), pointer :: t => null() + integer(HID_T) :: tallies_group + integer(HID_T) :: mesh_group + integer(HID_T) :: tally_group + integer(HID_T) :: filter_group + type(StructuredMesh), pointer :: m + type(TallyObject), pointer :: t + + tallies_group = create_group(file_id, "tallies") ! Write total number of meshes - call su % write_data(n_meshes, "n_meshes", group="tallies") + call write_dataset(tallies_group, "n_meshes", n_meshes) ! Write information for meshes MESH_LOOP: do i = 1, n_meshes m => meshes(i) + mesh_group = create_group(tallies_group, "mesh " // trim(to_str(m%id))) ! Write type and number of dimensions - call su % write_data(m % type, "type", & - group="tallies/mesh " // trim(to_str(m % id))) - - call su % write_data(m % n_dimension, "n_dimension", & - group="tallies/mesh " // trim(to_str(m % id))) + call write_dataset(mesh_group, "type", m%type) + call write_dataset(mesh_group, "n_dimension", m%n_dimension) ! Write mesh information - call su % write_data(m % dimension, "dimension", & - length=m % n_dimension, & - group="tallies/mesh " // trim(to_str(m % id))) - call su % write_data(m % lower_left, "lower_left", & - length=m % n_dimension, & - group="tallies/mesh " // trim(to_str(m % id))) - call su % write_data(m % upper_right, "upper_right", & - length=m % n_dimension, & - group="tallies/mesh " // trim(to_str(m % id))) - call su % write_data(m % width, "width", & - length=m % n_dimension, & - group="tallies/mesh " // trim(to_str(m % id))) + call write_dataset(mesh_group, "dimension", m%dimension) + call write_dataset(mesh_group, "lower_left", m%lower_left) + call write_dataset(mesh_group, "upper_right", m%upper_right) + call write_dataset(mesh_group, "width", m%width) + call close_group(mesh_group) end do MESH_LOOP ! Write number of tallies - call su % write_data(n_tallies, "n_tallies", group="tallies") + call write_dataset(tallies_group, "n_tallies", n_tallies) TALLY_METADATA: do i = 1, n_tallies ! Get pointer to tally t => tallies(i) + tally_group = create_group(tallies_group, "tally " // trim(to_str(t%id))) ! Write the name for this tally - call su % write_data(len(t % name), "name_size", & - group="tallies/tally " // trim(to_str(t % id))) - if (len(t % name) > 0) then - call su % write_data(t % name, "name", & - group="tallies/tally " // trim(to_str(t % id))) + call write_dataset(tally_group, "name_size", len(t%name)) + if (len(t%name) > 0) then + call write_dataset(tally_group, "name", t%name) endif ! Write size of each tally - call su % write_data(t % total_score_bins, "total_score_bins", & - group="tallies/tally " // trim(to_str(t % id))) - call su % write_data(t % total_filter_bins, "total_filter_bins", & - group="tallies/tally " // trim(to_str(t % id))) + call write_dataset(tally_group, "total_score_bins", t%total_score_bins) + call write_dataset(tally_group, "total_filter_bins", t%total_filter_bins) ! Write number of filters - call su % write_data(t % n_filters, "n_filters", & - group="tallies/tally " // trim(to_str(t % id))) + call write_dataset(tally_group, "n_filters", t%n_filters) + + FILTER_LOOP: do j = 1, t%n_filters + filter_group = create_group(tally_group, "filter " // trim(to_str(j))) - FILTER_LOOP: do j = 1, t % n_filters ! Write type of filter - call su % write_data(t % filters(j) % type, "type", & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + call write_dataset(filter_group, "type", t%filters(j)%type) ! Write number of bins for this filter - call su % write_data(t % filters(j) % n_bins, "n_bins", & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + call write_dataset(filter_group, "n_bins", t%filters(j)%n_bins) ! Write filter bins - if (t % filters(j) % type == FILTER_ENERGYIN .or. & - t % filters(j) % type == FILTER_ENERGYOUT) then - call su % write_data(t % filters(j) % real_bins, "bins", & - length=size(t % filters(j) % real_bins), & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + if (t%filters(j)%type == FILTER_ENERGYIN .or. & + t%filters(j)%type == FILTER_ENERGYOUT) then + call write_dataset(filter_group, "bins", t%filters(j)%real_bins) else - call su % write_data(t % filters(j) % int_bins, "bins", & - length=size(t % filters(j) % int_bins), & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + call write_dataset(filter_group, "bins", t%filters(j)%int_bins) end if ! Write name of type - select case (t % filters(j) % type) + select case (t%filters(j)%type) case(FILTER_UNIVERSE) - call su % write_data("universe", "type_name", & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + call write_dataset(filter_group, "type_name", "universe") case(FILTER_MATERIAL) - call su % write_data("material", "type_name", & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + call write_dataset(filter_group, "type_name", "material") case(FILTER_CELL) - call su % write_data("cell", "type_name", & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + call write_dataset(filter_group, "type_name", "cell") case(FILTER_CELLBORN) - call su % write_data("cellborn", "type_name", & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + call write_dataset(filter_group, "type_name", "cellborn") case(FILTER_SURFACE) - call su % write_data("surface", "type_name", & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + call write_dataset(filter_group, "type_name", "surface") case(FILTER_MESH) - call su % write_data("mesh", "type_name", & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + call write_dataset(filter_group, "type_name", "mesh") case(FILTER_ENERGYIN) - call su % write_data("energy", "type_name", & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + call write_dataset(filter_group, "type_name", "energy") case(FILTER_ENERGYOUT) - call su % write_data("energyout", "type_name", & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + call write_dataset(filter_group, "type_name", "energyout") end select + call close_group(filter_group) end do FILTER_LOOP ! Write number of nuclide bins - call su % write_data(t % n_nuclide_bins, "n_nuclide_bins", & - group="tallies/tally " // trim(to_str(t % id))) + call write_dataset(tally_group, "n_nuclide_bins", t%n_nuclide_bins) ! Create temporary array for nuclide bins - allocate(temp_array(t % n_nuclide_bins)) - NUCLIDE_LOOP: do j = 1, t % n_nuclide_bins - if (t % nuclide_bins(j) > 0) then - temp_array(j) = nuclides(t % nuclide_bins(j)) % zaid + allocate(temp_array(t%n_nuclide_bins)) + NUCLIDE_LOOP: do j = 1, t%n_nuclide_bins + if (t%nuclide_bins(j) > 0) then + temp_array(j) = nuclides(t%nuclide_bins(j))%zaid else - temp_array(j) = t % nuclide_bins(j) + temp_array(j) = t%nuclide_bins(j) end if end do NUCLIDE_LOOP ! Write and deallocate nuclide bins - call su % write_data(temp_array, "nuclide_bins", length=t % n_nuclide_bins, & - group="tallies/tally " // trim(to_str(t % id))) + call write_dataset(tally_group, "nuclide_bins", temp_array) deallocate(temp_array) ! Write number of score bins - call su % write_data(t % n_score_bins, "n_score_bins", & - group="tallies/tally " // trim(to_str(t % id))) - call su % write_data(t % score_bins, "score_bins", length=t % n_score_bins, & - group="tallies/tally " // trim(to_str(t % id))) + call write_dataset(tally_group, "n_score_bins", t%n_score_bins) + call write_dataset(tally_group, "score_bins", t%score_bins) + call close_group(tally_group) end do TALLY_METADATA + call close_group(tallies_group) + end subroutine hdf5_write_tallies !=============================================================================== ! HDF5_WRITE_NUCLIDES !=============================================================================== - subroutine hdf5_write_nuclides() + subroutine hdf5_write_nuclides(file_id) + integer(HID_T), intent(in) :: file_id integer :: i, j integer :: size_total integer :: size_xs integer :: size_angle integer :: size_energy - type(Nuclide), pointer :: nuc => null() - type(Reaction), pointer :: rxn => null() - type(UrrData), pointer :: urr => null() + integer(HID_T) :: nuclides_group, nuclide_group + integer(HID_T) :: reactions_group, rxn_group + type(Nuclide), pointer :: nuc + type(Reaction), pointer :: rxn + type(UrrData), pointer :: urr - ! Use H5LT interface to write number of nuclides - call su % write_data(n_nuclides_total, "n_nuclides", group="nuclides") + nuclides_group = create_group(file_id, "nuclides") + + ! write number of nuclides + call write_dataset(nuclides_group, "n_nuclides", n_nuclides_total) ! Write information on each nuclide NUCLIDE_LOOP: do i = 1, n_nuclides_total nuc => nuclides(i) + nuclide_group = create_group(nuclides_group, nuc%name) ! Write internal OpenMC index for this nuclide - call su % write_data(i, "index", & - group="nuclides/" // trim(nuc % name)) + call write_dataset(nuclide_group, "index", i) ! Determine size of cross-sections - size_xs = (5 + nuc % n_reaction) * nuc % n_grid * 8 + size_xs = (5 + nuc%n_reaction) * nuc%n_grid * 8 size_total = size_xs ! Write some basic attributes - call su % write_data(nuc % zaid, "zaid", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(xs_listings(nuc % listing) % alias, "alias", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(nuc % awr, "awr", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(nuc % kT, "kT", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(nuc % n_grid, "n_grid", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(nuc % n_reaction, "n_reactions", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(nuc % n_fission, "n_fission", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(size_xs, "size_xs", & - group="nuclides/" // trim(nuc % name)) + call write_dataset(nuclide_group, "zaid", nuc%zaid) + call write_dataset(nuclide_group, "alias", xs_listings(nuc%listing)%alias) + call write_dataset(nuclide_group, "awr", nuc%awr) + call write_dataset(nuclide_group, "kT", nuc%kT) + call write_dataset(nuclide_group, "n_grid", nuc%n_grid) + call write_dataset(nuclide_group, "n_reactions", nuc%n_reaction) + call write_dataset(nuclide_group, "n_fission", nuc%n_fission) + call write_dataset(nuclide_group, "size_xs", size_xs) ! ======================================================================= ! WRITE INFORMATION ON EACH REACTION ! Create overall group for reactions and close it - call su % open_group("nuclides/" // trim(nuc % name) // "/reactions") - call su % close_group() + reactions_group = create_group(nuclide_group, "reactions") - RXN_LOOP: do j = 1, nuc % n_reaction + RXN_LOOP: do j = 1, nuc%n_reaction ! Information on each reaction - rxn => nuc % reactions(j) + rxn => nuc%reactions(j) + rxn_group = create_group(reactions_group, trim(reaction_name(rxn%MT))) ! Determine size of angle distribution - if (rxn % has_angle_dist) then - size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8 + if (rxn%has_angle_dist) then + size_angle = rxn%adist%n_energy * 16 + size(rxn%adist%data) * 8 else size_angle = 0 end if ! Determine size of energy distribution - if (rxn % has_energy_dist) then - size_energy = size(rxn % edist % data) * 8 + if (rxn%has_energy_dist) then + size_energy = size(rxn%edist%data) * 8 else size_energy = 0 end if ! Write information on reaction - call su % write_data(rxn % Q_value, "Q_value", & - group="nuclides/" // trim(nuc % name) // "/reactions/" // & - trim(reaction_name(rxn % MT))) - call su % write_data(rxn % multiplicity, "multiplicity", & - group="nuclides/" // trim(nuc % name) // "/reactions/" // & - trim(reaction_name(rxn % MT))) - call su % write_data(rxn % threshold, "threshold", & - group="nuclides/" // trim(nuc % name) // "/reactions/" // & - trim(reaction_name(rxn % MT))) - call su % write_data(size_angle, "size_angle", & - group="nuclides/" // trim(nuc % name) // "/reactions/" // & - trim(reaction_name(rxn % MT))) - call su % write_data(size_energy, "size_energy", & - group="nuclides/" // trim(nuc % name) // "/reactions/" // & - trim(reaction_name(rxn % MT))) + call write_dataset(rxn_group, "Q_value", rxn%Q_value) + call write_dataset(rxn_group, "multiplicity", rxn%multiplicity) + call write_dataset(rxn_group, "threshold", rxn%threshold) + call write_dataset(rxn_group, "size_angle", size_angle) + call write_dataset(rxn_group, "size_energy", size_energy) ! Accumulate data size size_total = size_total + size_angle + size_energy + + call close_group(rxn_group) end do RXN_LOOP + call close_group(reactions_group) + ! ======================================================================= ! WRITE INFORMATION ON URR PROBABILITY TABLES - if (nuc % urr_present) then - urr => nuc % urr_data - call su % write_data(urr % n_energy, "urr_n_energy", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(urr % n_prob, "urr_n_prob", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(urr % interp, "urr_interp", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(urr % inelastic_flag, "urr_inelastic", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(urr % absorption_flag, "urr_absorption", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(urr % energy(1), "urr_min_E", & - group="nuclides/" // trim(nuc % name)) - call su % write_data(urr % energy(urr % n_energy), "urr_max_E", & - group="nuclides/" // trim(nuc % name)) + if (nuc%urr_present) then + urr => nuc%urr_data + call write_dataset(nuclide_group, "urr_n_energy", urr%n_energy) + call write_dataset(nuclide_group, "urr_n_prob", urr%n_prob) + call write_dataset(nuclide_group, "urr_interp", urr%interp) + call write_dataset(nuclide_group, "urr_inelastic", urr%inelastic_flag) + call write_dataset(nuclide_group, "urr_absorption", urr%absorption_flag) + call write_dataset(nuclide_group, "urr_min_E", urr%energy(1)) + call write_dataset(nuclide_group, "urr_max_E", urr%energy(urr%n_energy)) end if ! Write total memory used - call su % write_data(size_total, "size_total", & - group="nuclides/" // trim(nuc % name)) + call write_dataset(nuclide_group, "size_total", size_total) + call close_group(nuclide_group) end do NUCLIDE_LOOP + call close_group(nuclides_group) + end subroutine hdf5_write_nuclides !=============================================================================== ! HDF5_WRITE_TIMING !=============================================================================== - subroutine hdf5_write_timing() + subroutine hdf5_write_timing(file_id) + integer(HID_T), intent(in) :: file_id - integer(8) :: total_particles - real(8) :: speed + integer(8) :: total_particles + integer(HID_T) :: time_group + real(8) :: speed + + time_group = create_group(file_id, "timing") ! Write timing data - call su % write_data(time_initialize % elapsed, "time_initialize", & - group="timing") - call su % write_data(time_read_xs % elapsed, "time_read_xs", & - group="timing") - call su % write_data(time_transport % elapsed, "time_transport", & - group="timing") - call su % write_data(time_bank % elapsed, "time_bank", & - group="timing") - call su % write_data(time_bank_sample % elapsed, "time_bank_sample", & - group="timing") - call su % write_data(time_bank_sendrecv % elapsed, "time_bank_sendrecv", & - group="timing") - call su % write_data(time_tallies % elapsed, "time_tallies", & - group="timing") - call su % write_data(time_inactive % elapsed, "time_inactive", & - group="timing") - call su % write_data(time_active % elapsed, "time_active", & - group="timing") - call su % write_data(time_finalize % elapsed, "time_finalize", & - group="timing") - call su % write_data(time_total % elapsed, "time_total", & - group="timing") + call write_dataset(time_group, "time_initialize", time_initialize%elapsed) + call write_dataset(time_group, "time_read_xs", time_read_xs%elapsed) + call write_dataset(time_group, "time_transport", time_transport%elapsed) + call write_dataset(time_group, "time_bank", time_bank%elapsed) + call write_dataset(time_group, "time_bank_sample", time_bank_sample%elapsed) + call write_dataset(time_group, "time_bank_sendrecv", time_bank_sendrecv%elapsed) + call write_dataset(time_group, "time_tallies", time_tallies%elapsed) + call write_dataset(time_group, "time_inactive", time_inactive%elapsed) + call write_dataset(time_group, "time_active", time_active%elapsed) + call write_dataset(time_group, "time_finalize", time_finalize%elapsed) + call write_dataset(time_group, "time_total", time_total%elapsed) ! Add descriptions to timing data - call su % write_attribute_string("time_initialize", "description", & - "Total time elapsed for initialization (s)", group="timing") - call su % write_attribute_string("time_read_xs", "description", & - "Time reading cross-section libraries (s)", group="timing") - call su % write_attribute_string("time_transport", "description", & - "Time in transport only (s)", group="timing") - call su % write_attribute_string("time_bank", "description", & - "Total time synchronizing fission bank (s)", group="timing") - call su % write_attribute_string("time_bank_sample", "description", & - "Time between generations sampling source sites (s)", group="timing") - call su % write_attribute_string("time_bank_sendrecv", "description", & - "Time between generations SEND/RECVing source sites (s)", & - group="timing") - call su % write_attribute_string("time_tallies", "description", & - "Time between batches accumulating tallies (s)", group="timing") - call su % write_attribute_string("time_inactive", "description", & - "Total time in inactive batches (s)", group="timing") - call su % write_attribute_string("time_active", "description", & - "Total time in active batches (s)", group="timing") - call su % write_attribute_string("time_finalize", "description", & - "Total time for finalization (s)", group="timing") - call su % write_attribute_string("time_total", "description", & - "Total time elapsed (s)", group="timing") + call write_attribute_string(time_group, "time_initialize", "description", & + "Total time elapsed for initialization (s)") + call write_attribute_string(time_group, "time_read_xs", "description", & + "Time reading cross-section libraries (s)") + call write_attribute_string(time_group, "time_transport", "description", & + "Time in transport only (s)") + call write_attribute_string(time_group, "time_bank", "description", & + "Total time synchronizing fission bank (s)") + call write_attribute_string(time_group, "time_bank_sample", "description", & + "Time between generations sampling source sites (s)") + call write_attribute_string(time_group, "time_bank_sendrecv", "description", & + "Time between generations SEND/RECVing source sites (s)") + call write_attribute_string(time_group, "time_tallies", "description", & + "Time between batches accumulating tallies (s)") + call write_attribute_string(time_group, "time_inactive", "description", & + "Total time in inactive batches (s)") + call write_attribute_string(time_group, "time_active", "description", & + "Total time in active batches (s)") + call write_attribute_string(time_group, "time_finalize", "description", & + "Total time for finalization (s)") + call write_attribute_string(time_group, "time_total", "description", & + "Total time elapsed (s)") ! Write calculation rate total_particles = n_particles * n_batches * gen_per_batch - speed = real(total_particles) / (time_inactive % elapsed + & - time_active % elapsed) - call su % write_data(speed, "neutrons_per_second", group="timing") + speed = real(total_particles) / (time_inactive%elapsed + & + time_active%elapsed) + call write_dataset(time_group, "neutrons_per_second", speed) + call close_group(time_group) end subroutine hdf5_write_timing end module hdf5_summary diff --git a/src/initialize.F90 b/src/initialize.F90 index 22521b24da..29e5f8921f 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -12,12 +12,14 @@ module initialize use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,& &BASE_UNIVERSE use global + use hdf5_interface, only: file_open, read_dataset, file_close, hdf5_bank_t,& + hdf5_tallyresult_t, hdf5_integer8_t + use hdf5_summary, only: hdf5_write_summary use input_xml, only: read_input_xml, read_cross_sections_xml, & cells_in_univ_dict, read_plots_xml use material_header, only: Material use output, only: title, header, print_version, write_message, & print_usage, write_xs_summary, print_plot - use output_interface use random_lcg, only: initialize_prng use state_point, only: load_state_point use string, only: to_str, str_to_int, starts_with, ends_with @@ -32,8 +34,7 @@ module initialize use omp_lib #endif - use hdf5_interface - use hdf5_summary, only: hdf5_write_summary + use hdf5 implicit none @@ -49,8 +50,8 @@ contains subroutine initialize_run() ! Start total and initialization timer - call time_total % start() - call time_initialize % start() + call time_total%start() + call time_initialize%start() #ifdef MPI ! Setup MPI @@ -110,9 +111,9 @@ contains call normalize_ao() ! Read ACE-format cross sections - call time_read_xs % start() + call time_read_xs%start() call read_xs() - call time_read_xs % stop() + call time_read_xs%stop() ! Create linked lists for multiple instances of the same nuclide call same_nuclide_list() @@ -122,9 +123,9 @@ contains case (GRID_NUCLIDE) continue case (GRID_MAT_UNION) - call time_unionize % start() + call time_unionize%start() call unionized_grid() - call time_unionize % stop() + call time_unionize%stop() case (GRID_LOGARITHM) call logarithmic_grid() end select @@ -167,7 +168,7 @@ contains end if ! Stop initialization timer - call time_initialize % stop() + call time_initialize%stop() end subroutine initialize_run @@ -220,10 +221,10 @@ contains ! CREATE MPI_BANK TYPE ! Determine displacements for MPI_BANK type - call MPI_GET_ADDRESS(b % wgt, bank_disp(1), mpi_err) - call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err) - call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err) - call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err) + call MPI_GET_ADDRESS(b%wgt, bank_disp(1), mpi_err) + call MPI_GET_ADDRESS(b%xyz, bank_disp(2), mpi_err) + call MPI_GET_ADDRESS(b%uvw, bank_disp(3), mpi_err) + call MPI_GET_ADDRESS(b%E, bank_disp(4), mpi_err) ! Adjust displacements bank_disp = bank_disp - bank_disp(1) @@ -239,8 +240,8 @@ contains ! 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) + 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 @@ -274,6 +275,7 @@ contains 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 integer(HSIZE_T) :: dims(1) = (/3/) ! size of coordinates @@ -318,8 +320,8 @@ contains integer :: argc ! number of command line arguments integer :: last_flag ! index of last flag integer :: filetype + integer(HID_T) :: file_id character(MAX_WORD_LEN), allocatable :: argv(:) ! command line arguments - type(BinaryOutput) :: sp ! Check number of command line arguments and allocate argv argc = COMMAND_ARGUMENT_COUNT() @@ -356,9 +358,9 @@ contains i = i + 1 ! Check what type of file this is - call sp % file_open(argv(i), 'r', serial = .false.) - call sp % read_data(filetype, 'filetype') - call sp % file_close() + file_id = file_open(argv(i), 'r', parallel=.true.) + call read_dataset(file_id, 'filetype', filetype) + call file_close(file_id) ! Set path and flag for type of run select case (filetype) @@ -379,13 +381,12 @@ contains i = i + 1 ! Check if it has extension we can read - if ((ends_with(argv(i), '.binary') .or. & - ends_with(argv(i), '.h5'))) then + if (ends_with(argv(i), '.h5')) then ! Check file type is a source file - call sp % file_open(argv(i), 'r', serial = .false.) - call sp % read_data(filetype, 'filetype') - call sp % file_close() + file_id = file_open(argv(i), 'r', parallel=.true.) + call read_dataset(file_id, 'filetype', filetype) + call file_close(file_id) if (filetype /= FILETYPE_SOURCE) then call fatal_error("Second file after restart flag must be a & &source file") @@ -494,26 +495,26 @@ contains ! pairs are the id of the universe and the index in the array. In ! cells_in_univ_dict, it's the id of the universe and the number of cells. - pair_list => universe_dict % keys() + pair_list => universe_dict%keys() current => pair_list do while (associated(current)) ! Find index of universe in universes array - i_univ = current % value + i_univ = current%value univ => universes(i_univ) - univ % id = current % key + univ%id = current%key ! Check for lowest level universe - if (univ % id == 0) BASE_UNIVERSE = i_univ + if (univ%id == 0) BASE_UNIVERSE = i_univ ! Find cell count for this universe - n_cells_in_univ = cells_in_univ_dict % get_key(univ % id) + n_cells_in_univ = cells_in_univ_dict%get_key(univ%id) ! Allocate cell list for universe - allocate(univ % cells(n_cells_in_univ)) - univ % n_cells = n_cells_in_univ + allocate(univ%cells(n_cells_in_univ)) + univ%n_cells = n_cells_in_univ ! Move to next universe - next => current % next + next => current%next deallocate(current) current => next end do @@ -528,17 +529,17 @@ contains c => cells(i) ! Get pointer to corresponding universe - i_univ = universe_dict % get_key(c % universe) + i_univ = universe_dict%get_key(c%universe) univ => universes(i_univ) ! Increment the index for the cells array within the Universe object and ! then store the index of the Cell object in that array index_cell_in_univ(i_univ) = index_cell_in_univ(i_univ) + 1 - univ % cells(index_cell_in_univ(i_univ)) = i + univ%cells(index_cell_in_univ(i_univ)) = i end do ! Clear dictionary - call cells_in_univ_dict % clear() + call cells_in_univ_dict%clear() end subroutine prepare_universes @@ -568,15 +569,15 @@ contains ! ADJUST SURFACE LIST FOR EACH CELL c => cells(i) - do j = 1, c % n_surfaces - id = c % surfaces(j) + do j = 1, c%n_surfaces + id = c%surfaces(j) if (id < OP_DIFFERENCE) then - if (surface_dict % has_key(abs(id))) then - i_array = surface_dict % get_key(abs(id)) - c % surfaces(j) = sign(i_array, id) + if (surface_dict%has_key(abs(id))) then + i_array = surface_dict%get_key(abs(id)) + c%surfaces(j) = sign(i_array, id) else call fatal_error("Could not find surface " // trim(to_str(abs(id)))& - &// " specified on cell " // trim(to_str(c % id))) + &// " specified on cell " // trim(to_str(c%id))) end if end if end do @@ -584,40 +585,40 @@ contains ! ======================================================================= ! ADJUST UNIVERSE INDEX FOR EACH CELL - id = c % universe - if (universe_dict % has_key(id)) then - c % universe = universe_dict % get_key(id) + id = c%universe + if (universe_dict%has_key(id)) then + c%universe = universe_dict%get_key(id) else call fatal_error("Could not find universe " // trim(to_str(id)) & - &// " specified on cell " // trim(to_str(c % id))) + &// " specified on cell " // trim(to_str(c%id))) end if ! ======================================================================= ! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL - id = c % material + id = c%material if (id == MATERIAL_VOID) then - c % type = CELL_NORMAL + c%type = CELL_NORMAL elseif (id /= 0) then - if (material_dict % has_key(id)) then - c % type = CELL_NORMAL - c % material = material_dict % get_key(id) + if (material_dict%has_key(id)) then + c%type = CELL_NORMAL + c%material = material_dict%get_key(id) else call fatal_error("Could not find material " // trim(to_str(id)) & - &// " specified on cell " // trim(to_str(c % id))) + &// " specified on cell " // trim(to_str(c%id))) end if else - id = c % fill - if (universe_dict % has_key(id)) then - c % type = CELL_FILL - c % fill = universe_dict % get_key(id) - elseif (lattice_dict % has_key(id)) then - lid = lattice_dict % get_key(id) - c % type = CELL_LATTICE - c % fill = lid + id = c%fill + if (universe_dict%has_key(id)) then + c%type = CELL_FILL + c%fill = universe_dict%get_key(id) + elseif (lattice_dict%has_key(id)) then + lid = lattice_dict%get_key(id) + c%type = CELL_LATTICE + c%fill = lid else call fatal_error("Specified fill " // trim(to_str(id)) // " on cell "& - &// trim(to_str(c % id)) // " is neither a universe nor a & + &// trim(to_str(c%id)) // " is neither a universe nor a & &lattice.") end if end if @@ -627,41 +628,41 @@ contains ! ADJUST UNIVERSE INDICES FOR EACH LATTICE do i = 1, n_lattices - lat => lattices(i) % obj + lat => lattices(i)%obj select type (lat) type is (RectLattice) - do m = 1, lat % n_cells(3) - do k = 1, lat % n_cells(2) - do j = 1, lat % n_cells(1) - id = lat % universes(j,k,m) - if (universe_dict % has_key(id)) then - lat % universes(j,k,m) = universe_dict % get_key(id) + do m = 1, lat%n_cells(3) + do k = 1, lat%n_cells(2) + do j = 1, lat%n_cells(1) + id = lat%universes(j,k,m) + if (universe_dict%has_key(id)) then + lat%universes(j,k,m) = universe_dict%get_key(id) else call fatal_error("Invalid universe number " & &// trim(to_str(id)) // " specified on lattice " & - &// trim(to_str(lat % id))) + &// trim(to_str(lat%id))) end if end do end do end do type is (HexLattice) - do m = 1, lat % n_axial - do k = 1, 2*lat % n_rings - 1 - do j = 1, 2*lat % n_rings - 1 - if (j + k < lat % n_rings + 1) then + do m = 1, lat%n_axial + do k = 1, 2*lat%n_rings - 1 + do j = 1, 2*lat%n_rings - 1 + if (j + k < lat%n_rings + 1) then cycle - else if (j + k > 3*lat % n_rings - 1) then + else if (j + k > 3*lat%n_rings - 1) then cycle end if - id = lat % universes(j, k, m) - if (universe_dict % has_key(id)) then - lat % universes(j, k, m) = universe_dict % get_key(id) + id = lat%universes(j, k, m) + if (universe_dict%has_key(id)) then + lat%universes(j, k, m) = universe_dict%get_key(id) else call fatal_error("Invalid universe number " & &// trim(to_str(id)) // " specified on lattice " & - &// trim(to_str(lat % id))) + &// trim(to_str(lat%id))) end if end do end do @@ -669,13 +670,13 @@ contains end select - if (lat % outer /= NO_OUTER_UNIVERSE) then - if (universe_dict % has_key(lat % outer)) then - lat % outer = universe_dict % get_key(lat % outer) + if (lat%outer /= NO_OUTER_UNIVERSE) then + if (universe_dict%has_key(lat%outer)) then + lat%outer = universe_dict%get_key(lat%outer) else call fatal_error("Invalid universe number " & - &// trim(to_str(lat % outer)) & - &// " specified on lattice " // trim(to_str(lat % id))) + &// trim(to_str(lat%outer)) & + &// " specified on lattice " // trim(to_str(lat%id))) end if end if @@ -687,68 +688,68 @@ contains ! ======================================================================= ! ADJUST INDICES FOR EACH TALLY FILTER - FILTER_LOOP: do j = 1, t % n_filters + FILTER_LOOP: do j = 1, t%n_filters - select case (t % filters(j) % type) + select case (t%filters(j)%type) case (FILTER_DISTRIBCELL) - do k = 1, size(t % filters(j) % int_bins) - id = t % filters(j) % int_bins(k) - if (cell_dict % has_key(id)) then - t % filters(j) % int_bins(k) = cell_dict % get_key(id) + do k = 1, size(t%filters(j)%int_bins) + id = t%filters(j)%int_bins(k) + if (cell_dict%has_key(id)) then + t%filters(j)%int_bins(k) = cell_dict%get_key(id) else call fatal_error("Could not find cell " // trim(to_str(id)) // & - " specified on tally " // trim(to_str(t % id))) + " specified on tally " // trim(to_str(t%id))) end if end do case (FILTER_CELL, FILTER_CELLBORN) - do k = 1, t % filters(j) % n_bins - id = t % filters(j) % int_bins(k) - if (cell_dict % has_key(id)) then - t % filters(j) % int_bins(k) = cell_dict % get_key(id) + do k = 1, t%filters(j)%n_bins + id = t%filters(j)%int_bins(k) + if (cell_dict%has_key(id)) then + t%filters(j)%int_bins(k) = cell_dict%get_key(id) else call fatal_error("Could not find cell " // trim(to_str(id)) & - &// " specified on tally " // trim(to_str(t % id))) + &// " specified on tally " // trim(to_str(t%id))) end if end do case (FILTER_SURFACE) ! Check if this is a surface filter only for surface currents - if (any(t % score_bins == SCORE_CURRENT)) cycle FILTER_LOOP + if (any(t%score_bins == SCORE_CURRENT)) cycle FILTER_LOOP - do k = 1, t % filters(j) % n_bins - id = t % filters(j) % int_bins(k) - if (surface_dict % has_key(id)) then - t % filters(j) % int_bins(k) = surface_dict % get_key(id) + do k = 1, t%filters(j)%n_bins + id = t%filters(j)%int_bins(k) + if (surface_dict%has_key(id)) then + t%filters(j)%int_bins(k) = surface_dict%get_key(id) else call fatal_error("Could not find surface " // trim(to_str(id)) & - &// " specified on tally " // trim(to_str(t % id))) + &// " specified on tally " // trim(to_str(t%id))) end if end do case (FILTER_UNIVERSE) - do k = 1, t % filters(j) % n_bins - id = t % filters(j) % int_bins(k) - if (universe_dict % has_key(id)) then - t % filters(j) % int_bins(k) = universe_dict % get_key(id) + do k = 1, t%filters(j)%n_bins + id = t%filters(j)%int_bins(k) + if (universe_dict%has_key(id)) then + t%filters(j)%int_bins(k) = universe_dict%get_key(id) else call fatal_error("Could not find universe " // trim(to_str(id)) & - &// " specified on tally " // trim(to_str(t % id))) + &// " specified on tally " // trim(to_str(t%id))) end if end do case (FILTER_MATERIAL) - do k = 1, t % filters(j) % n_bins - id = t % filters(j) % int_bins(k) - if (material_dict % has_key(id)) then - t % filters(j) % int_bins(k) = material_dict % get_key(id) + do k = 1, t%filters(j)%n_bins + id = t%filters(j)%int_bins(k) + if (material_dict%has_key(id)) then + t%filters(j)%int_bins(k) = material_dict%get_key(id) else call fatal_error("Could not find material " // trim(to_str(id)) & - &// " specified on tally " // trim(to_str(t % id))) + &// " specified on tally " // trim(to_str(t%id))) end if end do @@ -786,46 +787,46 @@ contains do i = 1, n_materials mat => materials(i) - percent_in_atom = (mat % atom_density(1) > ZERO) - density_in_atom = (mat % density > ZERO) + percent_in_atom = (mat%atom_density(1) > ZERO) + density_in_atom = (mat%density > ZERO) sum_percent = ZERO - do j = 1, mat % n_nuclides + do j = 1, mat%n_nuclides ! determine atomic weight ratio - index_list = xs_listing_dict % get_key(mat % names(j)) - awr = xs_listings(index_list) % awr + index_list = xs_listing_dict%get_key(mat%names(j)) + awr = xs_listings(index_list)%awr ! if given weight percent, convert all values so that they are divided ! by awr. thus, when a sum is done over the values, it's actually ! sum(w/awr) if (.not. percent_in_atom) then - mat % atom_density(j) = -mat % atom_density(j) / awr + mat%atom_density(j) = -mat%atom_density(j) / awr end if end do ! determine normalized atom percents. if given atom percents, this is ! straightforward. if given weight percents, the value is w/awr and is ! divided by sum(w/awr) - sum_percent = sum(mat % atom_density) - mat % atom_density = mat % atom_density / sum_percent + sum_percent = sum(mat%atom_density) + mat%atom_density = mat%atom_density / sum_percent ! Change density in g/cm^3 to atom/b-cm. Since all values are now in atom ! percent, the sum needs to be re-evaluated as 1/sum(x*awr) if (.not. density_in_atom) then sum_percent = ZERO - do j = 1, mat % n_nuclides - index_list = xs_listing_dict % get_key(mat % names(j)) - awr = xs_listings(index_list) % awr - x = mat % atom_density(j) + do j = 1, mat%n_nuclides + index_list = xs_listing_dict%get_key(mat%names(j)) + awr = xs_listings(index_list)%awr + x = mat%atom_density(j) sum_percent = sum_percent + x*awr end do sum_percent = ONE / sum_percent - mat % density = -mat % density * N_AVOGADRO & + mat%density = -mat%density * N_AVOGADRO & / MASS_NEUTRON * sum_percent end if ! Calculate nuclide atom densities - mat % atom_density = mat % density * mat % atom_density + mat%atom_density = mat%density * mat%atom_density end do end subroutine normalize_ao @@ -942,16 +943,16 @@ contains ! Get pointer to tally tally => tallies(i) - n_filt = tally % n_filters + n_filt = tally%n_filters ! Loop over the filters to determine how many additional filters ! need to be added to this tally - do j = 1, tally % n_filters + do j = 1, tally%n_filters ! Determine type of filter - if (tally % filters(j) % type == FILTER_DISTRIBCELL) then + if (tally%filters(j)%type == FILTER_DISTRIBCELL) then count_all = .true. - if (size(tally % filters(j) % int_bins) > 1) then + if (size(tally%filters(j)%int_bins) > 1) then call fatal_error("A distribcell filter was specified with & &multiple bins. This feature is not supported.") end if @@ -975,12 +976,12 @@ contains tally => tallies(i) ! Initialize the filters - do j = 1, tally % n_filters + do j = 1, tally%n_filters ! Set the number of bins to the number of instances of the cell - if (tally % filters(j) % type == FILTER_DISTRIBCELL) then - c => cells(tally % filters(j) % int_bins(1)) - tally % filters(j) % n_bins = c % instances + if (tally%filters(j)%type == FILTER_DISTRIBCELL) then + c => cells(tally%filters(j)%int_bins(1)) + tally%filters(j)%n_bins = c%instances end if end do @@ -1031,12 +1032,12 @@ contains do i = 1, n_tallies tally => tallies(i) - do j = 1, tally % n_filters - filter => tally % filters(j) + do j = 1, tally%n_filters + filter => tally%filters(j) - if (filter % type == FILTER_DISTRIBCELL) then - if (.not. cell_list % contains(filter % int_bins(1))) then - call cell_list % add(filter % int_bins(1)) + if (filter%type == FILTER_DISTRIBCELL) then + if (.not. cell_list%contains(filter%int_bins(1))) then + call cell_list%add(filter%int_bins(1)) end if end if @@ -1047,8 +1048,8 @@ contains ! to determine the number of offset tables to allocate do i = 1, n_universes univ => universes(i) - do j = 1, univ % n_cells - if (cell_list % contains(univ % cells(j))) then + do j = 1, univ%n_cells + if (cell_list%contains(univ%cells(j))) then n_maps = n_maps + 1 end if end do @@ -1070,29 +1071,29 @@ contains do i = 1, n_universes univ => universes(i) - do j = 1, univ % n_cells + do j = 1, univ%n_cells - if (cell_list % contains(univ % cells(j))) then + if (cell_list%contains(univ%cells(j))) then ! Loop over all tallies do l = 1, n_tallies tally => tallies(l) - do m = 1, tally % n_filters - filter => tally % filters(m) + do m = 1, tally%n_filters + filter => tally%filters(m) ! Loop over only distribcell filters ! If filter points to cell we just found, set offset index - if (filter % type == FILTER_DISTRIBCELL) then - if (filter % int_bins(1) == univ % cells(j)) then - filter % offset = k + if (filter%type == FILTER_DISTRIBCELL) then + if (filter%int_bins(1) == univ%cells(j)) then + filter%offset = k end if end if end do end do - univ_list(k) = univ % id + univ_list(k) = univ%id k = k + 1 end if end do @@ -1100,26 +1101,26 @@ contains ! Allocate the offset tables for lattices do i = 1, n_lattices - lat => lattices(i) % obj + lat => lattices(i)%obj select type(lat) type is (RectLattice) - allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), & - lat % n_cells(3))) + allocate(lat%offset(n_maps, lat%n_cells(1), lat%n_cells(2), & + lat%n_cells(3))) type is (HexLattice) - allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, & - 2 * lat % n_rings - 1, lat % n_axial)) + allocate(lat%offset(n_maps, 2 * lat%n_rings - 1, & + 2 * lat%n_rings - 1, lat%n_axial)) end select - lat % offset(:, :, :, :) = 0 + lat%offset(:, :, :, :) = 0 end do ! Allocate offset table for fill cells do i = 1, n_cells - if (cells(i) % material == NONE) then - allocate(cells(i) % offset(n_maps)) + if (cells(i)%material == NONE) then + allocate(cells(i)%offset(n_maps)) end if end do diff --git a/src/output_interface.F90 b/src/output_interface.F90 deleted file mode 100644 index 1e08ce4911..0000000000 --- a/src/output_interface.F90 +++ /dev/null @@ -1,1725 +0,0 @@ -module output_interface - - use constants - use error, only: warning, fatal_error - use global - use tally_header, only: TallyResult - - use hdf5_interface - - implicit none - private - - type, public :: BinaryOutput - private - ! Compilation specific data - integer(HID_T) :: hdf5_fh - integer(HID_T) :: hdf5_grp - logical :: serial ! Serial I/O when using MPI/PHDF5 - contains - generic, public :: write_data => write_double, & - write_double_1Darray, & - write_double_2Darray, & - write_double_3Darray, & - write_double_4Darray, & - write_integer, & - write_integer_1Darray, & - write_integer_2Darray, & - write_integer_3Darray, & - write_integer_4Darray, & - write_long, & - write_string - generic, public :: read_data => read_double, & - read_double_1Darray, & - read_double_2Darray, & - read_double_3Darray, & - read_double_4Darray, & - read_integer, & - read_integer_1Darray, & - read_integer_2Darray, & - read_integer_3Darray, & - read_integer_4Darray, & - read_long, & - read_string - procedure :: write_double => write_double - procedure :: write_double_1Darray => write_double_1Darray - procedure :: write_double_2Darray => write_double_2Darray - procedure :: write_double_3Darray => write_double_3Darray - procedure :: write_double_4Darray => write_double_4Darray - procedure :: write_integer => write_integer - procedure :: write_integer_1Darray => write_integer_1Darray - procedure :: write_integer_2Darray => write_integer_2Darray - procedure :: write_integer_3Darray => write_integer_3Darray - procedure :: write_integer_4Darray => write_integer_4Darray - procedure :: write_long => write_long - procedure :: write_string => write_string - procedure :: read_double => read_double - procedure :: read_double_1Darray => read_double_1Darray - procedure :: read_double_2Darray => read_double_2Darray - procedure :: read_double_3Darray => read_double_3Darray - procedure :: read_double_4Darray => read_double_4Darray - procedure :: read_integer => read_integer - procedure :: read_integer_1Darray => read_integer_1Darray - procedure :: read_integer_2Darray => read_integer_2Darray - procedure :: read_integer_3Darray => read_integer_3Darray - procedure :: read_integer_4Darray => read_integer_4Darray - procedure :: read_long => read_long - procedure :: read_string => read_string - procedure, public :: file_create => file_create - procedure, public :: file_open => file_open - procedure, public :: file_close => file_close - procedure, public :: write_tally_result => write_tally_result - procedure, public :: read_tally_result => read_tally_result - procedure, public :: write_source_bank => write_source_bank - procedure, public :: read_source_bank => read_source_bank - procedure, public :: write_attribute_string => write_attribute_string - procedure, public :: open_group => open_group - procedure, public :: close_group => close_group - end type BinaryOutput - -contains - -!=============================================================================== -! FILE_CREATE creates a new file to write data to -!=============================================================================== - - subroutine file_create(self, filename, serial) - - character(*), intent(in) :: filename ! name of file to be created - logical, optional, intent(in) :: serial ! processor rank to write from - class(BinaryOutput) :: self - - ! Check for serial option - if (present(serial)) then - self % serial = serial - else - self % serial = .true. - end if - -#ifdef PHDF5 - if (self % serial) then - call hdf5_file_create(filename, self % hdf5_fh) - else - call hdf5_file_create_parallel(filename, self % hdf5_fh) - endif -#else - call hdf5_file_create(filename, self % hdf5_fh) -#endif - - end subroutine file_create - -!=============================================================================== -! FILE_OPEN opens an existing file for reading or read/writing -!=============================================================================== - - subroutine file_open(self, filename, mode, serial) - - character(*), intent(in) :: filename ! name of file to be opened - character(*), intent(in) :: mode ! file access mode - logical, optional, intent(in) :: serial ! processor rank to write from - class(BinaryOutput) :: self - - ! Check for serial option - if (present(serial)) then - self % serial = serial - else - self % serial = .true. - end if - -#ifdef PHDF5 - if (self % serial) then - call hdf5_file_open(filename, self % hdf5_fh, mode) - else - call hdf5_file_open_parallel(filename, self % hdf5_fh, mode) - endif -#else - call hdf5_file_open(filename, self % hdf5_fh, mode) -#endif - - end subroutine file_open - -!=============================================================================== -! FILE_CLOSE closes a file -!=============================================================================== - - subroutine file_close(self) - - class(BinaryOutput) :: self - - call hdf5_file_close(self % hdf5_fh) - - end subroutine file_close - -!=============================================================================== -! OPEN_GROUP call hdf5 routine to open a group within binary output context -!=============================================================================== - - subroutine open_group(self, group) - - character(*), intent(in) :: group ! HDF5 group name - class(BinaryOutput) :: self - - call hdf5_open_group(self % hdf5_fh, group, self % hdf5_grp) - - end subroutine open_group - -!=============================================================================== -! CLOSE_GROUP call hdf5 routine to close a group within binary output context -!=============================================================================== - - subroutine close_group(self) - - class(BinaryOutput) :: self - - call hdf5_close_group(self % hdf5_grp) - - end subroutine close_group - -!=============================================================================== -! WRITE_DOUBLE writes double precision scalar data -!=============================================================================== - - subroutine write_double(self, buffer, name, group, collect) - - real(8), intent(in) :: buffer ! data to write - character(*), intent(in) :: name ! name for data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_write_double(self % hdf5_grp, name_, buffer) - else - call hdf5_write_double_parallel(self % hdf5_grp, name_, buffer, collect_) - end if -#else - call hdf5_write_double(self % hdf5_grp, name_, buffer) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_double - -!=============================================================================== -! READ_DOUBLE reads double precision scalar data -!=============================================================================== - - subroutine read_double(self, buffer, name, group, collect) - - real(8), intent(inout) :: buffer ! read data to here - character(*), intent(in) :: name ! name for data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_read_double(self % hdf5_grp, name_, buffer) - else - call hdf5_read_double_parallel(self % hdf5_grp, name_, buffer, collect_) - end if -#else - call hdf5_read_double(self % hdf5_grp, name_, buffer) -#endif - ! Check if HDf5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_double - -!=============================================================================== -! WRITE_DOUBLE_1DARRAY writes double precision 1-D array data -!=============================================================================== - - subroutine write_double_1Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length ! length of array to write - real(8), intent(in) :: buffer(:) ! data to write - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_write_double_1Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_write_double_1Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - call hdf5_write_double_1Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_double_1Darray - -!=============================================================================== -! READ_DOUBLE_1DARRAY reads double precision 1-D array data -!=============================================================================== - - subroutine read_double_1Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length ! length of array to read - real(8), intent(inout) :: buffer(:) ! read data to here - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_read_double_1Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_read_double_1Darray_parallel(self % hdf5_grp, name_, buffer, & - length, collect_) - end if -#else - call hdf5_read_double_1Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_double_1Darray - -!=============================================================================== -! WRITE_DOUBLE_2DARRAY writes double precision 2-D array data -!=============================================================================== - - subroutine write_double_2Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length(2) ! dimension of array - real(8), intent(in) :: buffer(length(1),length(2)) ! the data - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_write_double_2Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_write_double_2Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - call hdf5_write_double_2Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_double_2Darray - -!=============================================================================== -! READ_DOUBLE_2DARRAY reads double precision 2-D array data -!=============================================================================== - - subroutine read_double_2Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length(2) ! dimension of array - real(8), intent(inout) :: buffer(length(1),length(2)) ! the data - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_read_double_2Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_read_double_2Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - call hdf5_read_double_2Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_double_2Darray - -!=============================================================================== -! WRITE_DOUBLE_3DARRAY writes double precision 3-D array data -!=============================================================================== - - subroutine write_double_3Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length(3) ! length of each dimension - real(8), intent(in) :: buffer(length(1),length(2),length(3)) - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_write_double_3Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_write_double_3Darray_parallel(self % hdf5_grp, name_, buffer, & - length, collect_) - end if -#else - call hdf5_write_double_3Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_double_3Darray - -!=============================================================================== -! READ_DOUBLE_3DARRAY reads double precision 3-D array data -!=============================================================================== - - subroutine read_double_3Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length(3) ! length of each dimension - real(8), intent(inout) :: buffer(length(1),length(2),length(3)) - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_read_double_3Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_read_double_3Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - call hdf5_read_double_3Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_double_3Darray - -!=============================================================================== -! WRITE_DOUBLE_4DARRAY writes double precision 4-D array data -!=============================================================================== - - subroutine write_double_4Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length(4) ! length of each dimension - real(8), intent(in) :: buffer(length(1),length(2),& - length(3),length(4)) - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_write_double_4Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_write_double_4Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - ! Write the data in serial - call hdf5_write_double_4Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_double_4Darray - -!=============================================================================== -! READ_DOUBLE_4DARRAY reads double precision 4-D array data -!=============================================================================== - - subroutine read_double_4Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length(4) ! length of each dimension - real(8), intent(inout) :: buffer(length(1),length(2),& - length(3),length(4)) - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_read_double_4Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_read_double_4Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - call hdf5_read_double_4Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_double_4Darray - -!=============================================================================== -! WRITE_INTEGER writes integer precision scalar data -!=============================================================================== - - subroutine write_integer(self, buffer, name, group, collect) - - integer, intent(in) :: buffer ! data to write - character(*), intent(in) :: name ! name for data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_write_integer(self % hdf5_grp, name_, buffer) - else - call hdf5_write_integer_parallel(self % hdf5_grp, name_, buffer, collect_) - end if -#else - call hdf5_write_integer(self % hdf5_grp, name_, buffer) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_integer - -!=============================================================================== -! READ_INTEGER reads integer precision scalar data -!=============================================================================== - - subroutine read_integer(self, buffer, name, group, collect) - - integer, intent(inout) :: buffer ! read data to here - character(*), intent(in) :: name ! name for data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_read_integer(self % hdf5_grp, name_, buffer) - else - call hdf5_read_integer_parallel(self % hdf5_grp, name_, buffer, collect_) - end if -#else - call hdf5_read_integer(self % hdf5_grp, name_, buffer) -#endif - ! Check if HDf5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_integer - -!=============================================================================== -! WRITE_INTEGER_1DARRAY writes integer precision 1-D array data -!=============================================================================== - - subroutine write_integer_1Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length ! length of array to write - integer, intent(in) :: buffer(:) ! data to write - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_write_integer_1Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_write_integer_1Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - call hdf5_write_integer_1Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_integer_1Darray - -!=============================================================================== -! READ_INTEGER_1DARRAY reads integer precision 1-D array data -!=============================================================================== - - subroutine read_integer_1Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length ! length of array to read - integer, intent(inout) :: buffer(:) ! read data to here - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_read_integer_1Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_read_integer_1Darray_parallel(self % hdf5_grp, name_, buffer, & - length, collect_) - end if -#else - ! Read the data in serial - call hdf5_read_integer_1Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_integer_1Darray - -!=============================================================================== -! WRITE_INTEGER_2DARRAY writes integer precision 2-D array data -!=============================================================================== - - subroutine write_integer_2Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length(2) ! dimension of array - integer, intent(in) :: buffer(length(1),length(2)) ! the data - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_write_integer_2Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_write_integer_2Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - call hdf5_write_integer_2Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_integer_2Darray - -!=============================================================================== -! READ_INTEGER_2DARRAY reads integer precision 2-D array data -!=============================================================================== - - subroutine read_integer_2Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length(2) ! dimension of array - integer, intent(inout) :: buffer(length(1),length(2)) ! the data - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_read_integer_2Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_read_integer_2Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - call hdf5_read_integer_2Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_integer_2Darray - -!=============================================================================== -! WRITE_INTEGER_3DARRAY writes integer precision 3-D array data -!=============================================================================== - - subroutine write_integer_3Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length(3) ! length of each dimension - integer, intent(in) :: buffer(length(1),length(2),length(3)) - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_write_integer_3Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_write_integer_3Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - call hdf5_write_integer_3Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_integer_3Darray - -!=============================================================================== -! READ_INTEGER_3DARRAY reads integer precision 3-D array data -!=============================================================================== - - subroutine read_integer_3Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length(3) ! length of each dimension - integer, intent(inout) :: buffer(length(1),length(2),length(3)) - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_read_integer_3Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_read_integer_3Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - call hdf5_read_integer_3Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_integer_3Darray - -!=============================================================================== -! WRITE_INTEGER_4DARRAY writes integer precision 4-D array data -!=============================================================================== - - subroutine write_integer_4Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length(4) ! length of each dimension - integer, intent(in) :: buffer(length(1),length(2),& - length(3),length(4)) - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_write_integer_4Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_write_integer_4Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - call hdf5_write_integer_4Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_integer_4Darray - -!=============================================================================== -! READ_INTEGER_4DARRAY reads integer precision 4-D array data -!=============================================================================== - - subroutine read_integer_4Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length(4) ! length of each dimension - integer, intent(inout) :: buffer(length(1),length(2),& - length(3),length(4)) - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_read_integer_4Darray(self % hdf5_grp, name_, buffer, length) - else - call hdf5_read_integer_4Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) - end if -#else - call hdf5_read_integer_4Darray(self % hdf5_grp, name_, buffer, length) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_integer_4Darray - -!=============================================================================== -! WRITE_LONG writes long integer scalar data -!=============================================================================== - - subroutine write_long(self, buffer, name, group, collect) - - integer(8), intent(in) :: buffer ! data to write - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_write_long(self % hdf5_grp, name_, buffer, hdf5_integer8_t) - else - call hdf5_write_long_parallel(self % hdf5_grp, name_, buffer, & - hdf5_integer8_t, collect_) - end if -#else - call hdf5_write_long(self % hdf5_grp, name_, buffer, hdf5_integer8_t) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_long - -!=============================================================================== -! READ_LONG reads long integer scalar data -!=============================================================================== - - subroutine read_long(self, buffer, name, group, collect) - - integer(8), intent(inout) :: buffer ! data to write - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - logical :: collect_ - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_read_long(self % hdf5_grp, name_, buffer, hdf5_integer8_t) - else - call hdf5_read_long_parallel(self % hdf5_grp, name_, buffer, & - hdf5_integer8_t, collect_) - end if -#else - call hdf5_read_long(self % hdf5_grp, name_, buffer, hdf5_integer8_t) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_long - -!=============================================================================== -! WRITE_STRING writes string data -!=============================================================================== - - subroutine write_string(self, buffer, name, group, collect) - - character(*), intent(in) :: buffer ! data to write - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - integer :: n - logical :: collect_ - - ! Get string length - n = len_trim(buffer) - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_write_string(self % hdf5_grp, name_, buffer, n) - else - call hdf5_write_string_parallel(self % hdf5_grp, name_, buffer, n, collect_) - end if -#else - ! Write the data - call hdf5_write_string(self % hdf5_grp, name_, buffer, n) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_string - -!=============================================================================== -! READ_STRING reads string data -!=============================================================================== - - subroutine read_string(self, buffer, name, group, collect) - - character(*), intent(inout) :: buffer ! data to write - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - integer :: n - logical :: collect_ - - ! Get string length - n = len(buffer) - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif -#ifdef PHDF5 - if (self % serial) then - call hdf5_read_string(self % hdf5_grp, name_, buffer, n) - else - call hdf5_read_string_parallel(self % hdf5_grp, name_, buffer, n, collect_) - end if -#else - call hdf5_read_string(self % hdf5_grp, name_, buffer, n) -#endif - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_string - -!=============================================================================== -! WRITE_ATTRIBUTE_STRING -!=============================================================================== - - subroutine write_attribute_string(self, var, attr_type, attr_str, group) - - character(*), intent(in) :: var ! variable name for attr - character(*), intent(in) :: attr_type ! attr identifier type - character(*), intent(in) :: attr_str ! string for attr id type - character(*), intent(in), optional :: group ! HDF5 group name - class(BinaryOutput) :: self - - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif - - ! Write the attribute string - call hdf5_write_attribute_string(self % hdf5_grp, var, attr_type, attr_str) - - ! Check if HDF5 group should be closed - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine write_attribute_string - -!=============================================================================== -! WRITE_TALLY_RESULT writes an OpenMC TallyResult type -!=============================================================================== - - subroutine write_tally_result(self, buffer, name, group, n1, n2) - - character(*), intent(in), optional :: group ! HDF5 group name - character(*), intent(in) :: name ! name of data - integer, intent(in) :: n1, n2 ! TallyResult dims - type(TallyResult), intent(in), target :: buffer(n1, n2) ! data to write - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Open up sub-group if present - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - end if - - ! Set overall size of vector to write - dims1(1) = n1*n2 - - ! Create up a dataspace for size - call h5screate_simple_f(1, dims1, dspace, hdf5_err) - - ! Create the dataset - call h5dcreate_f(self % hdf5_grp, name_, hdf5_tallyresult_t, dspace, dset, & - hdf5_err) - - ! Set pointer to first value and write - f_ptr = c_loc(buffer(1,1)) - call h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - - ! Close ids - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - if (present(group)) then - call hdf5_close_group(self % hdf5_grp) - end if - - end subroutine write_tally_result - -!=============================================================================== -! READ_TALLY_RESULT reads OpenMC TallyResult data -!=============================================================================== - - subroutine read_tally_result(self, buffer, name, group, n1, n2) - - character(*), intent(in), optional :: group ! HDF5 group name - character(*), intent(in) :: name ! name of data - integer, intent(in) :: n1, n2 ! TallyResult dims - type(TallyResult), intent(inout), target :: buffer(n1, n2) ! read data here - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Open up sub-group if present - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - end if - - ! Open the dataset - call h5dopen_f(self % hdf5_grp, name, dset, hdf5_err) - - ! Set pointer to first value and write - f_ptr = c_loc(buffer(1,1)) - call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - - ! Close ids - call h5dclose_f(dset, hdf5_err) - if (present(group)) call hdf5_close_group(self % hdf5_grp) - - end subroutine read_tally_result - -!=============================================================================== -! WRITE_SOURCE_BANK writes OpenMC source_bank data -!=============================================================================== - - subroutine write_source_bank(self) - - class(BinaryOutput) :: self - -#ifdef PHDF5 - integer(8) :: offset(1) ! source data offset -#endif - -#ifdef PHDF5 - - ! Set size of total dataspace for all procs and rank - dims1(1) = n_particles - hdf5_rank = 1 - - ! Create that dataspace - call h5screate_simple_f(hdf5_rank, dims1, dspace, hdf5_err) - - ! Create the dataset for that dataspace - call h5dcreate_f(self % hdf5_fh, "source_bank", hdf5_bank_t, dspace, dset, hdf5_err) - - ! Close the dataspace - call h5sclose_f(dspace, hdf5_err) - - ! Create another data space but for each proc individually - dims1(1) = work - call h5screate_simple_f(hdf5_rank, dims1, memspace, hdf5_err) - - ! Get the individual local proc dataspace - call h5dget_space_f(dset, dspace, hdf5_err) - - ! Select hyperslab for this dataspace - offset(1) = work_index(rank) - call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims1, hdf5_err) - - ! Set up the property list for parallel writing - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, H5FD_MPIO_COLLECTIVE_F, hdf5_err) - - ! Set up pointer to data - f_ptr = c_loc(source_bank(1)) - - ! Write data to file in parallel - call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err, & - file_space_id = dspace, mem_space_id = memspace, & - xfer_prp = plist) - - ! Close all ids - call h5sclose_f(dspace, hdf5_err) - call h5sclose_f(memspace, hdf5_err) - call h5dclose_f(dset, hdf5_err) - call h5pclose_f(plist, hdf5_err) - -#else - - ! Set size - dims1(1) = work - hdf5_rank = 1 - - ! Create dataspace - call h5screate_simple_f(hdf5_rank, dims1, dspace, hdf5_err) - - ! Create dataset - call h5dcreate_f(self % hdf5_fh, "source_bank", hdf5_bank_t, & - dspace, dset, hdf5_err) - - ! Set up pointer to data - f_ptr = c_loc(source_bank(1)) - - ! Write dataset to file - call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err) - - ! Close all ids - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - -#endif - - end subroutine write_source_bank - -!=============================================================================== -! READ_SOURCE_BANK reads OpenMC source_bank data -!=============================================================================== - - subroutine read_source_bank(self) - - class(BinaryOutput) :: self - -#ifdef PHDF5 - integer(8) :: offset(1) ! offset of data -#endif - -#ifdef PHDF5 - - ! Set size of total dataspace for all procs and rank - dims1(1) = n_particles - hdf5_rank = 1 - - ! Open the dataset - call h5dopen_f(self % hdf5_fh, "source_bank", dset, hdf5_err) - - ! Create another data space but for each proc individually - dims1(1) = work - call h5screate_simple_f(hdf5_rank, dims1, memspace, hdf5_err) - - ! Get the individual local proc dataspace - call h5dget_space_f(dset, dspace, hdf5_err) - - ! Select hyperslab for this dataspace - offset(1) = work_index(rank) - call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims1, hdf5_err) - - ! Set up the property list for parallel writing - call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) - call h5pset_dxpl_mpio_f(plist, H5FD_MPIO_COLLECTIVE_F, hdf5_err) - - ! Set up pointer to data - f_ptr = c_loc(source_bank(1)) - - ! Read data from file in parallel - call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err, & - file_space_id = dspace, mem_space_id = memspace, & - xfer_prp = plist) - - ! Close all ids - call h5sclose_f(dspace, hdf5_err) - call h5sclose_f(memspace, hdf5_err) - call h5dclose_f(dset, hdf5_err) - call h5pclose_f(plist, hdf5_err) - -#else - - ! Open dataset - call h5dopen_f(self % hdf5_fh, "source_bank", dset, hdf5_err) - - ! Set up pointer to data - f_ptr = c_loc(source_bank(1)) - - ! Read dataset from file - call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err) - - ! Close all ids - call h5dclose_f(dset, hdf5_err) - -#endif - - end subroutine read_source_bank - -end module output_interface diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 0ed25fdc80..c8874b0715 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -6,18 +6,18 @@ module particle_restart use constants use geometry_header, only: BASE_UNIVERSE use global + use hdf5_interface, only: file_open, file_close, read_dataset use output, only: write_message, print_particle - use output_interface, only: BinaryOutput use particle_header, only: Particle use random_lcg, only: set_particle_seed use tracking, only: transport + use hdf5, only: HID_T + implicit none private public :: run_particle_restart - type(BinaryOutput) :: pr ! Binary file - contains !=============================================================================== @@ -34,7 +34,7 @@ contains verbosity = 10 ! Initialize the particle to be tracked - call p % initialize() + call p%initialize() ! Read in the restart information call read_particle_restart(p, previous_run_mode) @@ -46,9 +46,9 @@ contains select case (previous_run_mode) case (MODE_EIGENVALUE) particle_seed = ((current_batch - 1)*gen_per_batch + & - current_gen - 1)*n_particles + p % id + current_gen - 1)*n_particles + p%id case (MODE_FIXEDSOURCE) - particle_seed = p % id + particle_seed = p%id end select call set_particle_seed(particle_seed) @@ -66,40 +66,41 @@ contains !=============================================================================== subroutine read_particle_restart(p, previous_run_mode) + type(Particle), intent(inout) :: p + integer, intent(inout) :: previous_run_mode integer :: int_scalar - integer, intent(inout) :: previous_run_mode - type(Particle), intent(inout) :: p + integer(HID_T) :: file_id ! Write meessage call write_message("Loading particle restart file " & &// trim(path_particle_restart) // "...", 1) ! Open file - call pr % file_open(path_particle_restart, 'r') + file_id = file_open(path_particle_restart, 'r') ! Read data from file - call pr % read_data(int_scalar, 'filetype') - call pr % read_data(int_scalar, 'revision') - call pr % read_data(current_batch, 'current_batch') - call pr % read_data(gen_per_batch, 'gen_per_batch') - call pr % read_data(current_gen, 'current_gen') - call pr % read_data(n_particles, 'n_particles') - call pr % read_data(previous_run_mode, 'run_mode') - call pr % read_data(p % id, 'id') - call pr % read_data(p % wgt, 'weight') - call pr % read_data(p % E, 'energy') - call pr % read_data(p % coord(1) % xyz, 'xyz', length=3) - call pr % read_data(p % coord(1) % uvw, 'uvw', length=3) + call read_dataset(file_id, 'filetype', int_scalar) + call read_dataset(file_id, 'revision', int_scalar) + call read_dataset(file_id, 'current_batch', current_batch) + call read_dataset(file_id, 'gen_per_batch', gen_per_batch) + call read_dataset(file_id, 'current_gen', current_gen) + call read_dataset(file_id, 'n_particles', n_particles) + call read_dataset(file_id, 'run_mode', previous_run_mode) + call read_dataset(file_id, 'id', p%id) + call read_dataset(file_id, 'weight', p%wgt) + call read_dataset(file_id, 'energy', p%E) + call read_dataset(file_id, 'xyz', p%coord(1)%xyz) + call read_dataset(file_id, 'uvw', p%coord(1)%uvw) ! Set particle last attributes - p % last_wgt = p % wgt - p % last_xyz = p % coord(1) % xyz - p % last_uvw = p % coord(1) % uvw - p % last_E = p % E + p%last_wgt = p%wgt + p%last_xyz = p%coord(1)%xyz + p%last_uvw = p%coord(1)%uvw + p%last_E = p%E ! Close hdf5 file - call pr % file_close() + call file_close(file_id) end subroutine read_particle_restart diff --git a/src/particle_restart_write.F90 b/src/particle_restart_write.F90 index 48b4af661e..0c010b64c2 100644 --- a/src/particle_restart_write.F90 +++ b/src/particle_restart_write.F90 @@ -2,17 +2,16 @@ module particle_restart_write use bank_header, only: Bank use global - use output_interface, only: BinaryOutput + use hdf5_interface use particle_header, only: Particle use string, only: to_str + use hdf5 + implicit none private public :: write_particle_restart - ! Binary output file - type(BinaryOutput) :: pr - contains !=============================================================================== @@ -20,43 +19,42 @@ contains !=============================================================================== subroutine write_particle_restart(p) - type(Particle), intent(in) :: p + integer(HID_T) :: file_id character(MAX_FILE_LEN) :: filename - type(Bank), pointer :: src => null() + type(Bank), pointer :: src ! Dont write another restart file if in particle restart mode if (run_mode == MODE_PARTICLE) return ! Set up file name filename = trim(path_output) // 'particle_' // trim(to_str(current_batch)) & - // '_' // trim(to_str(p % id)) - filename = trim(filename) // '.h5' + // '_' // trim(to_str(p%id)) // '.h5' !$omp critical (WriteParticleRestart) ! Create file - call pr % file_create(filename) + file_id = file_create(filename) ! Get information about source particle src => source_bank(current_work) ! Write data to file - call pr % write_data(FILETYPE_PARTICLE_RESTART, 'filetype') - call pr % write_data(REVISION_PARTICLE_RESTART, 'revision') - call pr % write_data(current_batch, 'current_batch') - call pr % write_data(gen_per_batch, 'gen_per_batch') - call pr % write_data(current_gen, 'current_gen') - call pr % write_data(n_particles, 'n_particles') - call pr % write_data(run_mode, 'run_mode') - call pr % write_data(p % id, 'id') - call pr % write_data(src % wgt, 'weight') - call pr % write_data(src % E, 'energy') - call pr % write_data(src % xyz, 'xyz', length = 3) - call pr % write_data(src % uvw, 'uvw', length = 3) + call write_dataset(file_id, 'filetype', FILETYPE_PARTICLE_RESTART) + call write_dataset(file_id, 'revision', REVISION_PARTICLE_RESTART) + call write_dataset(file_id, 'current_batch', current_batch) + call write_dataset(file_id, 'gen_per_batch', gen_per_batch) + call write_dataset(file_id, 'current_gen', current_gen) + call write_dataset(file_id, 'n_particles', n_particles) + call write_dataset(file_id, 'run_mode', run_mode) + call write_dataset(file_id, 'id', p%id) + call write_dataset(file_id, 'weight', src%wgt) + call write_dataset(file_id, 'energy', src%E) + call write_dataset(file_id, 'xyz', src%xyz) + call write_dataset(file_id, 'uvw', src%uvw) ! Close file - call pr % file_close() + call file_close(file_id) !$omp end critical (WriteParticleRestart) end subroutine write_particle_restart diff --git a/src/source.F90 b/src/source.F90 index 4b0de58b5a..285332ef8d 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -6,17 +6,20 @@ module source use geometry, only: find_cell use geometry_header, only: BASE_UNIVERSE use global + use hdf5_interface, only: file_create, file_open, file_close, read_dataset use math, only: maxwell_spectrum, watt_spectrum use output, only: write_message - use output_interface, only: BinaryOutput use particle_header, only: Particle use random_lcg, only: prn, set_particle_seed, prn_set_stream + use state_point, only: read_source_bank, write_source_bank use string, only: to_str #ifdef MPI use message_passing #endif + use hdf5, only: HID_T + implicit none contains @@ -27,12 +30,12 @@ contains subroutine initialize_source() - character(MAX_FILE_LEN) :: filename integer(8) :: i ! loop index over bank sites integer(8) :: id ! particle id integer(4) :: itmp ! temporary integer - type(Bank), pointer :: src => null() ! source bank site - type(BinaryOutput) :: sp ! statepoint/source binary file + integer(HID_T) :: file_id + character(MAX_FILE_LEN) :: filename + type(Bank), pointer :: src ! source bank site call write_message("Initializing source particles...", 6) @@ -44,10 +47,10 @@ contains &// '...', 6) ! Open the binary file - call sp % file_open(path_source, 'r', serial = .false.) + file_id = file_open(path_source, 'r', parallel=.true.) ! Read the file type - call sp % read_data(itmp, "filetype") + call read_dataset(file_id, "filetype", itmp) ! Check to make sure this is a source file if (itmp /= FILETYPE_SOURCE) then @@ -56,10 +59,10 @@ contains end if ! Read in the source bank - call sp % read_source_bank() + call read_source_bank(file_id) ! Close file - call sp % file_close() + call file_close(file_id) else ! Generation source sites from specified distribution in user input @@ -80,9 +83,9 @@ contains if (write_initial_source) then call write_message('Writing out initial source...', 1) filename = trim(path_output) // 'initial_source.h5' - call sp % file_create(filename, serial = .false.) - call sp % write_source_bank() - call sp % file_close() + file_id = file_create(filename, parallel=.true.) + call write_source_bank(file_id) + call file_close(file_id) end if end subroutine initialize_source @@ -109,28 +112,28 @@ contains integer, save :: num_resamples = 0 ! Number of resamples encountered ! Set weight to one by default - site % wgt = ONE + site%wgt = ONE ! Set the random number generator to the source stream. call prn_set_stream(STREAM_SOURCE) ! Sample position - select case (external_source % type_space) + select case (external_source%type_space) case (SRC_SPACE_BOX) ! Set particle defaults - call p % initialize() + call p%initialize() ! Repeat sampling source location until a good site has been found found = .false. do while (.not.found) ! Coordinates sampled uniformly over a box - p_min = external_source % params_space(1:3) - p_max = external_source % params_space(4:6) + p_min = external_source%params_space(1:3) + p_max = external_source%params_space(4:6) r = (/ (prn(), i = 1,3) /) - site % xyz = p_min + r*(p_max - p_min) + site%xyz = p_min + r*(p_max - p_min) ! Fill p with needed data - p % coord(1) % xyz = site % xyz - p % coord(1) % uvw = [ ONE, ZERO, ZERO ] + p%coord(1)%xyz = site%xyz + p%coord(1)%uvw = [ ONE, ZERO, ZERO ] ! Now search to see if location exists in geometry call find_cell(p, found) @@ -142,24 +145,24 @@ contains end if end if end do - call p % clear() + call p%clear() case (SRC_SPACE_FISSION) ! Repeat sampling source location until a good site has been found found = .false. do while (.not.found) ! Set particle defaults - call p % initialize() + call p%initialize() ! Coordinates sampled uniformly over a box - p_min = external_source % params_space(1:3) - p_max = external_source % params_space(4:6) + p_min = external_source%params_space(1:3) + p_max = external_source%params_space(4:6) r = (/ (prn(), i = 1,3) /) - site % xyz = p_min + r*(p_max - p_min) + site%xyz = p_min + r*(p_max - p_min) ! Fill p with needed data - p % coord(1) % xyz = site % xyz - p % coord(1) % uvw = [ ONE, ZERO, ZERO ] + p%coord(1)%xyz = site%xyz + p%coord(1)%uvw = [ ONE, ZERO, ZERO ] ! Now search to see if location exists in geometry call find_cell(p, found) @@ -171,66 +174,66 @@ contains end if cycle end if - if (p % material == MATERIAL_VOID) then + if (p%material == MATERIAL_VOID) then found = .false. cycle end if - if (.not. materials(p % material) % fissionable) found = .false. + if (.not. materials(p%material)%fissionable) found = .false. end do - call p % clear() + call p%clear() case (SRC_SPACE_POINT) ! Point source - site % xyz = external_source % params_space + site%xyz = external_source%params_space end select ! Sample angle - select case (external_source % type_angle) + select case (external_source%type_angle) case (SRC_ANGLE_ISOTROPIC) ! Sample isotropic distribution phi = TWO*PI*prn() mu = TWO*prn() - ONE - site % uvw(1) = mu - site % uvw(2) = sqrt(ONE - mu*mu) * cos(phi) - site % uvw(3) = sqrt(ONE - mu*mu) * sin(phi) + site%uvw(1) = mu + site%uvw(2) = sqrt(ONE - mu*mu) * cos(phi) + site%uvw(3) = sqrt(ONE - mu*mu) * sin(phi) case (SRC_ANGLE_MONO) ! Monodirectional source - site % uvw = external_source % params_angle + site%uvw = external_source%params_angle case default call fatal_error("No angle distribution specified for external source!") end select ! Sample energy distribution - select case (external_source % type_energy) + select case (external_source%type_energy) case (SRC_ENERGY_MONO) ! Monoenergtic source - site % E = external_source % params_energy(1) - if (site % E >= 20) then + site%E = external_source%params_energy(1) + if (site%E >= 20) then call fatal_error("Source energies above 20 MeV not allowed.") end if case (SRC_ENERGY_MAXWELL) - a = external_source % params_energy(1) + a = external_source%params_energy(1) do ! Sample Maxwellian fission spectrum - site % E = maxwell_spectrum(a) + site%E = maxwell_spectrum(a) ! resample if energy is >= 20 MeV - if (site % E < 20) exit + if (site%E < 20) exit end do case (SRC_ENERGY_WATT) - a = external_source % params_energy(1) - b = external_source % params_energy(2) + a = external_source%params_energy(1) + b = external_source%params_energy(2) do ! Sample Watt fission spectrum - site % E = watt_spectrum(a, b) + site%E = watt_spectrum(a, b) ! resample if energy is >= 20 MeV - if (site % E < 20) exit + if (site%E < 20) exit end do case default diff --git a/src/state_point.F90 b/src/state_point.F90 index 790fd58c57..084289bee0 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -15,9 +15,9 @@ module state_point use constants use error, only: fatal_error, warning use global + use hdf5_interface use output, only: write_message, time_stamp use string, only: to_str, zero_padded, count_digits - use output_interface use tally_header, only: TallyObject use mesh_header, only: StructuredMesh use dict_header, only: ElemKeyValueII, ElemKeyValueCI @@ -26,9 +26,9 @@ module state_point use message_passing #endif - implicit none + use hdf5 - type(BinaryOutput) :: sp ! Statepoint/source output file + implicit none contains @@ -38,23 +38,26 @@ contains subroutine write_state_point() - character(MAX_FILE_LEN) :: filename integer :: i, j, k + integer :: n_order ! loop index for moment orders + integer :: nm_order ! loop index for Ynm moment orders integer, allocatable :: id_array(:) integer, allocatable :: key_array(:) + integer(HID_T) :: file_id + integer(HID_T) :: cmfd_group + integer(HID_T) :: tallies_group, tally_group + integer(HID_T) :: meshes_group, mesh_group + integer(HID_T) :: filter_group, moments_group + character(8) :: moment_name ! name of moment (e.g, P3) + character(MAX_FILE_LEN) :: filename type(StructuredMesh), pointer :: mesh type(TallyObject), pointer :: tally type(ElemKeyValueII), pointer :: current type(ElemKeyValueII), pointer :: next - character(8) :: moment_name ! name of moment (e.g, P3) - integer :: n_order ! loop index for moment orders - integer :: nm_order ! loop index for Ynm moment orders ! Set filename for state point filename = trim(path_output) // 'statepoint.' // & & zero_padded(current_batch, count_digits(n_max_batches)) - - ! Append appropriate extension filename = trim(filename) // '.h5' ! Write message @@ -62,145 +65,126 @@ contains if (master) then ! Create statepoint file - call sp % file_create(filename) + file_id = file_create(filename) ! Write file type - call sp % write_data(FILETYPE_STATEPOINT, "filetype") + call write_dataset(file_id, "filetype", FILETYPE_STATEPOINT) ! Write revision number for state point file - call sp % write_data(REVISION_STATEPOINT, "revision") + call write_dataset(file_id, "revision", REVISION_STATEPOINT) ! Write OpenMC version - call sp % write_data(VERSION_MAJOR, "version_major") - call sp % write_data(VERSION_MINOR, "version_minor") - call sp % write_data(VERSION_RELEASE, "version_release") + call write_dataset(file_id, "version_major", VERSION_MAJOR) + call write_dataset(file_id, "version_minor", VERSION_MINOR) + call write_dataset(file_id, "version_release", VERSION_RELEASE) ! Write current date and time - call sp % write_data(time_stamp(), "date_and_time") + call write_dataset(file_id, "date_and_time", time_stamp()) ! Write path to input - call sp % write_data(path_input, "path") + call write_dataset(file_id, "path", path_input) ! Write out random number seed - call sp % write_data(seed, "seed") + call write_dataset(file_id, "seed", seed) ! Write run information - call sp % write_data(run_mode, "run_mode") - call sp % write_data(n_particles, "n_particles") - call sp % write_data(n_batches, "n_batches") + call write_dataset(file_id, "run_mode", run_mode) + call write_dataset(file_id, "n_particles", n_particles) + call write_dataset(file_id, "n_batches", n_batches) ! Write out current batch number - call sp % write_data(current_batch, "current_batch") + call write_dataset(file_id, "current_batch", current_batch) ! Indicate whether source bank is stored in statepoint if (source_separate) then - call sp % write_data(0, "source_present") + call write_dataset(file_id, "source_present", 0) else - call sp % write_data(1, "source_present") + call write_dataset(file_id, "source_present", 1) end if ! Write out information for eigenvalue run if (run_mode == MODE_EIGENVALUE) then - call sp % write_data(n_inactive, "n_inactive") - call sp % write_data(gen_per_batch, "gen_per_batch") - call sp % write_data(k_generation, "k_generation", & - length=current_batch*gen_per_batch) - call sp % write_data(entropy, "entropy", & - length=current_batch*gen_per_batch) - call sp % write_data(k_col_abs, "k_col_abs") - call sp % write_data(k_col_tra, "k_col_tra") - call sp % write_data(k_abs_tra, "k_abs_tra") - call sp % write_data(k_combined, "k_combined", length=2) + call write_dataset(file_id, "n_inactive", n_inactive) + call write_dataset(file_id, "gen_per_batch", gen_per_batch) + call write_dataset(file_id, "k_generation", k_generation) + call write_dataset(file_id, "entropy", entropy) + call write_dataset(file_id, "k_col_abs", k_col_abs) + call write_dataset(file_id, "k_col_tra", k_col_tra) + call write_dataset(file_id, "k_abs_tra", k_abs_tra) + call write_dataset(file_id, "k_combined", k_combined) ! Write out CMFD info if (cmfd_on) then - call sp % open_group("cmfd") - call sp % close_group() - call sp % write_data(1, "cmfd_on") - call sp % write_data(cmfd % indices, "indices", length=4, group="cmfd") - call sp % write_data(cmfd % k_cmfd, "k_cmfd", length=current_batch, & - group="cmfd") - call sp % write_data(cmfd % cmfd_src, "cmfd_src", & - length=(/cmfd % indices(4), cmfd % indices(1), & - cmfd % indices(2), cmfd % indices(3)/), & - group="cmfd") - call sp % write_data(cmfd % entropy, "cmfd_entropy", & - length=current_batch, group="cmfd") - call sp % write_data(cmfd % balance, "cmfd_balance", & - length=current_batch, group="cmfd") - call sp % write_data(cmfd % dom, "cmfd_dominance", & - length = current_batch, group="cmfd") - call sp % write_data(cmfd % src_cmp, "cmfd_srccmp", & - length = current_batch, group="cmfd") + call write_dataset(file_id, "cmfd_on", 1) + + cmfd_group = create_group(file_id, "cmfd") + call write_dataset(cmfd_group, "indices", cmfd%indices) + call write_dataset(cmfd_group, "k_cmfd", cmfd%k_cmfd) + call write_dataset(cmfd_group, "cmfd_src", cmfd%cmfd_src) + call write_dataset(cmfd_group, "cmfd_entropy", cmfd%entropy) + call write_dataset(cmfd_group, "cmfd_balance", cmfd%balance) + call write_dataset(cmfd_group, "cmfd_dominance", cmfd%dom) + call write_dataset(cmfd_group, "cmfd_srccmp", cmfd%src_cmp) + call close_group(cmfd_group) else - call sp % write_data(0, "cmfd_on") + call write_dataset(file_id, "cmfd_on", 0) end if end if - call sp % open_group("tallies") - call sp % close_group() + tallies_group = create_group(file_id, "tallies") ! Write number of meshes - call sp % write_data(n_meshes, "n_meshes", group="tallies/meshes") + meshes_group = create_group(tallies_group, "meshes") + call write_dataset(meshes_group, "n_meshes", n_meshes) if (n_meshes > 0) then ! Print list of mesh IDs - current => mesh_dict % keys() + current => mesh_dict%keys() allocate(id_array(n_meshes)) allocate(key_array(n_meshes)) i = 1 do while (associated(current)) - key_array(i) = current % key - id_array(i) = current % value + key_array(i) = current%key + id_array(i) = current%value ! Move to next mesh - next => current % next + next => current%next deallocate(current) current => next i = i + 1 end do - call sp % write_data(id_array, "ids", & - group="tallies/meshes", length=n_meshes) - call sp % write_data(key_array, "keys", & - group="tallies/meshes", length=n_meshes) + call write_dataset(meshes_group, "ids", id_array) + call write_dataset(meshes_group, "keys", key_array) deallocate(key_array) ! Write information for meshes MESH_LOOP: do i = 1, n_meshes - mesh => meshes(id_array(i)) + mesh_group = create_group(meshes_group, "mesh " // trim(to_str(mesh%id))) - call sp % write_data(mesh % id, "id", & - group="tallies/meshes/mesh " // trim(to_str(mesh % id))) - call sp % write_data(mesh % type, "type", & - group="tallies/meshes/mesh " // trim(to_str(mesh % id))) - call sp % write_data(mesh % n_dimension, "n_dimension", & - group="tallies/meshes/mesh " // trim(to_str(mesh % id))) - call sp % write_data(mesh % dimension, "dimension", & - group="tallies/meshes/mesh " // trim(to_str(mesh % id)), & - length=mesh % n_dimension) - call sp % write_data(mesh % lower_left, "lower_left", & - group="tallies/meshes/mesh " // trim(to_str(mesh % id)), & - length=mesh % n_dimension) - call sp % write_data(mesh % upper_right, "upper_right", & - group="tallies/meshes/mesh " // trim(to_str(mesh % id)), & - length=mesh % n_dimension) - call sp % write_data(mesh % width, "width", & - group="tallies/meshes/mesh " // trim(to_str(mesh % id)), & - length=mesh % n_dimension) + call write_dataset(mesh_group, "id", mesh%id) + call write_dataset(mesh_group, "type", mesh%type) + call write_dataset(mesh_group, "n_dimension", mesh%n_dimension) + call write_dataset(mesh_group, "dimension", mesh%dimension) + call write_dataset(mesh_group, "lower_left", mesh%lower_left) + call write_dataset(mesh_group, "upper_right", mesh%upper_right) + call write_dataset(mesh_group, "width", mesh%width) + + call close_group(mesh_group) end do MESH_LOOP deallocate(id_array) - end if + call close_group(meshes_group) + ! Write number of tallies - call sp % write_data(n_tallies, "n_tallies", group="tallies") + call write_dataset(tallies_group, "n_tallies", n_tallies) if (n_tallies > 0) then @@ -211,14 +195,12 @@ contains ! Write all tally information except results do i = 1, n_tallies tally => tallies(i) - key_array(i) = tally % id + key_array(i) = tally%id id_array(i) = i end do - call sp % write_data(id_array, "ids", & - group="tallies", length=n_tallies) - call sp % write_data(key_array, "keys", & - group="tallies", length=n_tallies) + call write_dataset(tallies_group, "ids", id_array) + call write_dataset(tallies_group, "keys", key_array) deallocate(key_array) @@ -227,110 +209,92 @@ contains ! Get pointer to tally tally => tallies(i) + tally_group = create_group(tallies_group, "tally " // & + trim(to_str(tally%id))) - call sp % write_data(tally % estimator, "estimator", & - group="tallies/tally " // trim(to_str(tally % id))) - call sp % write_data(tally % n_realizations, "n_realizations", & - group="tallies/tally " // trim(to_str(tally % id))) - call sp % write_data(tally % n_filters, "n_filters", & - group="tallies/tally " // trim(to_str(tally % id))) + call write_dataset(tally_group, "estimator", tally%estimator) + call write_dataset(tally_group, "n_realizations", tally%n_realizations) + call write_dataset(tally_group, "n_filters", tally%n_filters) ! Write filter information - FILTER_LOOP: do j = 1, tally % n_filters + FILTER_LOOP: do j = 1, tally%n_filters + filter_group = create_group(tally_group, "filter " // & + trim(to_str(j))) - call sp % write_data(tally % filters(j) % type, "type", & - group="tallies/tally " // trim(to_str(tally % id)) // & - "/filter " // to_str(j)) - call sp % write_data(tally % filters(j) % offset, "offset", & - group="tallies/tally " // trim(to_str(tally % id)) // & - "/filter " // to_str(j)) - call sp % write_data(tally % filters(j) % n_bins, "n_bins", & - group="tallies/tally " // trim(to_str(tally % id)) // & - "/filter " // to_str(j)) - if (tally % filters(j) % type == FILTER_ENERGYIN .or. & - tally % filters(j) % type == FILTER_ENERGYOUT) then - call sp % write_data(tally % filters(j) % real_bins, "bins", & - group="tallies/tally " // trim(to_str(tally % id)) // & - "/filter " // to_str(j), & - length=size(tally % filters(j) % real_bins)) + call write_dataset(filter_group, "type", tally%filters(j)%type) + call write_dataset(filter_group, "offset", tally%filters(j)%offset) + call write_dataset(filter_group, "n_bins", tally%filters(j)%n_bins) + if (tally%filters(j)%type == FILTER_ENERGYIN .or. & + tally%filters(j)%type == FILTER_ENERGYOUT) then + call write_dataset(filter_group, "bins", & + tally%filters(j)%real_bins) else - call sp % write_data(tally % filters(j) % int_bins, "bins", & - group="tallies/tally " // trim(to_str(tally % id)) // & - "/filter " // to_str(j), & - length=size(tally % filters(j) % int_bins)) + call write_dataset(filter_group, "bins", & + tally%filters(j)%int_bins) end if + call close_group(filter_group) end do FILTER_LOOP - call sp % write_data(tally % n_nuclide_bins, "n_nuclides", & - group="tallies/tally " // trim(to_str(tally % id))) + call write_dataset(tally_group, "n_nuclides", tally%n_nuclide_bins) ! Set up nuclide bin array and then write - allocate(key_array(tally % n_nuclide_bins)) - NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins - if (tally % nuclide_bins(j) > 0) then - key_array(j) = nuclides(tally % nuclide_bins(j)) % zaid + allocate(key_array(tally%n_nuclide_bins)) + NUCLIDE_LOOP: do j = 1, tally%n_nuclide_bins + if (tally%nuclide_bins(j) > 0) then + key_array(j) = nuclides(tally%nuclide_bins(j))%zaid else - key_array(j) = tally % nuclide_bins(j) + key_array(j) = tally%nuclide_bins(j) end if end do NUCLIDE_LOOP - call sp % write_data(key_array, "nuclides", & - group="tallies/tally " // trim(to_str(tally % id)), & - length=tally % n_nuclide_bins) + call write_dataset(tally_group, "nuclides", key_array) deallocate(key_array) - call sp % write_data(tally % n_score_bins, "n_score_bins", & - group="tallies/tally " // trim(to_str(tally % id))) - call sp % write_data(tally % score_bins, "score_bins", & - group="tallies/tally " // trim(to_str(tally % id)), & - length=tally % n_score_bins) - call sp % write_data(tally % n_user_score_bins, "n_user_score_bins", & - group="tallies/tally " // to_str(tally % id)) + call write_dataset(tally_group, "n_score_bins", tally%n_score_bins) + call write_dataset(tally_group, "score_bins", tally%score_bins) + call write_dataset(tally_group, "n_user_score_bins", tally%n_user_score_bins) ! Write explicit moment order strings for each score bin + moments_group = create_group(tally_group, "moments") k = 1 - MOMENT_LOOP: do j = 1, tally % n_user_score_bins - select case(tally % score_bins(k)) + MOMENT_LOOP: do j = 1, tally%n_user_score_bins + select case(tally%score_bins(k)) case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) - moment_name = 'P' // trim(to_str(tally % moment_order(k))) - call sp % write_data(moment_name, "order" // trim(to_str(k)), & - group="tallies/tally " // trim(to_str(tally % id)) // & - "/moments") + moment_name = 'P' // trim(to_str(tally%moment_order(k))) + call write_dataset(moments_group, "order" // trim(to_str(k)), moment_name) k = k + 1 case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN) - do n_order = 0, tally % moment_order(k) + do n_order = 0, tally%moment_order(k) moment_name = 'P' // trim(to_str(n_order)) - call sp % write_data(moment_name, "order" // trim(to_str(k)), & - group="tallies/tally " // trim(to_str(tally % id)) // & - "/moments") + call write_dataset(moments_group, "order" // trim(to_str(k)), moment_name) k = k + 1 end do case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, & SCORE_TOTAL_YN) - do n_order = 0, tally % moment_order(k) + do n_order = 0, tally%moment_order(k) do nm_order = -n_order, n_order moment_name = 'Y' // trim(to_str(n_order)) // ',' // & trim(to_str(nm_order)) - call sp % write_data(moment_name, "order" // & - trim(to_str(k)), & - group="tallies/tally " // trim(to_str(tally % id)) // & - "/moments") + call write_dataset(moments_group, "order" // & + trim(to_str(k)), moment_name) k = k + 1 end do end do case default moment_name = '' - call sp % write_data(moment_name, "order" // trim(to_str(k)), & - group="tallies/tally " // trim(to_str(tally % id)) // & - "/moments") + call write_dataset(moments_group, "order" // trim(to_str(k)), & + moment_name) k = k + 1 end select - end do MOMENT_LOOP + call close_group(moments_group) + call close_group(tally_group) end do TALLY_METADATA end if + + call close_group(tallies_group) end if ! Check for the no-tally-reduction method @@ -338,47 +302,41 @@ contains ! If using the no-tally-reduction method, we need to collect tally ! results before writing them to the state point file. - call write_tally_results_nr() + call write_tally_results_nr(file_id) elseif (master) then ! Write number of global realizations - call sp % write_data(n_realizations, "n_realizations") + call write_dataset(file_id, "n_realizations", n_realizations) ! Write global tallies - call sp % write_data(N_GLOBAL_TALLIES, "n_global_tallies") - call sp % write_tally_result(global_tallies, "global_tallies", & - n1=N_GLOBAL_TALLIES, n2=1) + call write_dataset(file_id, "n_global_tallies", N_GLOBAL_TALLIES) + call write_dataset(file_id, "global_tallies", global_tallies) ! Write tallies + tallies_group = open_group(file_id, "tallies") if (tallies_on) then - ! Indicate that tallies are on - call sp % write_data(1, "tallies_present", group="tallies") + call write_dataset(tallies_group, "tallies_present", 1) ! Write all tally results TALLY_RESULTS: do i = 1, n_tallies - ! Set point to current tally tally => tallies(i) ! Write sum and sum_sq for each bin - call sp % write_tally_result(tally % results, "results", & - group="tallies/tally " // trim(to_str(tally % id)), & - n1=size(tally % results, 1), n2=size(tally % results, 2)) - + tally_group = open_group(tallies_group, "tally " // to_str(tally%id)) + call write_dataset(tally_group, "results", tally%results) + call close_group(tally_group) end do TALLY_RESULTS else - ! Indicate tallies are off - call sp % write_data(0, "tallies_present", group="tallies") - + call write_dataset(tallies_group, "tallies_present", 0) end if - ! Close the file for serial writing - call sp % file_close() - + call close_group(tallies_group) + call file_close(file_id) end if if (master .and. n_tallies > 0) then @@ -393,71 +351,64 @@ contains subroutine write_source_point() - type(BinaryOutput) :: sp + integer(HID_T) :: file_id character(MAX_FILE_LEN) :: filename ! Check to write out source for a specified batch - if (sourcepoint_batch % contains(current_batch)) then + if (sourcepoint_batch%contains(current_batch)) then ! Create or open up file if (source_separate) then - ! Set filename filename = trim(path_output) // 'source.' // & & zero_padded(current_batch, count_digits(n_max_batches)) - filename = trim(filename) // '.h5' ! Write message for new file creation call write_message("Creating source file " // trim(filename) // "...", & - &1) + 1) ! Create separate source file - call sp % file_create(filename, serial = .false.) + file_id = file_create(filename, parallel=.true.) ! Write file type - call sp % write_data(FILETYPE_SOURCE, "filetype") - + call write_dataset(file_id, "filetype", FILETYPE_SOURCE) else - ! Set filename for state point filename = trim(path_output) // 'statepoint.' // & zero_padded(current_batch, count_digits(n_max_batches)) filename = trim(filename) // '.h5' ! Reopen statepoint file in parallel - call sp % file_open(filename, 'w', serial = .false.) - + file_id = file_open(filename, 'w', parallel=.true.) end if ! Write out source - call sp % write_source_bank() + call write_source_bank(file_id) ! Close file - call sp % file_close() - + call file_close(file_id) end if ! Also check to write source separately in overwritten file if (source_latest) then ! Set filename - filename = trim(path_output) // 'source' - filename = trim(filename) // '.h5' + filename = trim(path_output) // 'source' // '.h5' ! Write message for new file creation call write_message("Creating source file " // trim(filename) // "...", 1) ! Always create this file because it will be overwritten - call sp % file_create(filename, serial = .false.) + file_id = file_create(filename, parallel=.true.) ! Write file type - call sp % write_data(FILETYPE_SOURCE, "filetype") + call write_dataset(file_id, "filetype", FILETYPE_SOURCE) ! Write out source - call sp % write_source_bank() + call write_source_bank(file_id) ! Close file - call sp % file_close() + call file_close(file_id) end if @@ -467,12 +418,14 @@ contains ! WRITE_TALLY_RESULTS_NR !=============================================================================== - subroutine write_tally_results_nr() + subroutine write_tally_results_nr(file_id) + integer(HID_T), intent(in) :: file_id integer :: i ! loop index integer :: n ! number of filter bins integer :: m ! number of score bins 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) #ifdef MPI @@ -489,16 +442,18 @@ contains if (master) then ! Write number of realizations - call sp % write_data(n_realizations, "n_realizations") + call write_dataset(file_id, "n_realizations", n_realizations) ! Write number of global tallies - call sp % write_data(N_GLOBAL_TALLIES, "n_global_tallies") + call write_dataset(file_id, "n_global_tallies", N_GLOBAL_TALLIES) + + tallies_group = open_group(file_id, "tallies") 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 + 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 @@ -510,19 +465,17 @@ 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(:)%sum = global_temp(1,:) + global_tallies(:)%sum_sq = global_temp(2,:) 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,:) - + tallyresult_temp(:,1)%sum = global_temp(1,:) + tallyresult_temp(:,1)%sum_sq = global_temp(2,:) ! Write out global tallies sum and sum_sq - call sp % write_tally_result(tallyresult_temp, "global_tallies", & - n1=N_GLOBAL_TALLIES, n2=1) + call write_dataset(file_id, "global_tallies", tallyresult_temp) ! Deallocate temporary tally result deallocate(tallyresult_temp) @@ -537,17 +490,17 @@ contains if (tallies_on) then ! Indicate that tallies are on if (master) then - call sp % write_data(1, "tallies_present", group="tallies") + call write_dataset(tallies_group, "tallies_present", 1) ! Build list of tally IDs - current => tally_dict % keys() + current => tally_dict%keys() allocate(id_array(n_tallies)) i = 1 do while (associated(current)) - id_array(i) = current % value + id_array(i) = current%value ! Move to next tally - next => current % next + next => current%next deallocate(current) current => next i = i + 1 @@ -561,17 +514,20 @@ contains tally => tallies(i) ! Determine size of tally results array - m = size(tally % results, 1) - n = size(tally % results, 2) + m = size(tally%results, 1) + n = size(tally%results, 2) 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(:,:)%sum + tally_temp(2,:,:) = tally%results(:,:)%sum_sq if (master) then + tally_group = open_group(tallies_group, "tally " // & + trim(to_str(tally%id))) + ! The MPI_IN_PLACE specifier allows the master to copy values into ! a receive buffer without having a temporary variable #ifdef MPI @@ -582,18 +538,17 @@ 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(:,:)%sum = tally_temp(1,:,:) + tally%results(:,:)%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,:,:) + tallyresult_temp(:,:)%sum = tally_temp(1,:,:) + tallyresult_temp(:,:)%sum_sq = tally_temp(2,:,:) ! Write reduced tally results to file - call sp % write_tally_result(tally % results, "results", & - group="tallies/tally " // trim(to_str(tally % id)), n1=m, n2=n) + call write_dataset(tally_group, "results", tally%results) ! Deallocate temporary tally result deallocate(tallyresult_temp) @@ -607,6 +562,8 @@ contains ! Deallocate temporary copy of tally results deallocate(tally_temp) + + if (master) call close_group(tally_group) end do TALLY_RESULTS deallocate(id_array) @@ -614,10 +571,12 @@ contains else if (master) then ! Indicate that tallies are off - call sp % write_data(0, "tallies_present", group="tallies") + call write_dataset(tallies_group, "tallies_present", 0) end if end if + if (master) call close_group(tallies_group) + end subroutine write_tally_results_nr !=============================================================================== @@ -626,45 +585,49 @@ contains subroutine load_state_point() - character(MAX_FILE_LEN) :: path_temp - character(19) :: current_time integer :: i, j, k - integer :: length(4) integer :: int_array(3) - integer, allocatable :: id_array(:) - integer, allocatable :: key_array(:) integer :: curr_key - integer, allocatable :: temp_array(:) - logical :: source_present - real(8) :: real_array(3) - type(StructuredMesh), pointer :: mesh - type(TallyObject), pointer :: tally integer :: n_order ! loop index for moment orders integer :: nm_order ! loop index for Ynm moment orders + integer, allocatable :: id_array(:) + integer, allocatable :: key_array(:) + integer, allocatable :: temp_array(:) + integer(HID_T) :: file_id + integer(HID_T) :: cmfd_group + integer(HID_T) :: tallies_group, tally_group + integer(HID_T) :: meshes_group, mesh_group + integer(HID_T) :: filter_group, moments_group + real(8) :: real_array(3) + logical :: source_present + character(MAX_FILE_LEN) :: path_temp + character(19) :: current_time character(8) :: moment_name ! name of moment (e.g, P3, Y-1,1) + type(StructuredMesh), pointer :: mesh + type(TallyObject), pointer :: tally ! Write message call write_message("Loading state point " // trim(path_state_point) & - &// "...", 1) + // "...", 1) ! Open file for reading - call sp % file_open(path_state_point, 'r', serial = .false.) + file_id = file_open(path_state_point, 'r', parallel=.true.) ! Read filetype - call sp % read_data(int_array(1), "filetype") + call read_dataset(file_id, "filetype", int_array(1)) ! Read revision number for state point file and make sure it matches with ! current version - call sp % read_data(int_array(1), "revision") + call read_dataset(file_id, "revision", int_array(1)) if (int_array(1) /= REVISION_STATEPOINT) then call fatal_error("State point version does not match current version & &in OpenMC.") end if ! Read OpenMC version - call sp % read_data(int_array(1), "version_major") - call sp % read_data(int_array(2), "version_minor") - call sp % read_data(int_array(3), "version_release") + call read_dataset(file_id, "version_major", int_array(1)) + call read_dataset(file_id, "version_minor", int_array(2)) + call read_dataset(file_id, "version_release", int_array(3)) if (int_array(1) /= VERSION_MAJOR .or. int_array(2) /= VERSION_MINOR & .or. int_array(3) /= VERSION_RELEASE) then if (master) call warning("State point file was created with a different & @@ -672,27 +635,27 @@ contains end if ! Read date and time - call sp % read_data(current_time, "date_and_time") + call read_dataset(file_id, "date_and_time", current_time) ! Read path to input - call sp % read_data(path_temp, "path") + call read_dataset(file_id, "path", path_temp) ! Read and overwrite random number seed - call sp % read_data(seed, "seed") + call read_dataset(file_id, "seed", seed) ! Read and overwrite run information except number of batches - call sp % read_data(run_mode, "run_mode") - call sp % read_data(n_particles, "n_particles") - call sp % read_data(int_array(1), "n_batches") + call read_dataset(file_id, "run_mode", run_mode) + call read_dataset(file_id, "n_particles", n_particles) + call read_dataset(file_id, "n_batches", int_array(1)) ! Take maximum of statepoint n_batches and input n_batches n_batches = max(n_batches, int_array(1)) ! Read batch number to restart at - call sp % read_data(restart_batch, "current_batch") + call read_dataset(file_id, "current_batch", restart_batch) ! Check for source in statepoint if needed - call sp % read_data(int_array(1), "source_present") + call read_dataset(file_id, "source_present", int_array(1)) if (int_array(1) == 1) then source_present = .true. else @@ -706,43 +669,45 @@ contains ! Read information specific to eigenvalue run if (run_mode == MODE_EIGENVALUE) then - call sp % read_data(int_array(1), "n_inactive") - call sp % read_data(gen_per_batch, "gen_per_batch") - call sp % read_data(k_generation, "k_generation", & - length=restart_batch*gen_per_batch) - call sp % read_data(entropy, "entropy", length=restart_batch*gen_per_batch) - call sp % read_data(k_col_abs, "k_col_abs") - call sp % read_data(k_col_tra, "k_col_tra") - call sp % read_data(k_abs_tra, "k_abs_tra") - call sp % read_data(real_array(1:2), "k_combined", length=2) + call read_dataset(file_id, "n_inactive", int_array(1)) + call read_dataset(file_id, "gen_per_batch", gen_per_batch) + call read_dataset(file_id, "k_generation", & + k_generation(1:restart_batch*gen_per_batch)) + call read_dataset(file_id, "entropy", & + entropy(1:restart_batch*gen_per_batch)) + call read_dataset(file_id, "k_col_abs", k_col_abs) + call read_dataset(file_id, "k_col_tra", k_col_tra) + call read_dataset(file_id, "k_abs_tra", k_abs_tra) + call read_dataset(file_id, "k_combined", real_array(1:2)) ! Take maximum of statepoint n_inactive and input n_inactive n_inactive = max(n_inactive, int_array(1)) ! Read in to see if CMFD was on - call sp % read_data(int_array(1), "cmfd_on") + call read_dataset(file_id, "cmfd_on", int_array(1)) ! Read in CMFD info if (int_array(1) == 1) then - call sp % read_data(cmfd % indices, "indices", length=4, group="cmfd") - call sp % read_data(cmfd % k_cmfd, "k_cmfd", length=restart_batch, & - group="cmfd") - length = cmfd % indices([4,1,2,3]) - call sp % read_data(cmfd % cmfd_src, "cmfd_src", & - length=length, group="cmfd") - call sp % read_data(cmfd % entropy, "cmfd_entropy", & - length=restart_batch, group="cmfd") - call sp % read_data(cmfd % balance, "cmfd_balance", & - length=restart_batch, group="cmfd") - call sp % read_data(cmfd % dom, "cmfd_dominance", & - length = restart_batch, group="cmfd") - call sp % read_data(cmfd % src_cmp, "cmfd_srccmp", & - length = restart_batch, group="cmfd") + cmfd_group = open_group(file_id, "cmfd") + call read_dataset(cmfd_group, "indices", cmfd%indices) + call read_dataset(cmfd_group, "k_cmfd", cmfd%k_cmfd(1:restart_batch)) + call read_dataset(cmfd_group, "cmfd_src", cmfd%cmfd_src) + call read_dataset(cmfd_group, "cmfd_entropy", & + cmfd%entropy(1:restart_batch)) + call read_dataset(cmfd_group, "cmfd_balance", & + cmfd%balance(1:restart_batch)) + call read_dataset(cmfd_group, "cmfd_dominance", & + cmfd%dom(1:restart_batch)) + call read_dataset(cmfd_group, "cmfd_srccmp", & + cmfd%src_cmp(1:restart_batch)) + call close_group(cmfd_group) end if end if ! Read number of meshes - call sp % read_data(n_meshes, "n_meshes", group="tallies/meshes") + tallies_group = open_group(file_id, "tallies") + meshes_group = open_group(tallies_group, "meshes") + call read_dataset(meshes_group, "n_meshes", n_meshes) if (n_meshes > 0) then @@ -750,10 +715,8 @@ contains allocate(id_array(n_meshes)) allocate(key_array(n_meshes)) - call sp % read_data(id_array, "ids", & - group="tallies/meshes", length=n_meshes) - call sp % read_data(key_array, "keys", & - group="tallies/meshes", length=n_meshes) + call read_dataset(meshes_group, "ids", id_array) + call read_dataset(meshes_group, "keys", key_array) ! Read and overwrite mesh information MESH_LOOP: do i = 1, n_meshes @@ -761,25 +724,16 @@ contains mesh => meshes(id_array(i)) curr_key = key_array(id_array(i)) - call sp % read_data(mesh % id, "id", & - group="tallies/meshes/mesh " // trim(to_str(curr_key))) - call sp % read_data(mesh % type, "type", & - group="tallies/meshes/mesh " // trim(to_str(curr_key))) - call sp % read_data(mesh % n_dimension, "n_dimension", & - group="tallies/meshes/mesh " // trim(to_str(meshes(i) % id))) - call sp % read_data(mesh % dimension, "dimension", & - group="tallies/meshes/mesh " // trim(to_str(curr_key)), & - length=mesh % n_dimension) - call sp % read_data(mesh % lower_left, "lower_left", & - group="tallies/meshes/mesh " // trim(to_str(curr_key)), & - length=mesh % n_dimension) - call sp % read_data(mesh % upper_right, "upper_right", & - group="tallies/meshes/mesh " // trim(to_str(curr_key)), & - length=mesh % n_dimension) - call sp % read_data(mesh % width, "width", & - group="tallies/meshes/mesh " // trim(to_str(curr_key)), & - length=meshes(i) % n_dimension) - + mesh_group = open_group(meshes_group, "mesh " // & + trim(to_str(curr_key))) + call read_dataset(mesh_group, "id", mesh%id) + call read_dataset(mesh_group, "type", mesh%type) + call read_dataset(mesh_group, "n_dimension", mesh%n_dimension) + call read_dataset(mesh_group, "dimension", mesh%dimension) + call read_dataset(mesh_group, "lower_left", mesh%lower_left) + call read_dataset(mesh_group, "upper_right", mesh%upper_right) + call read_dataset(mesh_group, "width", mesh%width) + call close_group(mesh_group) end do MESH_LOOP deallocate(id_array) @@ -787,15 +741,17 @@ contains end if + call close_group(meshes_group) + ! Read and overwrite number of tallies - call sp % read_data(n_tallies, "n_tallies", group="tallies") + call read_dataset(tallies_group, "n_tallies", n_tallies) ! Read list of tally keys-> IDs allocate(id_array(n_tallies)) allocate(key_array(n_tallies)) - call sp % read_data(id_array, "ids", group="tallies", length=n_tallies) - call sp % read_data(key_array, "keys", group="tallies", length=n_tallies) + call read_dataset(tallies_group, "ids", id_array) + call read_dataset(tallies_group, "keys", key_array) ! Read in tally metadata TALLY_METADATA: do i = 1, n_tallies @@ -803,99 +759,85 @@ contains ! Get pointer to tally tally => tallies(i) curr_key = key_array(id_array(i)) + tally_group = open_group(tallies_group, "tally " // & + trim(to_str(curr_key))) - call sp % read_data(tally % estimator, "estimator", & - group="tallies/tally " // trim(to_str(curr_key))) - call sp % read_data(tally % n_realizations, "n_realizations", & - group="tallies/tally " // trim(to_str(curr_key))) - call sp % read_data(tally % n_filters, "n_filters", & - group="tallies/tally " // trim(to_str(curr_key))) + call read_dataset(tally_group, "estimator", tally%estimator) + call read_dataset(tally_group, "n_realizations", tally%n_realizations) + call read_dataset(tally_group, "n_filters", tally%n_filters) - FILTER_LOOP: do j = 1, tally % n_filters - call sp % read_data(tally % filters(j) % type, "type", & - group="tallies/tally " // trim(to_str(curr_key)) // & - "/filter " // to_str(j)) - call sp % read_data(tally % filters(j) % offset, "offset", & - group="tallies/tally " // trim(to_str(curr_key)) // & - "/filter " // to_str(j)) - call sp % read_data(tally % filters(j) % n_bins, "n_bins", & - group="tallies/tally " // trim(to_str(curr_key)) // & - "/filter " // to_str(j)) - if (tally % filters(j) % type == FILTER_ENERGYIN .or. & - tally % filters(j) % type == FILTER_ENERGYOUT) then - call sp % read_data(tally % filters(j) % real_bins, "bins", & - group="tallies/tally " // trim(to_str(curr_key)) // & - "/filter " // to_str(j), & - length=size(tally % filters(j) % real_bins)) + FILTER_LOOP: do j = 1, tally%n_filters + filter_group = open_group(tally_group, "filter " // trim(to_str(j))) + + call read_dataset(filter_group, "type", tally%filters(j)%type) + call read_dataset(filter_group, "offset", tally%filters(j)%offset) + call read_dataset(filter_group, "n_bins", tally%filters(j)%n_bins) + if (tally%filters(j)%type == FILTER_ENERGYIN .or. & + tally%filters(j)%type == FILTER_ENERGYOUT) then + call read_dataset(filter_group, "bins", tally%filters(j)%real_bins) else - call sp % read_data(tally % filters(j) % int_bins, "bins", & - group="tallies/tally " // trim(to_str(curr_key)) // & - "/filter " // to_str(j), & - length=size(tally % filters(j) % int_bins)) + call read_dataset(filter_group, "bins", tally%filters(j)%int_bins) end if + call close_group(filter_group) end do FILTER_LOOP - call sp % read_data(tally % n_nuclide_bins, "n_nuclides", & - group="tallies/tally " // trim(to_str(curr_key))) + call read_dataset(tally_group, "n_nuclides", tally%n_nuclide_bins) ! Set up nuclide bin array and then read - allocate(temp_array(tally % n_nuclide_bins)) - call sp % read_data(temp_array, "nuclides", & - group="tallies/tally " // trim(to_str(curr_key)), & - length=tally % n_nuclide_bins) + allocate(temp_array(tally%n_nuclide_bins)) + call read_dataset(tally_group, "nuclides", temp_array) - NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins + NUCLIDE_LOOP: do j = 1, tally%n_nuclide_bins if (temp_array(j) > 0) then - tally % nuclide_bins(j) = temp_array(j) + tally%nuclide_bins(j) = temp_array(j) else - tally % nuclide_bins(j) = temp_array(j) + tally%nuclide_bins(j) = temp_array(j) end if end do NUCLIDE_LOOP deallocate(temp_array) ! Write number of score bins, score bins, user score bins - call sp % read_data(tally % n_score_bins, "n_score_bins", & - group="tallies/tally " // trim(to_str(curr_key))) - call sp % read_data(tally % score_bins, "score_bins", & - group="tallies/tally " // trim(to_str(curr_key)), & - length=tally % n_score_bins) - call sp % read_data(tally % n_user_score_bins, "n_user_score_bins", & - group="tallies/tally " // trim(to_str(curr_key))) + call read_dataset(tally_group, "n_score_bins", tally%n_score_bins) + call read_dataset(tally_group, "score_bins", tally%score_bins) + call read_dataset(tally_group, "n_user_score_bins", tally%n_user_score_bins) ! Read explicit moment order strings for each score bin k = 1 - MOMENT_LOOP: do j = 1, tally % n_user_score_bins - select case(tally % score_bins(k)) + moments_group = open_group(tally_group, "moments") + MOMENT_LOOP: do j = 1, tally%n_user_score_bins + select case(tally%score_bins(k)) case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) - call sp % read_data(moment_name, "order" // trim(to_str(k)), & - group="tallies/tally " // trim(to_str(curr_key)) // "/moments") + call read_dataset(moments_group, "order" // trim(to_str(k)), & + moment_name) k = k + 1 case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN) - do n_order = 0, tally % moment_order(k) - call sp % read_data(moment_name, "order" // trim(to_str(k)), & - group="tallies/tally " // trim(to_str(curr_key)) // "/moments") + do n_order = 0, tally%moment_order(k) + call read_dataset(moments_group, "order" // trim(to_str(k)), & + moment_name) k = k + 1 end do case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, & SCORE_TOTAL_YN) - do n_order = 0, tally % moment_order(k) + do n_order = 0, tally%moment_order(k) do nm_order = -n_order, n_order - call sp % read_data(moment_name, "order" // trim(to_str(k)), & - group="tallies/tally " // trim(to_str(curr_key)) // & - "/moments") + call read_dataset(moments_group, "order" // trim(to_str(k)), & + moment_name) k = k + 1 end do end do case default - call sp % read_data(moment_name, "order" // trim(to_str(k)), & - group="tallies/tally " // trim(to_str(curr_key)) // "/moments") + call read_dataset(moments_group, "order" // trim(to_str(k)), & + moment_name) k = k + 1 end select end do MOMENT_LOOP + call close_group(moments_group) + call close_group(tally_group) + end do TALLY_METADATA ! Check to make sure source bank is present @@ -908,22 +850,21 @@ contains if (master) then ! Read number of realizations for global tallies - call sp % read_data(n_realizations, "n_realizations", collect=.false.) + call read_dataset(file_id, "n_realizations", n_realizations, indep=.true.) ! Read number of global tallies - call sp % read_data(int_array(1), "n_global_tallies", collect=.false.) + call read_dataset(file_id, "n_global_tallies", int_array(1), indep=.false.) if (int_array(1) /= N_GLOBAL_TALLIES) then call fatal_error("Number of global tallies does not match in state & &point.") end if ! Read global tally data - call sp % read_tally_result(global_tallies, "global_tallies", & - n1=N_GLOBAL_TALLIES, n2=1) + call read_dataset(file_id, "global_tallies", global_tallies) ! Check if tally results are present - call sp % read_data(int_array(1), "tallies_present", & - group="tallies", collect=.false.) + tallies_group = open_group(file_id, "tallies") + call read_dataset(file_id, "tallies_present", int_array(1), indep=.true.) ! Read in sum and sum squared if (int_array(1) == 1) then @@ -934,13 +875,14 @@ contains curr_key = key_array(id_array(i)) ! Read sum and sum_sq for each bin - call sp % read_tally_result(tally % results, "results", & - group="tallies/tally " // trim(to_str(curr_key)), & - n1=size(tally % results, 1), n2=size(tally % results, 2)) - + tally_group = open_group(tallies_group, "tally " // & + trim(to_str(curr_key))) + call read_dataset(tally_group, "results", tally%results) + call close_group(tally_group) end do TALLY_RESULTS - end if + + call close_group(tallies_group) end if deallocate(id_array) @@ -953,30 +895,186 @@ contains if (.not. source_present) then ! Close statepoint file - call sp % file_close() + call file_close(file_id) ! Write message call write_message("Loading source file " // trim(path_source_point) & - &// "...", 1) + // "...", 1) ! Open source file - call sp % file_open(path_source_point, 'r', serial = .false.) + file_id = file_open(path_source_point, 'r', parallel=.true.) ! Read file type - call sp % read_data(int_array(1), "filetype") + call read_dataset(file_id, "filetype", int_array(1)) end if ! Write out source - call sp % read_source_bank() + call read_source_bank(file_id) end if ! Close file - call sp % file_close() + call file_close(file_id) end subroutine load_state_point +!=============================================================================== +! WRITE_SOURCE_BANK writes OpenMC source_bank data +!=============================================================================== + + subroutine write_source_bank(group_id) + use bank_header, only: Bank + + integer(HID_T), intent(in) :: group_id + + integer :: hdf5_err + integer(HID_T) :: dset ! data set handle + integer(HID_T) :: dspace ! data or file space handle + integer(HSIZE_T) :: dims(1) + type(c_ptr) :: f_ptr +#ifdef PHDF5 + integer :: data_xfer_mode + integer(HID_T) :: plist ! property list + integer(HID_T) :: memspace ! memory space handle + integer(HSIZE_T) :: offset(1) ! source data offset +#endif + +#ifdef PHDF5 + ! Set size of total dataspace for all procs and rank + dims(1) = n_particles + call h5screate_simple_f(1, dims, dspace, hdf5_err) + call h5dcreate_f(group_id, "source_bank", hdf5_bank_t, dspace, dset, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + + ! Create another data space but for each proc individually + dims(1) = work + call h5screate_simple_f(1, dims, memspace, hdf5_err) + + ! Get the individual local proc dataspace + call h5dget_space_f(dset, dspace, hdf5_err) + + ! Select hyperslab for this dataspace + offset(1) = work_index(rank) + call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims, hdf5_err) + + ! Set up the property list for parallel writing + call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) + call h5pset_dxpl_mpio_f(plist, H5FD_MPIO_COLLECTIVE_F, hdf5_err) + + ! Set up pointer to data + f_ptr = c_loc(source_bank) + + ! Write data to file in parallel + call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err, & + file_space_id=dspace, mem_space_id=memspace, & + xfer_prp=plist) + + ! Close all ids + call h5sclose_f(dspace, hdf5_err) + call h5sclose_f(memspace, hdf5_err) + call h5dclose_f(dset, hdf5_err) + call h5pclose_f(plist, hdf5_err) + +#else + + ! Set size + dims(1) = work + + ! Create dataspace + call h5screate_simple_f(1, dims, dspace, hdf5_err) + + ! Create dataset + call h5dcreate_f(group_id, "source_bank", hdf5_bank_t, & + dspace, dset, hdf5_err) + + ! Set up pointer to data + f_ptr = c_loc(source_bank) + + ! Write dataset to file + call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err) + + ! Close all ids + call h5dclose_f(dset, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + +#endif + + end subroutine write_source_bank + +!=============================================================================== +! READ_SOURCE_BANK reads OpenMC source_bank data +!=============================================================================== + + subroutine read_source_bank(group_id) + use bank_header, only: Bank + + integer(HID_T), intent(in) :: group_id + + integer :: hdf5_err + integer(HID_T) :: dset ! data set handle + type(c_ptr) :: f_ptr +#ifdef PHDF5 + integer :: data_xfer_mode + integer(HID_T) :: plist ! property list + integer(HID_T) :: dspace ! data space handle + integer(HID_T) :: memspace ! memory space handle + integer(HSIZE_T) :: offset(1) ! offset of data + integer(HSIZE_T) :: dims(1) +#endif + +#ifdef PHDF5 + + ! Open the dataset + call h5dopen_f(group_id, "source_bank", dset, hdf5_err) + + ! Create another data space but for each proc individually + dims(1) = work + call h5screate_simple_f(1, dims, memspace, hdf5_err) + + ! Get the individual local proc dataspace + call h5dget_space_f(dset, dspace, hdf5_err) + + ! Select hyperslab for this dataspace + offset(1) = work_index(rank) + call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims, hdf5_err) + + ! Set up the property list for parallel writing + call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) + call h5pset_dxpl_mpio_f(plist, H5FD_MPIO_COLLECTIVE_F, hdf5_err) + + ! Set up pointer to data + f_ptr = c_loc(source_bank) + + ! Read data from file in parallel + call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err, & + file_space_id=dspace, mem_space_id=memspace, & + xfer_prp=plist) + + ! Close all ids + call h5sclose_f(dspace, hdf5_err) + call h5sclose_f(memspace, hdf5_err) + call h5dclose_f(dset, hdf5_err) + call h5pclose_f(plist, hdf5_err) + +#else + + ! Open dataset + call h5dopen_f(group_id, "source_bank", dset, hdf5_err) + + ! Set up pointer to data + f_ptr = c_loc(source_bank) + + ! Read dataset from file + call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err) + + ! Close all ids + call h5dclose_f(dset, hdf5_err) + +#endif + + end subroutine read_source_bank + subroutine read_source ! TODO write this routine ! TODO what if n_particles does not match source bank diff --git a/src/track_output.F90 b/src/track_output.F90 index d3716e9b05..d4b2878009 100644 --- a/src/track_output.F90 +++ b/src/track_output.F90 @@ -6,19 +6,27 @@ module track_output use global - use output_interface, only: BinaryOutput + use hdf5_interface use particle_header, only: Particle use string, only: to_str - implicit none + use hdf5 - type, private :: TrackCoordinates + implicit none + private + + type TrackCoordinates real(8), allocatable :: coords(:,:) end type TrackCoordinates - type(TrackCoordinates), private, allocatable :: tracks(:) + type(TrackCoordinates), allocatable :: tracks(:) !$omp threadprivate(tracks) + public :: initialize_particle_track + public :: write_particle_track + public :: add_particle_track + public :: finalize_particle_track + contains !=============================================================================== @@ -43,19 +51,19 @@ contains ! Add another column to coords i = size(tracks) - if (allocated(tracks(i) % coords)) then - n_tracks = size(tracks(i) % coords, 2) + if (allocated(tracks(i)%coords)) then + n_tracks = size(tracks(i)%coords, 2) allocate(new_coords(3, n_tracks + 1)) - new_coords(:, 1:n_tracks) = tracks(i) % coords - call move_alloc(FROM=new_coords, TO=tracks(i) % coords) + new_coords(:, 1:n_tracks) = tracks(i)%coords + call move_alloc(FROM=new_coords, TO=tracks(i)%coords) else n_tracks = 0 - allocate(tracks(i) % coords(3, 1)) + allocate(tracks(i)%coords(3, 1)) end if ! Write current coordinates into the newest column. n_tracks = n_tracks + 1 - tracks(i) % coords(:, n_tracks) = p % coord(1) % xyz + tracks(i)%coords(:, n_tracks) = p%coord(1)%xyz end subroutine write_particle_track !=============================================================================== @@ -87,35 +95,32 @@ contains subroutine finalize_particle_track(p) type(Particle), intent(in) :: p - integer :: length(2) - character(MAX_FILE_LEN) :: fname - type(BinaryOutput) :: binout - integer :: i - integer, allocatable :: n_coords(:) integer :: n_particle_tracks + integer(HID_T) :: file_id + character(MAX_FILE_LEN) :: fname + integer, allocatable :: n_coords(:) fname = trim(path_output) // 'track_' // trim(to_str(current_batch)) & - // '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) & + // '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p%id)) & // '.h5' ! Determine total number of particles and number of coordinates for each n_particle_tracks = size(tracks) allocate(n_coords(n_particle_tracks)) do i = 1, n_particle_tracks - n_coords(i) = size(tracks(i) % coords, 2) + n_coords(i) = size(tracks(i)%coords, 2) end do !$omp critical (FinalizeParticleTrack) - call binout % file_create(fname) - call binout % write_data(n_particle_tracks, 'n_particles') - call binout % write_data(n_coords, 'n_coords', length=n_particle_tracks) + file_id = file_create(fname) + call write_dataset(file_id, 'n_particles', n_particle_tracks) + call write_dataset(file_id, 'n_coords', n_coords) do i = 1, n_particle_tracks - length(:) = [3, n_coords(i)] - call binout % write_data(tracks(i) % coords, 'coordinates_' // & - trim(to_str(i)), length=length) + call write_dataset(file_id, 'coordinates_' // trim(to_str(i)), & + tracks(i)%coords) end do - call binout % file_close() + call file_close(file_id) !$omp end critical (FinalizeParticleTrack) deallocate(tracks) end subroutine finalize_particle_track