!!****************************************************************************** !! !! This file is part of the AMUN source code, a program to perform !! Newtonian or relativistic magnetohydrodynamical simulations on uniform or !! adaptive mesh. !! !! Copyright (C) 2008-2021 Grzegorz Kowal !! !! This program is free software: you can redistribute it and/or modify !! it under the terms of the GNU General Public License as published by !! the Free Software Foundation, either version 3 of the License, or !! (at your option) any later version. !! !! This program is distributed in the hope that it will be useful, !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! GNU General Public License for more details. !! !! You should have received a copy of the GNU General Public License !! along with this program. If not, see . !! !!****************************************************************************** !! !! module: IO !! !! This module handles data storage and job restart from restart files. !! !! !!****************************************************************************** ! module io ! import external subroutines ! use blocks, only : pointer_meta #ifdef HDF5 use hdf5 , only : hid_t #endif /* HDF5 */ use timers, only : set_timer, start_timer, stop_timer ! module variables are not implicit by default ! implicit none ! subroutine interfaces ! interface read_snapshot_parameter module procedure read_snapshot_parameter_string module procedure read_snapshot_parameter_integer module procedure read_snapshot_parameter_double end interface interface write_attribute_xml module procedure write_attribute_xml_string module procedure write_attribute_xml_integer module procedure write_attribute_xml_double module procedure write_attribute_xml_file end interface #ifdef HDF5 interface read_snapshot_parameter_h5 module procedure read_snapshot_parameter_string_h5 module procedure read_snapshot_parameter_integer_h5 module procedure read_snapshot_parameter_double_h5 end interface interface restore_attribute_h5 module procedure restore_attribute_string_h5 module procedure restore_attribute_number_h5 end interface interface store_attribute_h5 module procedure store_attribute_string_h5 module procedure store_attribute_number_h5 end interface #endif /* HDF5 */ interface read_attribute #ifdef HDF5 module procedure read_scalar_attribute_integer_h5 module procedure read_scalar_attribute_double_h5 module procedure read_array_attribute_long_h5 module procedure read_array_attribute_complex_h5 #endif /* HDF5 */ end interface interface read_array #ifdef HDF5 module procedure read_1d_array_integer_h5 module procedure read_2d_array_integer_h5 #if NDIMS == 2 module procedure read_3d_array_integer_h5 #endif /* NDIMS == 2 */ module procedure read_4d_array_integer_h5 #if NDIMS == 3 module procedure read_5d_array_integer_h5 #endif /* NDIMS == 3 */ module procedure read_1d_array_double_h5 module procedure read_4d_array_double_h5 #endif /* HDF5 */ end interface ! timer indices ! integer , save :: iio ! MODULE PARAMETERS: ! ================= ! ! snapshot formats ! integer, parameter :: snapshot_xml = 0 #ifdef HDF5 integer, parameter :: snapshot_hdf5 = 1 #endif /* HDF5 */ ! snapshot_format - the format of snapshots; ! restart_format - the format of restart snapshots; ! integer, save :: snapshot_format = snapshot_xml integer, save :: restart_format = snapshot_xml ! respath - the directory from which the restart snapshots should be read; ! nrest - for job restarting, this is the number of restart snapshot; ! irest - the local counter for the restart snapshots; ! isnap - the local counter for the regular snapshots; ! ishift - the shift of the snapshot counter for restarting job with ! different snapshot interval; ! hrest - the execution time interval for restart snapshot storing ! (in hours); the minimum allowed value is 3 minutes; ! hsnap - the problem time interval for regular snapshot storing; ! tsnap - the next snapshot time; ! character(len=255), save :: respath = "./" integer , save :: nrest = -1 integer(kind=4) , save :: irest = 1 integer(kind=4) , save :: isnap = 0 integer(kind=4) , save :: ishift = 0 real(kind=8) , save :: hrest = 6.0d+00 real(kind=8) , save :: hsnap = 1.0d+00 real(kind=8) , save :: tsnap = 0.0d+00 ! flag indicating to store snapshots at exact intervals ! logical , save :: precise_snapshots = .false. ! a flag to determine if XDMF files should be generated ! logical , save :: with_xdmf = .false. ! the compression format and level of the XML+binary files ! character(len=255), save :: cformat = "none" ! compression format integer , save :: clevel = 3 ! compression level ! the suffix of binary files in the XML+binary format ! character(len=8) , save :: binary_file_suffix = ".bin" ! the type of digest to use ! integer , save :: hash_type = 1 character(len=8) , save :: hash_str = 'xxh64' #ifdef HDF5 ! compression type ! integer , parameter :: H5Z_DEFLATE = 1, H5Z_ZSTANDARD = 32015 ! compression type (0 for no compressions, 1 for deflate, 32015 for zstandard) ! integer , save :: compression = 0, hclevel = 3 ! HDF5 property object identifier ! integer(hid_t) , save :: prp_id ! local variables to store the number of processors ! integer(kind=4) , save :: nfiles = 1 ! array of pointer used during job restart ! type(pointer_meta), dimension(:), allocatable, save :: block_array #endif /* HDF5 */ ! by default everything is private ! private ! declare public subroutines ! public :: initialize_io, finalize_io, print_io public :: restart_snapshot_number, restart_from_snapshot public :: read_snapshot_parameter public :: read_restart_snapshot, write_restart_snapshot, write_snapshot public :: update_dtp !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine INITIALIZE_IO: ! ------------------------ ! ! Subroutine initializes module IO by setting its parameters. ! ! Arguments: ! ! verbose - flag determining if the subroutine should be verbose; ! status - return flag of the procedure execution status; ! !=============================================================================== ! subroutine initialize_io(verbose, status) use compression, only : set_compression, get_compression use hash , only : hash_xxh64 #ifdef XXHASH use hash , only : hash_xxh3 #endif /* XXHASH */ #ifdef HDF5 use hdf5 , only : hsize_t use hdf5 , only : H5P_DATASET_CREATE_F, H5Z_FLAG_OPTIONAL_F use hdf5 , only : h5open_f, h5zfilter_avail_f, h5pcreate_f use hdf5 , only : h5pset_deflate_f, h5pset_filter_f #endif /* HDF5 */ use helpers , only : print_message use mpitools , only : nproc use parameters , only : get_parameter implicit none logical, intent(in) :: verbose integer, intent(out) :: status logical :: test character(len=255) :: dname character(len=255) :: sformat = "xml" character(len=255) :: precise = "off" character(len=255) :: ghosts = "on" character(len=255) :: xdmf = "off" character(len=255) :: suffix = "" ! compression file suffix character(len=8) :: dtype = "xxh64" #ifdef HDF5 logical :: cmpstatus = .false. integer(hsize_t) :: cd_nelmts = 1 integer, dimension(1) :: cd_values = 3 #endif /* HDF5 */ #ifdef HDF5 character(len=*), parameter :: loc = 'IO::initialize_io()' #endif /* HDF5 */ !------------------------------------------------------------------------------- ! call set_timer('SNAPSHOTS I/O' , iio) status = 0 ! get module parameters ! call get_parameter("restart_path" , respath) call get_parameter("restart_number" , nrest ) call get_parameter("restart_interval" , hrest ) call get_parameter("snapshot_interval", hsnap ) call get_parameter("precise_snapshots", precise) call get_parameter("include_ghosts" , ghosts ) call get_parameter("generate_xdmf" , xdmf ) ! add slash at the end of respath if not present ! if (index(respath, '/', back = .true.) /= len(trim(respath))) then write(respath,"(a)") trim(adjustl(respath)) // '/' end if ! check the snapshot format ! call get_parameter("snapshot_format", sformat) select case(sformat) #ifdef HDF5 case('h5', 'hdf5', 'H5', 'HDF5') snapshot_format = snapshot_hdf5 #endif /* HDF5 */ case default snapshot_format = snapshot_xml end select ! check the restart snapshot format ! call get_parameter("restart_format", sformat) select case(sformat) #ifdef HDF5 case('h5', 'hdf5', 'H5', 'HDF5') restart_format = snapshot_hdf5 #endif /* HDF5 */ case default restart_format = snapshot_xml end select ! check the last available restart snapshot ! if (nrest == 0) then test = .true. nrest = 0 select case(restart_format) #ifdef HDF5 case(snapshot_hdf5) do while (test) nrest = nrest + 1 write(dname, "(a,'r',i6.6,'_',i5.5,'.h5')") & trim(respath), nrest, nproc inquire(file = dname, exist = test) end do #endif /* HDF5 */ case default do while (test) nrest = nrest + 1 write(dname, "(a,'restart-',i5.5)") trim(respath), nrest #ifdef __INTEL_COMPILER inquire(directory = dname, exist = test) #else /* __INTEL_COMPILER */ inquire(file = dname, exist = test) #endif /* __INTEL_COMPILER */ end do end select nrest = nrest - 1 end if ! get the compression format and level for XML+binary files ! call get_parameter("compression_format", cformat) call get_parameter("compression_level" , clevel) call set_compression(cformat, clevel, suffix) if (get_compression() > 0) & binary_file_suffix = ".bin" // trim(adjustl(suffix)) ! get the hash type ! call get_parameter("digest_type", dtype) select case(dtype) #ifdef XXHASH case('xxh3', 'XXH3') hash_type = hash_xxh3 hash_str = 'xxh3' #endif /* XXHASH */ case default hash_type = hash_xxh64 hash_str = 'xxh64' end select if (status == 0) then ! check if the snapshots should be stored at precise moments ! select case(trim(precise)) case ("off", "OFF", "n", "N", "false", "FALSE", "no", "NO") precise_snapshots = .false. case default precise_snapshots = .true. end select #ifdef HDF5 ! check if the XDMF files should be generated too ! select case(trim(xdmf)) case ("off", "OFF", "n", "N", "false", "FALSE", "no", "NO") with_xdmf = .false. case default with_xdmf = .true. end select ! initialize the HDF5 Fortran interface ! call h5open_f(status) if (status < 0) then call print_message(loc, & "Cannot initialize the HDF5 Fortran interface!") else ! prepare the property object for compression ! call h5pcreate_f(H5P_DATASET_CREATE_F, prp_id, status) if (status < 0) then call print_message(loc, & "Cannot create the compression property for datasets!") else ! detect available compression formats ! cmpstatus = .false. if (.not. cmpstatus) then call h5zfilter_avail_f(H5Z_ZSTANDARD, cmpstatus, status) if (cmpstatus) compression = H5Z_ZSTANDARD end if if (.not. cmpstatus) then call h5zfilter_avail_f(H5Z_DEFLATE, cmpstatus, status) if (cmpstatus) compression = H5Z_DEFLATE end if ! get the compression_level ! call get_parameter("compression_level", hclevel) ! initialize the compressor ! if (status == 0) then select case(compression) case(H5Z_ZSTANDARD) hclevel = max(1, min(20, hclevel)) cd_values(:) = hclevel call h5pset_filter_f(prp_id, H5Z_ZSTANDARD, H5Z_FLAG_OPTIONAL_F, & cd_nelmts, cd_values, status) case(H5Z_DEFLATE) hclevel = max(1, min(9, hclevel)) call h5pset_deflate_f(prp_id, hclevel, status) case default end select end if end if end if #endif /* HDF5 */ end if !------------------------------------------------------------------------------- ! end subroutine initialize_io ! !=============================================================================== ! ! subroutine FINALIZE_IO: ! ---------------------- ! ! Subroutine releases memory used by the module. ! ! Arguments: ! ! status - the subroutine call status; ! !=============================================================================== ! subroutine finalize_io(status) #ifdef HDF5 use hdf5 , only : h5pclose_f, h5close_f use helpers, only : print_message #endif /* HDF5 */ implicit none integer, intent(out) :: status #ifdef HDF5 character(len=*), parameter :: loc = 'IO::finalize_io()' #endif /* HDF5 */ !------------------------------------------------------------------------------- ! status = 0 #ifdef HDF5 call h5pclose_f(prp_id, status) if (status < 0) & call print_message(loc, "Could not close the HDF5 compression property!") call h5close_f(status) if (status < 0) & call print_message(loc, "Could not close the HDF5 Fortran interface!") #endif /* HDF5 */ !------------------------------------------------------------------------------- ! end subroutine finalize_io ! !=============================================================================== ! ! subroutine PRINT_IO: ! ------------------- ! ! Subroutine prints IO parameters. ! ! Arguments: ! ! verbose - flag determining if the subroutine should be verbose; ! !=============================================================================== ! subroutine print_io(verbose) use helpers, only : print_section, print_parameter implicit none logical, intent(in) :: verbose character(len=80) :: sfmt, msg integer :: dd, hh, mm, ss !------------------------------------------------------------------------------- ! if (verbose) then call print_section(verbose, "Snapshots") select case(snapshot_format) #ifdef HDF5 case(snapshot_hdf5) call print_parameter(verbose, "snapshot format", "HDF5") #endif /* HDF5 */ case default call print_parameter(verbose, "snapshot format", "XML+binary") end select select case(restart_format) #ifdef HDF5 case(snapshot_hdf5) call print_parameter(verbose, "restart snapshot format", "HDF5") #endif /* HDF5 */ case default call print_parameter(verbose, "restart snapshot format", "XML+binary") call print_parameter(verbose, "digest type", hash_str) call print_parameter(verbose, "compression format", cformat) call print_parameter(verbose, "compression level", clevel) end select call print_parameter(verbose, "precise snapshot intervals", & precise_snapshots, "on") #ifdef HDF5 select case(compression) case(H5Z_ZSTANDARD) call print_parameter(verbose, "HDF5 compression" , "zstd" ) call print_parameter(verbose, "compression level", hclevel ) case(H5Z_DEFLATE) call print_parameter(verbose, "HDF5 compression" , "deflate") call print_parameter(verbose, "compression level", hclevel ) case default call print_parameter(verbose, "HDF5 compression" , "none" ) end select call print_parameter(verbose, "generate XDMF files", with_xdmf, "on") #endif /* HDF5 */ call print_parameter(verbose, "snapshot interval" , hsnap) if (hrest > 0.0d+00) then dd = int(hrest / 2.4d+01) hh = int(mod(hrest, 2.4d+01)) mm = int(mod(6.0d+01 * hrest, 6.0d+01)) ss = int(mod(3.6d+03 * hrest, 6.0d+01)) sfmt = "(i2.2,'d',i2.2,'h',i2.2,'m',i2.2,'s')" write(msg,sfmt) dd, hh, mm, ss call print_parameter(verbose, "restart interval" , msg ) end if if (restart_from_snapshot()) then call print_parameter(verbose, "restart from path" , respath) call print_parameter(verbose, "restart from snapshot", nrest ) end if end if !------------------------------------------------------------------------------- ! end subroutine print_io ! !=============================================================================== ! ! function RESTART_SNAPSHOT_NUMBER: ! -------------------------------- ! ! Subroutine returns the number of restart snapshot. ! !=============================================================================== ! integer function restart_snapshot_number() implicit none !------------------------------------------------------------------------------- ! restart_snapshot_number = nrest !------------------------------------------------------------------------------- ! end function restart_snapshot_number ! !=============================================================================== ! ! function RESTART_FROM_SNAPSHOT: ! ------------------------------ ! ! Subroutine returns true if the current job is the restarted one. ! !=============================================================================== ! logical function restart_from_snapshot() implicit none !------------------------------------------------------------------------------- ! restart_from_snapshot = (nrest > 0) !------------------------------------------------------------------------------- ! end function restart_from_snapshot ! !=============================================================================== ! ! subroutine READ_RESTART_SNAPSHOT: ! -------------------------------- ! ! Subroutine reads restart snapshot files in order to resume the job. ! This is a wrapper calling specific format subroutine. ! ! Arguments: ! ! status - the status flag to inform if subroutine succeeded or failed; ! !=============================================================================== ! subroutine read_restart_snapshot(status) use evolution, only : time implicit none integer, intent(out) :: status !------------------------------------------------------------------------------- ! call start_timer(iio) status = 0 select case(restart_format) #ifdef HDF5 case(snapshot_hdf5) call read_restart_snapshot_h5(status) #endif /* HDF5 */ case default call read_restart_snapshot_xml(status) end select ! calculate the shift of the snapshot counter, and the next snapshot time ! ishift = int(time / hsnap) - isnap + 1 tsnap = (ishift + isnap) * hsnap call stop_timer(iio) !------------------------------------------------------------------------------- ! end subroutine read_restart_snapshot ! !=============================================================================== ! ! subroutine WRITE_RESTART_SNAPSHOT: ! --------------------------------- ! ! Subroutine stores current restart snapshot files. This is a wrapper ! calling specific format subroutine. ! ! Arguments: ! ! thrs - the current execution time in hours; ! problem - the problem's name; ! nrun - the run number; ! status - the status flag; ! !=============================================================================== ! subroutine write_restart_snapshot(thrs, problem, nrun, status) implicit none real(kind=8) , intent(in) :: thrs character(len=*), intent(in) :: problem integer , intent(in) :: nrun integer , intent(out) :: status !------------------------------------------------------------------------------- ! status = 0 ! check if conditions for storing the restart snapshot have been met ! if (hrest < 5.0d-02 .or. thrs < irest * hrest) return call start_timer(iio) select case(snapshot_format) #ifdef HDF5 case(snapshot_hdf5) call store_restart_snapshot_h5(problem, nrun, status) #endif /* HDF5 */ case default call write_restart_snapshot_xml(problem, nrun, status) end select ! increase the restart snapshot counter ! irest = irest + 1 call stop_timer(iio) !------------------------------------------------------------------------------- ! end subroutine write_restart_snapshot ! !=============================================================================== ! ! subroutine WRITE_SNAPSHOT: ! ------------------------- ! ! Subroutine stores block data in snapshots. Block variables are grouped ! together and stored in big 4D arrays separately. This is a wrapper for ! specific format storing. ! ! Arguments: ! ! problem - the problem's name; ! !=============================================================================== ! subroutine write_snapshot(problem) use evolution, only : time #ifdef HDF5 use mpitools , only : master #endif /* HDF5 */ implicit none character(len=*), intent(in) :: problem integer :: status !------------------------------------------------------------------------------- ! if (hsnap <= 0.0e+00 .or. time < tsnap) return call start_timer(iio) select case(snapshot_format) #ifdef HDF5 case(snapshot_hdf5) call store_snapshot_h5(problem, status) if (with_xdmf) then call write_snapshot_xdmf() if (master) call write_snapshot_xdmf_master() end if #endif /* HDF5 */ case default call write_snapshot_xml(problem, status) end select ! increase the snapshot counter and calculate the next snapshot time ! isnap = isnap + 1 tsnap = (ishift + isnap) * hsnap call stop_timer(iio) !------------------------------------------------------------------------------- ! end subroutine write_snapshot ! !=============================================================================== ! ! subroutine UPDATE_DTP: ! --------------------- ! ! Subroutine updates dtp from module EVOLUTION. ! !=============================================================================== ! subroutine update_dtp() use evolution, only : time, dtp implicit none !------------------------------------------------------------------------------- ! if (precise_snapshots) then if (tsnap > time) then dtp = tsnap - time else dtp = hsnap endif end if !------------------------------------------------------------------------------- ! end subroutine update_dtp ! !=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine READ_SNAPSHOT_PARAMETER_STRING: ! ----------------------------------------- ! ! Subroutine reads a string parameter from the restart snapshot. ! ! Arguments: ! ! pname - the parameter name; ! pvalue - the parameter value; ! status - the status flag (the success is 0, failure otherwise); ! !=============================================================================== ! subroutine read_snapshot_parameter_string(pname, pvalue, status) use helpers, only : print_message implicit none character(len=*), intent(in) :: pname character(len=*), intent(inout) :: pvalue integer , intent(inout) :: status logical :: test character(len=255) :: dname, fname, line integer(kind=4) :: lun = 103 integer :: l, u character(len=*), parameter :: loc = 'IO::read_snapshot_parameter_string' !------------------------------------------------------------------------------- ! status = 0 select case(restart_format) #ifdef HDF5 case(snapshot_hdf5) call read_snapshot_parameter_h5(pname, pvalue, status) #endif /* HDF5 */ case default ! check if the snapshot directory and metafile exist ! write(dname, "(a,'restart-',i5.5)") trim(respath), nrest #ifdef __INTEL_COMPILER inquire(directory = dname, exist = test) #else /* __INTEL_COMPILER */ inquire(file = dname, exist = test) #endif /* __INTEL_COMPILER */ if (.not. test) then call print_message(loc, trim(dname) // " does not exist!") status = 121 return end if write(fname,"(a,'/metadata.xml')") trim(dname) inquire(file = fname, exist = test) if (.not. test) then call print_message(loc, trim(fname) // " does not exist!") status = 121 return end if ! read requested parameter from the file ! open(newunit = lun, file = fname, status = 'old') 10 read(lun, fmt = "(a)", end = 20) line l = index(line, trim(adjustl(pname))) if (l > 0) then l = index(line, '>') + 1 u = index(line, '<', back = .true.) - 1 pvalue = trim(adjustl(line(l:u))) end if go to 10 20 close(lun) end select !------------------------------------------------------------------------------- ! end subroutine read_snapshot_parameter_string ! !=============================================================================== ! ! subroutine READ_SNAPSHOT_PARAMETER_INTEGER: ! ------------------------------------------ ! ! Subroutine reads an integer parameter from the restart snapshot. ! ! Arguments: ! ! pname - the parameter name; ! pvalue - the parameter value; ! status - the status flag (the success is 0, failure otherwise); ! !=============================================================================== ! subroutine read_snapshot_parameter_integer(pname, pvalue, status) use helpers, only : print_message implicit none character(len=*), intent(in) :: pname integer , intent(inout) :: pvalue integer , intent(inout) :: status logical :: test character(len=255) :: dname, fname, line, svalue integer(kind=4) :: lun = 103 integer :: l, u character(len=*), parameter :: loc = 'IO::read_snapshot_parameter_integer' !------------------------------------------------------------------------------- ! status = 0 select case(restart_format) #ifdef HDF5 case(snapshot_hdf5) call read_snapshot_parameter_h5(pname, pvalue, status) #endif /* HDF5 */ case default ! check if the snapshot directory and metafile exist ! write(dname, "(a,'restart-',i5.5)") trim(respath), nrest #ifdef __INTEL_COMPILER inquire(directory = dname, exist = test) #else /* __INTEL_COMPILER */ inquire(file = dname, exist = test) #endif /* __INTEL_COMPILER */ if (.not. test) then call print_message(loc, trim(dname) // " does not exist!") status = 121 return end if write(fname,"(a,'/metadata.xml')") trim(dname) inquire(file = fname, exist = test) if (.not. test) then call print_message(loc, trim(fname) // " does not exist!") status = 121 return end if ! read parameter from the file ! open(newunit = lun, file = fname, status = 'old') 10 read(lun, fmt = "(a)", end = 20) line l = index(line, trim(adjustl(pname))) if (l > 0) then l = index(line, '>') + 1 u = index(line, '<', back = .true.) - 1 svalue = trim(adjustl(line(l:u))) read(svalue, fmt = *) pvalue end if go to 10 20 close(lun) end select !------------------------------------------------------------------------------- ! end subroutine read_snapshot_parameter_integer ! !=============================================================================== ! ! subroutine READ_SNAPSHOT_PARAMETER_DOUBLE: ! ----------------------------------------- ! ! Subroutine reads a floating point parameter from the restart snapshot. ! ! Arguments: ! ! pname - the parameter name; ! pvalue - the parameter value; ! status - the status flag (the success is 0, failure otherwise); ! !=============================================================================== ! subroutine read_snapshot_parameter_double(pname, pvalue, status) use helpers, only : print_message implicit none character(len=*), intent(in) :: pname real(kind=8) , intent(inout) :: pvalue integer , intent(inout) :: status logical :: test character(len=255) :: dname, fname, line, svalue integer(kind=4) :: lun = 103 integer :: l, u character(len=*), parameter :: loc = 'IO::read_snapshot_parameter_double' !------------------------------------------------------------------------------- ! status = 0 select case(restart_format) #ifdef HDF5 case(snapshot_hdf5) call read_snapshot_parameter_h5(pname, pvalue, status) #endif /* HDF5 */ case default ! check if the snapshot directory and metafile exist ! write(dname, "(a,'restart-',i5.5)") trim(respath), nrest #ifdef __INTEL_COMPILER inquire(directory = dname, exist = test) #else /* __INTEL_COMPILER */ inquire(file = dname, exist = test) #endif /* __INTEL_COMPILER */ if (.not. test) then call print_message(loc, trim(dname) // " does not exist!") status = 121 return end if write(fname,"(a,'/metadata.xml')") trim(dname) inquire(file = fname, exist = test) if (.not. test) then call print_message(loc, trim(fname) // " does not exist!") status = 121 return end if ! read parameter from the file ! open(newunit = lun, file = fname, status = 'old') 10 read(lun, fmt = "(a)", end = 20) line l = index(line, trim(adjustl(pname))) if (l > 0) then l = index(line, '>') + 1 u = index(line, '<', back = .true.) - 1 svalue = trim(adjustl(line(l:u))) read(svalue, fmt = *) pvalue end if go to 10 20 close(lun) end select !------------------------------------------------------------------------------- ! end subroutine read_snapshot_parameter_double ! !=============================================================================== ! ! subroutine READ_RESTART_SNAPSHOT_XML: ! ------------------------------------ ! ! Subroutine reads restart snapshot, i.e. parameters, meta and data blocks ! stored in the XML+binary format restart files and reconstructs ! the data structure in order to resume a terminated job. ! ! Arguments: ! ! status - the return flag to inform if subroutine succeeded or failed; ! !=============================================================================== ! subroutine read_restart_snapshot_xml(status) use blocks , only : block_meta, block_data, pointer_meta, list_meta use blocks , only : ns => nsides, nc => nchildren, nregs use blocks , only : append_metablock, append_datablock, link_blocks use blocks , only : get_mblocks use blocks , only : set_last_id, get_last_id use blocks , only : metablock_set_id, metablock_set_process use blocks , only : metablock_set_refinement use blocks , only : metablock_set_configuration use blocks , only : metablock_set_level, metablock_set_position use blocks , only : metablock_set_coordinates, metablock_set_bounds use blocks , only : metablock_set_leaf use blocks , only : change_blocks_process use coordinates, only : nn => bcells, ncells, nghosts use coordinates, only : xmin, xmax, ymin, ymax, zmin, zmax use equations , only : cmax, cmax2 use evolution , only : step, time, dt, dth, dte use evolution , only : niterations, nrejections, errs use forcing , only : nmodes, fcoefs, einj use hash , only : check_hash, hash_xxh64 use helpers , only : print_message #ifdef MPI use mesh , only : redistribute_blocks #endif /* MPI */ use mpitools , only : nprocs, nproc use random , only : gentype, set_seeds implicit none integer, intent(out) :: status ! local variables ! logical :: test character(len=255) :: dname, fname, line, sname, svalue integer :: il, iu, nl, nx, nm, nd, nv, i, j, k, l, n, p, nu integer(kind=4) :: lndims, lnprocs, lnproc, lmblocks, lnleafs, llast_id integer(kind=4) :: ldblocks, lncells, lnghosts, lnseeds, lnmodes real(kind=8) :: deinj ! local pointers ! type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata ! local variables ! integer(kind=4) :: lun = 104, version = 1 integer(kind=8) :: digest, bytes ! hash strings ! character(len=16) :: hfield, hchild, hface, hedge, hcorner, hbound character(len=16) :: hids, harray, hseed, hforce ! local arrays ! integer(kind=8), dimension(:) , allocatable :: hprim, hcons integer(kind=4), dimension(:) , allocatable :: ids integer(kind=4), dimension(:,:) , allocatable :: fields integer(kind=4), dimension(:,:) , allocatable :: children #if NDIMS == 2 integer(kind=4), dimension(:,:,:,:) , allocatable :: edges integer(kind=4), dimension(:,:,:) , allocatable :: corners #endif /* NDIMS == 2 */ #if NDIMS == 3 integer(kind=4), dimension(:,:,:,:,:) , allocatable :: faces integer(kind=4), dimension(:,:,:,:,:) , allocatable :: edges integer(kind=4), dimension(:,:,:,:) , allocatable :: corners #endif /* NDIMS == 3 */ integer(kind=8), dimension(:,:) , allocatable :: seeds real(kind=8) , dimension(:,:,:) , allocatable :: bounds real(kind=8) , dimension(:,:,:,:,:,:), allocatable :: arrays real(kind=8) , dimension(:,:,:,:) , allocatable :: array ! array of pointer used during job restart ! type(pointer_meta), dimension(:), allocatable :: barray ! local parameters ! character(len=*), parameter :: loc = 'IO::read_restart_snapshot_xml()' character(len=*), parameter :: fmt = "(a,a,'_',i6.6,'.',a)" ! !------------------------------------------------------------------------------- ! status = 0 ! check if the snapshot directory exists ! write(dname, "(a,'restart-',i5.5)") trim(respath), nrest #ifdef __INTEL_COMPILER inquire(directory = dname, exist = test) #else /* __INTEL_COMPILER */ inquire(file = dname, exist = test) #endif /* __INTEL_COMPILER */ if (.not. test) then call print_message(loc, trim(dname) // " does not exist!") status = 121 return end if dname = trim(dname) // "/" write(fname,"(a,'metadata.xml')") trim(dname) inquire(file = fname, exist = test) if (.not. test) then call print_message(loc, trim(fname) // " does not exist!") status = 121 return end if ! read attributes from the metadata file ! open(newunit = lun, file = fname, status = 'old') 10 read(lun, fmt = "(a)", end = 20) line if (index(line, ' 0) then il = index(line, 'name="') if (il > 0) then il = il + 6 iu = index(line(il:), '"') + il - 2 write(sname,*) line(il:iu) il = index(line, '>') + 1 iu = index(line, '<', back = .true.) - 1 write(svalue,*) line(il:iu) select case(trim(adjustl(sname))) case('ndims') read(svalue, fmt = *) lndims case('nprocs') read(svalue, fmt = *) lnprocs case('nproc') read(svalue, fmt = *) lnproc case('mblocks') read(svalue, fmt = *) lmblocks case('dblocks') read(svalue, fmt = *) ldblocks case('nleafs') read(svalue, fmt = *) lnleafs case('last_id') read(svalue, fmt = *) llast_id case('ncells') read(svalue, fmt = *) lncells case('nghosts') read(svalue, fmt = *) lnghosts case('nseeds') read(svalue, fmt = *) lnseeds case('step') read(svalue, fmt = *) step case('isnap') read(svalue, fmt = *) isnap case('nvars') read(svalue, fmt = *) nv case('nmodes') read(svalue, fmt = *) lnmodes case('xmin') read(svalue, fmt = *) xmin case('xmax') read(svalue, fmt = *) xmax case('ymin') read(svalue, fmt = *) ymin case('ymax') read(svalue, fmt = *) ymax case('zmin') read(svalue, fmt = *) zmin case('zmax') read(svalue, fmt = *) zmax case('time') read(svalue, fmt = *) time case('dt') read(svalue, fmt = *) dt case('dth') read(svalue, fmt = *) dth case('dte') read(svalue, fmt = *) dte case('cmax') read(svalue, fmt = *) cmax cmax2 = cmax * cmax case('niterations') read(svalue, fmt = *) niterations case('nrejections') read(svalue, fmt = *) nrejections case('errs(1)') read(svalue, fmt = *) errs(1) case('errs(2)') read(svalue, fmt = *) errs(2) case('errs(3)') read(svalue, fmt = *) errs(3) case('fields') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) hfield case('children') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) hchild case('faces') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) hface case('edges') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) hedge case('corners') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) hcorner case('bounds') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) hbound case('forcing') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) hforce end select end if end if go to 10 20 close(lun) ! check the number of dimensions ! if (lndims /= NDIMS) then call print_message(loc, "The number of dimensions does not match!") return end if ! check the block dimensions ! if (lncells /= ncells) then call print_message(loc, "The block dimensions do not match!") return end if ! allocate all metablocks ! do l = 1, lmblocks call append_metablock(pmeta, status) end do ! check if the number of created metablocks is equal to lbmcloks ! if (lmblocks /= get_mblocks()) then call print_message(loc, "Number of metablocks does not match!") end if ! set the last_id ! call set_last_id(llast_id) ! get numbers of meta and data blocks ! nx = llast_id nm = lmblocks ! prepare and store metablocks ! allocate(barray(nx), fields(nm,14), children(nm,nc), bounds(nm,3,2), & #if NDIMS == 3 faces(nm,NDIMS,ns,ns,ns), & edges(nm,NDIMS,ns,ns,ns), corners(nm,ns,ns,ns), & #else /* NDIMS == 3 */ edges(nm,NDIMS,ns,ns), corners(nm,ns,ns), & #endif /* NDIMS == 3 */ stat = status) if (status == 0) then fields(:,:) = -1 children(:,:) = -1 #if NDIMS == 3 faces(:,:,:,:,:) = -1 edges(:,:,:,:,:) = -1 corners(:,:,:,:) = -1 #else /* NDIMS == 3 */ edges(:,:,:,:) = -1 corners(:,:,:) = -1 #endif /* NDIMS == 3 */ bounds(:,:,:) = 0.0d+00 ! read metablocks from binary files and check hashes ! bytes = size(transfer(fields, [ 0_1 ]), kind=8) write(fname,"(a,'metablock_fields.bin')") trim(dname) open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) fields close(lun) read(hfield, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(fields, 1_1, bytes), & digest, hash_xxh64) write(fname,"(a,'metablock_children.bin')") trim(dname) bytes = size(transfer(children, [ 0_1 ]), kind=8) open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) children close(lun) read(hchild, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(children, 1_1, bytes), & digest, hash_xxh64) #if NDIMS == 3 bytes = size(transfer(faces, [ 0_1 ]), kind=8) write(fname,"(a,'metablock_faces.bin')") trim(dname) open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) faces close(lun) read(hface, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(faces, 1_1, bytes), & digest, hash_xxh64) #endif /* NDIMS == 3 */ bytes = size(transfer(edges, [ 0_1 ]), kind=8) write(fname,"(a,'metablock_edges.bin')") trim(dname) open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) edges close(lun) read(hedge, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(edges, 1_1, bytes), & digest, hash_xxh64) bytes = size(transfer(corners, [ 0_1 ]), kind=8) write(fname,"(a,'metablock_corners.bin')") trim(dname) open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) corners close(lun) read(hcorner, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(corners, 1_1, bytes), & digest, hash_xxh64) bytes = size(transfer(bounds, [ 0_1 ]), kind=8) write(fname,"(a,'metablock_bounds.bin')") trim(dname) open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) bounds close(lun) read(hbound, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(bounds, 1_1, bytes), & digest, hash_xxh64) ! iterate over all meta blocks and restore their fields ! l = 0 pmeta => list_meta do while(associated(pmeta)) l = l + 1 barray(fields(l,1))%ptr => pmeta call metablock_set_id (pmeta, fields(l, 1)) call metablock_set_process (pmeta, fields(l, 2)) call metablock_set_level (pmeta, fields(l, 3)) call metablock_set_configuration(pmeta, fields(l, 4)) call metablock_set_refinement (pmeta, fields(l, 5)) call metablock_set_position (pmeta, fields(l, 6:8)) call metablock_set_coordinates (pmeta, fields(l, 9:11)) call metablock_set_bounds (pmeta, bounds(l,1,1), bounds(l,1,2), & bounds(l,2,1), bounds(l,2,2), & bounds(l,3,1), bounds(l,3,2)) if (fields(l,12) == 1) call metablock_set_leaf(pmeta) pmeta => pmeta%next end do ! over all meta blocks ! iterate over all meta blocks and restore their pointers ! l = 0 pmeta => list_meta do while(associated(pmeta)) l = l + 1 if (fields(l,14) > 0) pmeta%parent => barray(fields(l,14))%ptr do p = 1, nc if (children(l,p) > 0) then pmeta%child(p)%ptr => barray(children(l,p))%ptr end if end do ! p = 1, nc #if NDIMS == 2 do j = 1, ns do i = 1, ns do n = 1, NDIMS if (edges(l,n,i,j) > 0) & pmeta%edges(i,j,n)%ptr => barray(edges(l,n,i,j))%ptr end do ! NDIMS if (corners(l,i,j) > 0) & pmeta%corners(i,j)%ptr => barray(corners(l,i,j))%ptr end do ! i = 1, ns end do ! j = 1, ns #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, ns do j = 1, ns do i = 1, ns do n = 1, NDIMS if (faces(l,n,i,j,k) > 0) & pmeta%faces(i,j,k,n)%ptr => barray(faces(l,n,i,j,k))%ptr if (edges(l,n,i,j,k) > 0) & pmeta%edges(i,j,k,n)%ptr => barray(edges(l,n,i,j,k))%ptr end do ! NDIMS if (corners(l,i,j,k) > 0) & pmeta%corners(i,j,k)%ptr => barray(corners(l,i,j,k))%ptr end do ! i = 1, ns end do ! j = 1, ns end do ! k = 1, ns #endif /* NDIMS == 3 */ pmeta => pmeta%next end do ! over all meta blocks if (allocated(fields)) deallocate(fields) if (allocated(children)) deallocate(children) if (allocated(bounds)) deallocate(bounds) #if NDIMS == 3 if (allocated(faces)) deallocate(faces) #endif /* NDIMS == 3 */ if (allocated(edges)) deallocate(edges) if (allocated(corners)) deallocate(corners) end if ! allocation ! check the number of forcing modes ! if (lnmodes == nmodes) then if (lnmodes > 0) then bytes = size(transfer(fcoefs, [ 0_1 ]), kind=8) write(fname,"(a,'forcing_coefficients.bin')") trim(dname) open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) fcoefs close(lun) read(hforce, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(fcoefs, 1_1, bytes), & digest, hash_xxh64) end if else call print_message(loc, "The number of forcing modes does not match!") end if ! if the number of processes is bigger after the restart than before ! if (nprocs >= lnprocs) then if (nproc < lnprocs) then write(fname,fmt) trim(dname), "datablocks", nproc, "xml" inquire(file = fname, exist = test) if (.not. test) then write(*,*) trim(fname) // " does not exist!" status = 121 return end if ! read attributes from the metadata file ! open(newunit = lun, file = fname, status = 'old') 30 read(lun, fmt = "(a)", end = 40) line if (index(line, ' 0) then if (index(line, 'version="2.0"') > 0) version = 2 end if if (index(line, ' 0) then il = index(line, 'name="') if (il > 0) then il = il + 6 iu = index(line(il:), '"') + il - 2 write(sname,*) line(il:iu) il = index(line, '>') + 1 iu = index(line, '<', back = .true.) - 1 write(svalue,*) line(il:iu) select case(trim(adjustl(sname))) case('dblocks') read(svalue, fmt = *) nd if (version > 1) then if (nd > 0) then allocate(hprim(nd), hcons(nd), stat = status) end if end if case('einj') read(svalue, fmt = *) einj case('ids') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) hids case('arrays') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) harray case('seeds') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) hseed end select if (index(sname, 'prim') > 0) then read(sname(7:), fmt = *) l il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = "(1z16)") hprim(l) end if if (index(sname, 'cons') > 0) then read(sname(7:), fmt = *) l il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = "(1z16)") hcons(l) end if end if end if go to 30 40 close(lun) nm = lncells + 2 * lnghosts if (lnghosts >= nghosts) then il = 1 + (lnghosts - nghosts) iu = nm - (lnghosts - nghosts) else il = 1 + (nghosts - lnghosts) iu = nn - (nghosts - lnghosts) end if if (nd > 0) then if (version > 1) then #if NDIMS == 3 allocate(ids(nd), array(nv,nm,nm,nm), stat = status) #else /* NDIMS == 3 */ allocate(ids(nd), array(nv,nm,nm, 1), stat = status) #endif /* NDIMS == 3 */ if (status == 0) then bytes = size(transfer(ids, [ 0_1 ]), kind=8) write(fname, fmt) trim(dname), "datablock_ids", nproc, "bin" open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) ids close(lun) read(hids, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(ids, 1_1, bytes), & digest, hash_xxh64) bytes = size(transfer(array, [ 0_1 ]), kind=8) do l = 1, nd call append_datablock(pdata, status) call link_blocks(barray(ids(l))%ptr, pdata) write(fname,"(a,'datablock_prim_',i6.6,'_',i6.6,'.bin')") & trim(dname), nproc, l open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) array close(lun) call check_hash(loc, fname, transfer(array, 1_1, bytes), & hprim(l), hash_xxh64) if (lnghosts >= nghosts) then #if NDIMS == 3 pdata%q(:,:,:,:) = array(:,il:iu,il:iu,il:iu) #else /* NDIMS == 3 */ pdata%q(:,:,:,:) = array(:,il:iu,il:iu,:) #endif /* NDIMS == 3 */ else #if NDIMS == 3 pdata%q(:,il:iu,il:iu,il:iu) = array(:,:,:,:) #else /* NDIMS == 3 */ pdata%q(:,il:iu,il:iu,:) = array(:,:,:,:) #endif /* NDIMS == 3 */ end if write(fname,"(a,'datablock_cons_',i6.6,'_',i6.6,'.bin')") & trim(dname), nproc, l open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) array close(lun) call check_hash(loc, fname, transfer(array, 1_1, bytes), & hcons(l), hash_xxh64) if (lnghosts >= nghosts) then #if NDIMS == 3 pdata%uu(:,:,:,:,1) = array(:,il:iu,il:iu,il:iu) #else /* NDIMS == 3 */ pdata%uu(:,:,:,:,1) = array(:,il:iu,il:iu,:) #endif /* NDIMS == 3 */ else #if NDIMS == 3 pdata%uu(:,il:iu,il:iu,il:iu,1) = array(:,:,:,:) #else /* NDIMS == 3 */ pdata%uu(:,il:iu,il:iu,:,1) = array(:,:,:,:) #endif /* NDIMS == 3 */ end if end do if (allocated(ids)) deallocate(ids) if (allocated(array)) deallocate(array) if (allocated(hprim)) deallocate(hprim) if (allocated(hcons)) deallocate(hcons) end if ! allocation else #if NDIMS == 3 allocate(ids(nd), arrays(nd,3,nv,nm,nm,nm), stat = status) #else /* NDIMS == 3 */ allocate(ids(nd), arrays(nd,3,nv,nm,nm, 1), stat = status) #endif /* NDIMS == 3 */ if (status == 0) then bytes = size(transfer(ids, [ 0_1 ]), kind=8) write(fname, fmt) trim(dname), "datablock_ids", nproc, "bin" open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) ids close(lun) read(hids, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(ids, 1_1, bytes), & digest, hash_xxh64) bytes = size(transfer(arrays, [ 0_1 ]), kind=8) write(fname, fmt) trim(dname), "datablock_arrays", nproc, "bin" open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) arrays close(lun) read(harray, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(arrays, 1_1, bytes), & digest, hash_xxh64) do l = 1, nd call append_datablock(pdata, status) call link_blocks(barray(ids(l))%ptr, pdata) if (lnghosts >= nghosts) then #if NDIMS == 3 pdata%q (:,:,:,:) = arrays(l,1,:,il:iu,il:iu,il:iu) pdata%uu(:,:,:,:,1) = arrays(l,2,:,il:iu,il:iu,il:iu) if (nregs > 1) pdata%uu(:,:,:,:,2) = arrays(l,3,:,il:iu,il:iu,il:iu) #else /* NDIMS == 3 */ pdata%q (:,:,:,:) = arrays(l,1,:,il:iu,il:iu,:) pdata%uu(:,:,:,:,1) = arrays(l,2,:,il:iu,il:iu,:) if (nregs > 1) pdata%uu(:,:,:,:,2) = arrays(l,3,:,il:iu,il:iu,:) #endif /* NDIMS == 3 */ else #if NDIMS == 3 pdata%q (:,il:iu,il:iu,il:iu) = arrays(l,1,:,:,:,:) pdata%uu(:,il:iu,il:iu,il:iu,1) = arrays(l,2,:,:,:,:) if (nregs > 1) pdata%uu(:,il:iu,il:iu,il:iu,2) = arrays(l,3,:,:,:,:) #else /* NDIMS == 3 */ pdata%q (:,il:iu,il:iu,:) = arrays(l,1,:,:,:,:) pdata%uu(:,il:iu,il:iu,:,1) = arrays(l,2,:,:,:,:) if (nregs > 1) pdata%uu(:,il:iu,il:iu,:,2) = arrays(l,3,:,:,:,:) #endif /* NDIMS == 3 */ end if end do if (allocated(ids)) deallocate(ids) if (allocated(arrays)) deallocate(arrays) end if ! allocation end if ! version == 1 end if ! nd > 0 ! restore PRNG seeds ! allocate(seeds(4,lnseeds), stat = status) if (status == 0) then bytes = size(transfer(seeds, [ 0_1 ]), kind=8) write(fname, fmt) trim(dname), "random_seeds", nproc, "bin" open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) seeds close(lun) read(hseed, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(seeds, 1_1, bytes), & digest, hash_xxh64) call set_seeds(lnseeds, seeds(:,:), .false.) if (allocated(seeds)) deallocate(seeds) end if ! allocation else ! nproc < lnprocs write(fname,fmt) trim(dname), "datablocks", 0, "xml" inquire(file = fname, exist = test) if (.not. test) then write(*,*) trim(fname) // " does not exist!" status = 121 return end if ! read attributes from the metadata file ! open(newunit = lun, file = fname, status = 'old') 50 read(lun, fmt = "(a)", end = 60) line if (index(line, ' 0) then il = index(line, 'name="') if (il > 0) then il = il + 6 iu = index(line(il:), '"') + il - 2 write(sname,*) line(il:iu) il = index(line, '>') + 1 iu = index(line, '<', back = .true.) - 1 write(svalue,*) line(il:iu) select case(trim(adjustl(sname))) case('seeds') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) hseed end select end if end if go to 50 60 close(lun) ! restore PRNG seeds for remaining processes ! if (trim(gentype) == "same") then allocate(seeds(4,lnseeds), stat = status) if (status == 0) then bytes = size(transfer(seeds, [ 0_1 ]), kind=8) write(fname, fmt) trim(dname), "random_seeds", 0, "bin" open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) seeds close(lun) read(hseed, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(seeds, 1_1, bytes), & digest, hash_xxh64) call set_seeds(lnseeds, seeds(:,:), .false.) if (allocated(seeds)) deallocate(seeds) end if ! allocation end if ! gentype == "same" end if ! nproc < nprocs else ! nprocs < lnprocs ! divide files between processes ! nl = 0 i = mod(lnprocs, nprocs) j = lnprocs / nprocs do p = 0, nprocs k = 0 do n = 0, p nl = k if (n < i) then nu = k + j else nu = k + j - 1 end if k = nu + 1 end do do n = nl, nu call change_blocks_process(n, p) end do end do k = 0 do n = 0, nproc nl = k if (n < i) then nu = k + j else nu = k + j - 1 end if k = nu + 1 end do do n = nl, nu write(fname,fmt) trim(dname), "datablocks", n, "xml" inquire(file = fname, exist = test) if (.not. test) then write(*,*) trim(fname) // " does not exist!" status = 121 return end if ! read attributes from the metadata file ! open(newunit = lun, file = fname, status = 'old') 70 read(lun, fmt = "(a)", end = 80) line if (index(line, ' 0) then il = index(line, 'name="') if (il > 0) then il = il + 6 iu = index(line(il:), '"') + il - 2 write(sname,*) line(il:iu) il = index(line, '>') + 1 iu = index(line, '<', back = .true.) - 1 write(svalue,*) line(il:iu) select case(trim(adjustl(sname))) case('dblocks') read(svalue, fmt = *) nd case('einj') read(svalue, fmt = *) deinj case('ids') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) hids case('arrays') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) harray case('seeds') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 read(line(il:iu), fmt = *) hseed end select end if end if go to 70 80 close(lun) einj = einj + deinj nm = lncells + 2 * lnghosts if (lnghosts >= nghosts) then il = 1 + (lnghosts - nghosts) iu = nm - (lnghosts - nghosts) else il = 1 + (nghosts - lnghosts) iu = nn - (nghosts - lnghosts) end if #if NDIMS == 3 allocate(ids(nd), arrays(nd,3,nv,nm,nm,nm), stat = status) #else /* NDIMS == 3 */ allocate(ids(nd), arrays(nd,3,nv,nm,nm, 1), stat = status) #endif /* NDIMS == 3 */ if (status == 0) then bytes = size(transfer(ids, [ 0_1 ]), kind=8) write(fname, fmt) trim(dname), "datablock_ids", n, "bin" open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) ids close(lun) read(hids, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(ids, 1_1, bytes), & digest, hash_xxh64) bytes = size(transfer(arrays, [ 0_1 ]), kind=8) write(fname, fmt) trim(dname), "datablock_arrays", n, "bin" open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) arrays close(lun) read(harray, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(arrays, 1_1, bytes), & digest, hash_xxh64) do l = 1, nd call append_datablock(pdata, status) call link_blocks(barray(ids(l))%ptr, pdata) if (lnghosts >= nghosts) then #if NDIMS == 3 pdata%q (:,:,:,:) = arrays(l,1,:,il:iu,il:iu,il:iu) pdata%uu(:,:,:,:,1) = arrays(l,2,:,il:iu,il:iu,il:iu) if (nregs > 1) pdata%uu(:,:,:,:,2) = arrays(l,3,:,il:iu,il:iu,il:iu) #else /* NDIMS == 3 */ pdata%q (:,:,:,:) = arrays(l,1,:,il:iu,il:iu,:) pdata%uu(:,:,:,:,1) = arrays(l,2,:,il:iu,il:iu,:) if (nregs > 1) pdata%uu(:,:,:,:,2) = arrays(l,3,:,il:iu,il:iu,:) #endif /* NDIMS == 3 */ else #if NDIMS == 3 pdata%q (:,il:iu,il:iu,il:iu) = arrays(l,1,:,:,:,:) pdata%uu(:,il:iu,il:iu,il:iu,1) = arrays(l,2,:,:,:,:) if (nregs > 1) pdata%uu(:,il:iu,il:iu,il:iu,2) = arrays(l,3,:,:,:,:) #else /* NDIMS == 3 */ pdata%q (:,il:iu,il:iu,:) = arrays(l,1,:,:,:,:) pdata%uu(:,il:iu,il:iu,:,1) = arrays(l,2,:,:,:,:) if (nregs > 1) pdata%uu(:,il:iu,il:iu,:,2) = arrays(l,3,:,:,:,:) #endif /* NDIMS == 3 */ end if end do if (allocated(ids)) deallocate(ids) if (allocated(arrays)) deallocate(arrays) end if ! allocation end do ! n = nl, nu ! restore seeds ! allocate(seeds(4,lnseeds), stat = status) if (status == 0) then bytes = size(transfer(seeds, [ 0_1 ]), kind=8) write(fname, fmt) trim(dname), "random_seeds", nproc, "bin" open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream') read(lun) seeds close(lun) read(hseed, fmt = "(1z16)") digest call check_hash(loc, fname, transfer(seeds, 1_1, bytes), & digest, hash_xxh64) call set_seeds(lnseeds, seeds(:,:), .false.) if (allocated(seeds)) deallocate(seeds) end if ! allocation end if ! nprocs >= lnprocs if (allocated(barray)) deallocate(barray) #ifdef MPI ! redistribute blocks between processors ! call redistribute_blocks() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine read_restart_snapshot_xml ! !=============================================================================== ! ! subroutine WRITE_RESTART_SNAPSHOT_XML: ! ------------------------------------- ! ! Subroutine saves a restart snapshot, i.e. parameters, meta and data blocks ! using the XML format for parameters and binary format for meta and data ! block fields. ! ! Arguments: ! ! problem - the problem's name; ! nrun - the snapshot number; ! status - the status flag to inform if subroutine succeeded or failed; ! !=============================================================================== ! subroutine write_restart_snapshot_xml(problem, nrun, status) use blocks , only : block_meta, block_data, list_meta, list_data use blocks , only : get_mblocks, get_dblocks, get_nleafs use blocks , only : get_last_id use blocks , only : ns => nsides, nc => nchildren use coordinates, only : nn => bcells, ncells, nghosts, minlev, maxlev use coordinates, only : xmin, xmax, ymin, ymax #if NDIMS == 3 use coordinates, only : zmin, zmax #endif /* NDIMS == 3 */ use coordinates, only : bdims => domain_base_dims use equations , only : eqsys, eos, nv, cmax use evolution , only : step, time, dt, dth, dte, cfl, glm_alpha, errs use evolution , only : atol, rtol, mrej, niterations, nrejections use forcing , only : nmodes, fcoefs, einj use hash , only : hash_xxh64, hash_string use helpers , only : print_message use mpitools , only : nprocs, nproc use parameters , only : get_parameter_file use random , only : gentype, nseeds, get_seeds implicit none ! input and output arguments ! character(len=*), intent(in) :: problem integer , intent(in) :: nrun integer , intent(out) :: status ! local variables ! logical :: test character(len=64) :: dname, fname, aname, hname integer(kind=8) :: digest, bytes integer(kind=4) :: lun = 103 integer :: nd, nl, nm, nx, i, j, l, n, p #if NDIMS == 3 integer :: k #endif /* NDIMS == 3 */ ! local pointers ! type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata ! local arrays ! integer(kind=4), dimension(:) , allocatable :: ids integer(kind=4), dimension(:,:) , allocatable :: fields integer(kind=4), dimension(:,:) , allocatable :: children #if NDIMS == 2 integer(kind=4), dimension(:,:,:,:) , allocatable :: edges integer(kind=4), dimension(:,:,:) , allocatable :: corners #endif /* NDIMS == 2 */ #if NDIMS == 3 integer(kind=4), dimension(:,:,:,:,:) , allocatable :: faces integer(kind=4), dimension(:,:,:,:,:) , allocatable :: edges integer(kind=4), dimension(:,:,:,:) , allocatable :: corners #endif /* NDIMS == 3 */ integer(kind=8), dimension(:,:) , allocatable :: seeds real(kind=8) , dimension(:,:,:) , allocatable :: bounds ! local parameters ! character(len=*), parameter :: loc = "IO::write_restart_snapshot_xml()" character(len=*), parameter :: fmt = "(a,a,'_',i6.6,'.',a)" ! !------------------------------------------------------------------------------- ! status = 0 hname = hash_string(hash_xxh64) ! create snapshot directory, if it does not exist ! write(dname, "('restart-',i5.5)") nrun #ifdef __INTEL_COMPILER inquire(directory = dname, exist = test) do while(.not. test) if (.not. test) call system("mkdir -p " // trim(dname)) inquire(directory = dname, exist = test) end do #else /* __INTEL_COMPILER */ inquire(file = dname, exist = test) do while(.not. test) if (.not. test) call system("mkdir -p " // trim(dname)) inquire(file = dname, exist = test) end do #endif /* __INTEL_COMPILER */ dname = trim(dname) // "/" ! get numbers of meta and data blocks ! nx = get_last_id() nm = get_mblocks() nd = get_dblocks() nl = get_nleafs() ! only process 0 stores the metadata ! if (nproc == 0) then ! copy the parameter file to the restart directory ! call get_parameter_file(fname, status) if (status == 0) then call system("cp -a " // trim(fname) // " " // trim(dname)) else call print_message(loc, "Cannot get the location of parameter file!") return end if ! store metadata (parameters and attributes) ! write(fname,"(a,'metadata.xml')") trim(dname) open(newunit = lun, file = fname, status = 'replace') write(lun,"(a)") "" write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "problem" , problem) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "nprocs" , nprocs) call write_attribute_xml(lun, "nproc" , nproc) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "eqsys" , eqsys) call write_attribute_xml(lun, "eos" , eos) call write_attribute_xml(lun, "nvars" , nv) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "ndims" , NDIMS) call write_attribute_xml(lun, "xblocks" , bdims(1)) call write_attribute_xml(lun, "yblocks" , bdims(2)) #if NDIMS == 3 call write_attribute_xml(lun, "zblocks" , bdims(3)) #endif /* NDIMS */ call write_attribute_xml(lun, "xmin" , xmin) call write_attribute_xml(lun, "xmax" , xmax) call write_attribute_xml(lun, "ymin" , ymin) call write_attribute_xml(lun, "ymax" , ymax) #if NDIMS == 3 call write_attribute_xml(lun, "zmin" , zmin) call write_attribute_xml(lun, "zmax" , zmax) #endif /* NDIMS */ write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "minlev" , minlev) call write_attribute_xml(lun, "maxlev" , maxlev) call write_attribute_xml(lun, "ncells" , ncells) call write_attribute_xml(lun, "nghosts" , nghosts) call write_attribute_xml(lun, "bcells" , nn) call write_attribute_xml(lun, "nchildren", nc) call write_attribute_xml(lun, "mblocks" , nm) call write_attribute_xml(lun, "nleafs" , nl) call write_attribute_xml(lun, "last_id" , nx) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "step" , step) call write_attribute_xml(lun, "time" , time) call write_attribute_xml(lun, "dt" , dt) call write_attribute_xml(lun, "dth" , dth) call write_attribute_xml(lun, "dte" , dte) call write_attribute_xml(lun, "cfl" , cfl) call write_attribute_xml(lun, "cmax" , cmax) call write_attribute_xml(lun, "glm_alpha", glm_alpha) call write_attribute_xml(lun, "absolute_tolerance", atol) call write_attribute_xml(lun, "relative_tolerance", rtol) call write_attribute_xml(lun, "maximum_rejections", mrej) call write_attribute_xml(lun, "niterations", niterations) call write_attribute_xml(lun, "nrejections", nrejections) call write_attribute_xml(lun, "errs(1)", errs(1)) call write_attribute_xml(lun, "errs(2)", errs(2)) call write_attribute_xml(lun, "errs(3)", errs(3)) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "nmodes" , nmodes) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "gentype" , gentype) call write_attribute_xml(lun, "nseeds" , nseeds) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "isnap" , isnap) write(lun,"(a)") '' write(lun,"(a)") '' ! prepare and store metablocks ! allocate(fields(nm,14), children(nm,nc), bounds(nm,3,2), & #if NDIMS == 3 faces(nm,NDIMS,ns,ns,ns), & edges(nm,NDIMS,ns,ns,ns), corners(nm,ns,ns,ns), & #else /* NDIMS == 3 */ edges(nm,NDIMS,ns,ns), corners(nm,ns,ns), & #endif /* NDIMS == 3 */ stat = status) if (status == 0) then fields(:,:) = -1 children(:,:) = -1 #if NDIMS == 3 faces(:,:,:,:,:) = -1 edges(:,:,:,:,:) = -1 corners(:,:,:,:) = -1 #else /* NDIMS == 3 */ edges(:,:,:,:) = -1 corners(:,:,:) = -1 #endif /* NDIMS == 3 */ bounds(:,:,:) = 0.0d+00 l = 0 pmeta => list_meta do while(associated(pmeta)) l = l + 1 fields(l, 1) = pmeta%id fields(l, 2) = pmeta%process fields(l, 3) = pmeta%level fields(l, 4) = pmeta%conf fields(l, 5) = pmeta%refine fields(l, 6) = pmeta%pos(1) fields(l, 7) = pmeta%pos(2) #if NDIMS == 3 fields(l, 8) = pmeta%pos(3) #endif /* NDIMS == 3 */ fields(l, 9) = pmeta%coords(1) fields(l,10) = pmeta%coords(2) #if NDIMS == 3 fields(l,11) = pmeta%coords(3) #endif /* NDIMS == 3 */ if (pmeta%leaf) fields(l,12) = 1 if (associated(pmeta%data) ) fields(l,13) = 1 if (associated(pmeta%parent)) fields(l,14) = pmeta%parent%id do p = 1, nc if (associated(pmeta%child(p)%ptr)) & children(l,p) = pmeta%child(p)%ptr%id end do #if NDIMS == 2 do j = 1, ns do i = 1, ns do n = 1, NDIMS if (associated(pmeta%edges(i,j,n)%ptr)) & edges(l,n,i,j) = pmeta%edges(i,j,n)%ptr%id end do ! NDIMS if (associated(pmeta%corners(i,j)%ptr)) & corners(l,i,j) = pmeta%corners(i,j)%ptr%id end do ! i = 1, ns end do ! j = 1, ns #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, ns do j = 1, ns do i = 1, ns do n = 1, NDIMS if (associated(pmeta%faces(i,j,k,n)%ptr)) & faces(l,n,i,j,k) = pmeta%faces(i,j,k,n)%ptr%id if (associated(pmeta%edges(i,j,k,n)%ptr)) & edges(l,n,i,j,k) = pmeta%edges(i,j,k,n)%ptr%id end do ! NDIMS if (associated(pmeta%corners(i,j,k)%ptr)) & corners(l,i,j,k) = pmeta%corners(i,j,k)%ptr%id end do ! i = 1, ns end do ! j = 1, ns end do ! k = 1, ns #endif /* NDIMS == 3 */ bounds(l,1,1) = pmeta%xmin bounds(l,1,2) = pmeta%xmax bounds(l,2,1) = pmeta%ymin bounds(l,2,2) = pmeta%ymax #if NDIMS == 3 bounds(l,3,1) = pmeta%zmin bounds(l,3,2) = pmeta%zmax #endif /* NDIMS == 3 */ pmeta => pmeta%next end do ! metablocks ! store metablock data ! write(fname,"(a,'.bin')") "metablock_fields" call write_binary_xml(dname, fname, transfer(fields, [ 0_1 ]), & hash_xxh64, bytes, digest) call write_attribute_xml(lun, "fields", fname, hname, bytes, digest) write(fname,"(a,'.bin')") "metablock_children" call write_binary_xml(dname, fname, transfer(children, [ 0_1 ]), & hash_xxh64, bytes, digest) call write_attribute_xml(lun, "children", fname, hname, bytes, digest) #if NDIMS == 3 write(fname,"(a,'.bin')") "metablock_faces" call write_binary_xml(dname, fname, transfer(faces, [ 0_1 ]), & hash_xxh64, bytes, digest) call write_attribute_xml(lun, "faces", fname, hname, bytes, digest) #endif /* NDIMS == 3 */ write(fname,"(a,'.bin')") "metablock_edges" call write_binary_xml(dname, fname, transfer(edges, [ 0_1 ]), & hash_xxh64, bytes, digest) call write_attribute_xml(lun, "edges", fname, hname, bytes, digest) write(fname,"(a,'.bin')") "metablock_corners" call write_binary_xml(dname, fname, transfer(corners, [ 0_1 ]), & hash_xxh64, bytes, digest) call write_attribute_xml(lun, "corners", trim(fname), trim(hname), & bytes, digest) write(fname,"(a,'.bin')") "metablock_bounds" call write_binary_xml(trim(dname), trim(fname), & transfer(bounds, [ 0_1 ]), hash_xxh64, bytes, digest) call write_attribute_xml(lun, "bounds", trim(fname), trim(hname), & bytes, digest) if (nmodes > 0) then write(fname,"(a,'.bin')") "forcing_coefficients" call write_binary_xml(trim(dname), trim(fname), & transfer(fcoefs, [ 0_1 ]), & hash_xxh64, bytes, digest) call write_attribute_xml(lun, "forcing", trim(fname), trim(hname), & bytes, digest) end if if (allocated(fields)) deallocate(fields) if (allocated(children)) deallocate(children) if (allocated(bounds)) deallocate(bounds) #if NDIMS == 3 if (allocated(faces)) deallocate(faces) #endif /* NDIMS == 3 */ if (allocated(edges)) deallocate(edges) if (allocated(corners)) deallocate(corners) else call print_message(loc, "Cannot allocate space for metablocks!") status = 1001 return end if ! allocation #if NDIMS == 3 #endif /* NDIMS == 3 */ write(lun,"(a)") '' write(lun,"(a)") '' close(lun) end if ! meta data file is stored only on the master process ! prepare and store data block info ! write(fname,fmt) trim(dname), "datablocks", nproc, "xml" open(newunit = lun, file = fname, status = 'replace') write(lun,"(a)") "" write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "nprocs" , nprocs) call write_attribute_xml(lun, "nproc" , nproc) call write_attribute_xml(lun, "ndims" , NDIMS) call write_attribute_xml(lun, "ncells" , ncells) call write_attribute_xml(lun, "nghosts", nghosts) call write_attribute_xml(lun, "bcells" , nn) call write_attribute_xml(lun, "dblocks", nd) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "einj" , einj) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "gentype", gentype) call write_attribute_xml(lun, "nseeds" , nseeds) write(lun,"(a)") '' write(lun,"(a)") '' if (nd > 0) then allocate(ids(nd), stat = status) if (status == 0) then l = 0 pdata => list_data do while(associated(pdata)) l = l + 1 ids(l) = pdata%meta%id write(aname,"('_',i6.6)") l write(fname,"('datablock_prim_',i6.6,a,'.bin')") nproc, trim(aname) call write_binary_xml(trim(dname), trim(fname), & transfer(pdata%q (:,:,:,:), [ 0_1 ]), & hash_xxh64, bytes, digest) call write_attribute_xml(lun, "prim" // trim(aname), trim(fname), & trim(hname), bytes, digest) write(fname,"('datablock_cons_',i6.6,a,'.bin')") nproc, trim(aname) call write_binary_xml(trim(dname), trim(fname), & transfer(pdata%uu(:,:,:,:,1), [ 0_1 ]), & hash_xxh64, bytes, digest) call write_attribute_xml(lun, "cons" // trim(aname), trim(fname), & trim(hname), bytes, digest) pdata => pdata%next end do ! data blocks ! store block IDs and arrays ! write(fname,"(a,'_',a,'_',i6.6,'.bin')") "datablock", "ids", nproc call write_binary_xml(trim(dname), trim(fname), & transfer(ids, [ 0_1 ]), hash_xxh64, bytes, digest) call write_attribute_xml(lun, "ids", trim(fname), trim(hname), & bytes, digest) if (allocated(ids)) deallocate(ids) else call print_message(loc, "Cannot allocate space for datablocks!") status = 1001 return end if ! allocation end if ! store PRNG seeds ! allocate(seeds(4,nseeds), stat = status) if (status == 0) then call get_seeds(seeds(:,:)) ! store seeds ! write(fname,"(a,'_',a,'_',i6.6,'.bin')") "random", "seeds", nproc call write_binary_xml(trim(dname), trim(fname), & transfer(seeds, [ 0_1 ]), hash_xxh64, bytes, digest) call write_attribute_xml(lun, "seeds", trim(fname), trim(hname), & bytes, digest) if (allocated(seeds)) deallocate(seeds) else call print_message(loc, "Cannot allocate space for random generator seeds!") status = 1001 return end if write(lun,"(a)") '' write(lun,"(a)") '' close(lun) !------------------------------------------------------------------------------- ! end subroutine write_restart_snapshot_xml ! !=============================================================================== ! ! subroutine WRITE_SNAPSHOT_XML: ! ----------------------------- ! ! Subroutine saves a regular snapshot, i.e. parameters, leafs and data blocks ! using the XML format for metadata and binary format for meta and data ! block fields. ! ! Arguments: ! ! problem - the problem's name; ! status - the status flag to inform if subroutine succeeded or failed; ! !=============================================================================== ! subroutine write_snapshot_xml(problem, status) use blocks , only : block_meta, block_data, list_meta, list_data use blocks , only : get_dblocks, get_nleafs use coordinates, only : nn => bcells, ncells, nghosts, minlev, maxlev use coordinates, only : xmin, xmax, ymin, ymax #if NDIMS == 3 use coordinates, only : zmin, zmax #endif /* NDIMS == 3 */ use coordinates, only : bdims => domain_base_dims use equations , only : eqsys, eos, nv, pvars, adiabatic_index, csnd use evolution , only : step, time, dt, cfl, glm_alpha use helpers , only : print_message use mpitools , only : nprocs, nproc use parameters , only : get_parameter_file use sources , only : viscosity, resistivity implicit none ! input and output arguments ! character(len=*), intent(in) :: problem integer , intent(out) :: status ! local variables ! logical :: test character(len=64) :: dname, fname character(len=256) :: vars integer(kind=8) :: dbytes = 0_8, ddigest = 0_8 integer(kind=8) :: cbytes = 0_8, cdigest = 0_8 integer(kind=4) :: lun = 103 integer :: nd, nl, l, p ! local pointers ! type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata ! local arrays ! integer(kind=4), dimension(:) , allocatable :: ids integer(kind=4), dimension(:,:) , allocatable :: fields real(kind=8) , dimension(:,:,:) , allocatable :: bounds real(kind=8) , dimension(:,:,:,:), allocatable :: array ! local parameters ! character(len=*), parameter :: loc = "IO::write_snapshot_xml()" character(len=*), parameter :: fmt = "(a,a,'_',i6.6,'.',a)" ! !------------------------------------------------------------------------------- ! status = 0 ! create snapshot directory, if it does not exist ! write(dname, "('snapshot-',i9.9)") isnap #ifdef __INTEL_COMPILER inquire(directory = dname, exist = test) do while(.not. test) if (.not. test) call system("mkdir -p " // trim(dname)) inquire(directory = dname, exist = test) end do #else /* __INTEL_COMPILER */ inquire(file = dname, exist = test) do while(.not. test) if (.not. test) call system("mkdir -p " // trim(dname)) inquire(file = dname, exist = test) end do #endif /* __INTEL_COMPILER */ dname = trim(dname) // "/" ! get numbers of meta and data blocks, leafs, and prepare stored variables ! nd = get_dblocks() nl = get_nleafs() vars = "" do l = 1, nv vars = trim(vars) // " " // trim(pvars(l)) end do ! only process 0 stores the metadata ! if (nproc == 0) then ! copy the parameter file to the restart directory ! call get_parameter_file(fname, status) if (status == 0) then call system("cp -a " // trim(fname) // " " // trim(dname)) else call print_message(loc, "Cannot get the location of parameter file!") return end if ! store metadata (parameters and attributes) ! write(fname,"(a,'metadata.xml')") trim(dname) open(newunit = lun, file = fname, status = 'replace') write(lun,"(a)") "" write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "problem" , problem) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "nprocs" , nprocs) call write_attribute_xml(lun, "nproc" , nproc) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "eqsys" , eqsys) call write_attribute_xml(lun, "eos" , eos) call write_attribute_xml(lun, "nvars" , nv) call write_attribute_xml(lun, "adiabatic_index", adiabatic_index) call write_attribute_xml(lun, "sound_speed" , csnd) call write_attribute_xml(lun, "viscosity" , viscosity) call write_attribute_xml(lun, "resistivity" , resistivity) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "ndims" , NDIMS) call write_attribute_xml(lun, "xblocks" , bdims(1)) call write_attribute_xml(lun, "yblocks" , bdims(2)) #if NDIMS == 3 call write_attribute_xml(lun, "zblocks" , bdims(3)) #endif /* NDIMS */ call write_attribute_xml(lun, "xmin" , xmin) call write_attribute_xml(lun, "xmax" , xmax) call write_attribute_xml(lun, "ymin" , ymin) call write_attribute_xml(lun, "ymax" , ymax) #if NDIMS == 3 call write_attribute_xml(lun, "zmin" , zmin) call write_attribute_xml(lun, "zmax" , zmax) #endif /* NDIMS */ write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "minlev" , minlev) call write_attribute_xml(lun, "maxlev" , maxlev) call write_attribute_xml(lun, "ncells" , ncells) call write_attribute_xml(lun, "nghosts" , nghosts) call write_attribute_xml(lun, "bcells" , nn) call write_attribute_xml(lun, "nleafs" , nl) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "step" , step) call write_attribute_xml(lun, "time" , time) call write_attribute_xml(lun, "dt" , dt) call write_attribute_xml(lun, "cfl" , cfl) call write_attribute_xml(lun, "glm_alpha", glm_alpha) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "isnap" , isnap) call write_attribute_xml(lun, "variables", trim(vars)) write(lun,"(a)") '' write(lun,"(a)") '' ! prepare and store metablocks ! allocate(fields(nl,5), bounds(nl,3,2), stat = status) if (status == 0) then fields(:,:) = -1 bounds(:,:,:) = 0.0d+00 l = 0 pmeta => list_meta do while(associated(pmeta)) if (pmeta%leaf) then l = l + 1 fields(l,1) = pmeta%id fields(l,2) = pmeta%level fields(l,3) = pmeta%coords(1) fields(l,4) = pmeta%coords(2) #if NDIMS == 3 fields(l,5) = pmeta%coords(3) #endif /* NDIMS == 3 */ bounds(l,1,1) = pmeta%xmin bounds(l,1,2) = pmeta%xmax bounds(l,2,1) = pmeta%ymin bounds(l,2,2) = pmeta%ymax #if NDIMS == 3 bounds(l,3,1) = pmeta%zmin bounds(l,3,2) = pmeta%zmax #endif /* NDIMS == 3 */ end if ! leaf pmeta => pmeta%next end do ! metablocks ! store metablock data ! write(fname,"(a,a)") "metablock_fields", trim(binary_file_suffix) call write_binary_xml(dname, fname, transfer(fields, [ 0_1 ]), & hash_type, dbytes, ddigest, cbytes, cdigest) call write_attribute_xml(lun, "fields", fname, hash_str, & dbytes, ddigest, cbytes, cdigest) write(fname,"(a,a)") "metablock_bounds", trim(binary_file_suffix) call write_binary_xml(dname, fname, transfer(bounds, [ 0_1 ]), & hash_type, dbytes, ddigest, cbytes, cdigest) call write_attribute_xml(lun, "bounds", fname, hash_str, & dbytes, ddigest, cbytes, cdigest) if (allocated(fields)) deallocate(fields) if (allocated(bounds)) deallocate(bounds) else call print_message(loc, "Cannot allocate space for metablocks!") status = 1001 return end if ! allocation write(lun,"(a)") '' write(lun,"(a)") '' close(lun) end if ! meta data file is stored only on the master process ! prepare and store data block info ! write(fname,fmt) trim(dname), "datablocks", nproc, "xml" open(newunit = lun, file = fname, status = 'replace') write(lun,"(a)") "" write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "nprocs" , nprocs) call write_attribute_xml(lun, "nproc" , nproc) call write_attribute_xml(lun, "ndims" , NDIMS) call write_attribute_xml(lun, "ncells" , ncells) call write_attribute_xml(lun, "nghosts" , nghosts) call write_attribute_xml(lun, "bcells" , nn) call write_attribute_xml(lun, "nvars" , nv) call write_attribute_xml(lun, "dblocks" , nd) call write_attribute_xml(lun, "variables", trim(vars)) write(lun,"(a)") '' write(lun,"(a)") '' if (nd > 0) then #if NDIMS == 3 allocate(ids(nd), array(nd,nn,nn,nn), stat = status) #else /* NDIMS == 3 */ allocate(ids(nd), array(nd,nn,nn, 1), stat = status) #endif /* NDIMS == 3 */ if (status == 0) then l = 0 pdata => list_data do while(associated(pdata)) l = l + 1 ids(l) = pdata%meta%id pdata => pdata%next end do ! data blocks write(fname,"(a,'_',a,'_',i6.6,a)") "datablock", "ids", & nproc, trim(binary_file_suffix) call write_binary_xml(dname, fname, transfer(ids, [ 0_1 ]), & hash_type, dbytes, ddigest, cbytes, cdigest) call write_attribute_xml(lun, "ids", fname, hash_str, & dbytes, ddigest, cbytes, cdigest) do p = 1, nv l = 0 pdata => list_data do while(associated(pdata)) l = l + 1 array(l,:,:,:) = pdata%q(p,:,:,:) pdata => pdata%next end do ! data blocks write(fname,"(a,'_',a,'_',i6.6,a)") "datablock", trim(pvars(p)), & nproc, trim(binary_file_suffix) call write_binary_xml(dname, fname, & transfer(array(:,:,:,:), [ 0_1 ]), & hash_type, dbytes, ddigest, cbytes, cdigest) call write_attribute_xml(lun, pvars(p), fname, hash_str, & dbytes, ddigest, cbytes, cdigest) end do if (allocated(ids)) deallocate(ids) if (allocated(array)) deallocate(array) else call print_message(loc, "Cannot allocate space for datablocks!") status = 1001 return end if ! allocation end if write(lun,"(a)") '' write(lun,"(a)") '' close(lun) !------------------------------------------------------------------------------- ! end subroutine write_snapshot_xml ! !=============================================================================== ! ! subroutine WRITE_ATTRIBUTE_XML_STRING: ! ------------------------------------- ! ! Subroutine writes a string attribute in XML format to specified ! file handler. ! ! Arguments: ! ! lun - the file handler to write to; ! aname - the name of attribute; ! avalue - the value of attribute; ! !=============================================================================== ! subroutine write_attribute_xml_string(lun, aname, avalue) implicit none ! input and output arguments ! integer , intent(in) :: lun character(len=*), intent(in) :: aname character(len=*), intent(in) :: avalue ! local parameters ! character(len=*), parameter :: afmt = "('',a,'')" ! !------------------------------------------------------------------------------- ! write(lun,afmt) "string", trim(adjustl(aname)), trim(adjustl(avalue)) !------------------------------------------------------------------------------- ! end subroutine write_attribute_xml_string ! !=============================================================================== ! ! subroutine WRITE_ATTRIBUTE_XML_INTEGER: ! -------------------------------------- ! ! Subroutine writes an integer attribute in XML format to specified ! file handler. ! ! Arguments: ! ! lun - the file handler to write to; ! aname - the name of attribute; ! avalue - the value of attribute; ! !=============================================================================== ! subroutine write_attribute_xml_integer(lun, aname, avalue) implicit none ! input and output arguments ! integer , intent(in) :: lun character(len=*), intent(in) :: aname integer(kind=4) , intent(in) :: avalue ! local variables ! character(len=32) :: svalue ! local parameters ! character(len=*), parameter :: afmt = "('',a,'')" ! !------------------------------------------------------------------------------- ! write(svalue,"(1i32)") avalue write(lun,afmt) "integer", trim(adjustl(aname)), trim(adjustl(svalue)) !------------------------------------------------------------------------------- ! end subroutine write_attribute_xml_integer ! !=============================================================================== ! ! subroutine WRITE_ATTRIBUTE_XML_DOUBLE: ! -------------------------------------- ! ! Subroutine writes a double precision attribute in XML format to specified ! file handler. ! ! Arguments: ! ! lun - the file handler to write to; ! aname - the name of attribute; ! avalue - the value of attribute; ! !=============================================================================== ! subroutine write_attribute_xml_double(lun, aname, avalue) implicit none ! input and output arguments ! integer , intent(in) :: lun character(len=*), intent(in) :: aname real(kind=8) , intent(in) :: avalue ! local variables ! character(len=32) :: svalue ! local parameters ! character(len=*), parameter :: afmt = "('',a,'')" ! !------------------------------------------------------------------------------- ! write(svalue,"(1es32.20)") avalue write(lun,afmt) "double", trim(adjustl(aname)), trim(adjustl(svalue)) !------------------------------------------------------------------------------- ! end subroutine write_attribute_xml_double ! !=============================================================================== ! ! subroutine WRITE_ATTRIBUTE_XML_FILE: ! ----------------------------------- ! ! Subroutine writes a file attribute in XML format to specified file handler. ! ! Arguments: ! ! lun - the file handler to write to; ! aname - the file attribute name; ! filename - the file name; ! digest_string - the digest type string; ! data_bytes - the file size in bytes; ! data_digest - the digest of the file content; ! compressed_bytes - the size of the compressed data in bytes; ! compressed_digest - the digest of the compressed data; ! !=============================================================================== ! subroutine write_attribute_xml_file(lun, aname, filename, & hash_string, data_bytes, data_digest, & compressed_bytes, compressed_digest) implicit none ! input and output arguments ! integer , intent(in) :: lun character(len=*) , intent(in) :: aname, filename, hash_string integer(kind=8) , intent(in) :: data_bytes, data_digest integer(kind=8) , optional, intent(in) :: compressed_bytes, compressed_digest ! local variables ! character(len=256) :: fname character(len=32) :: digest_string, bytes_string character(len=1024) :: string integer :: l ! !------------------------------------------------------------------------------- ! fname = filename string = ' 0) then write(bytes_string,"(1i32)") compressed_bytes string = trim(string) // & ' compression_format="' // trim(adjustl(cformat)) // '"' // & ' compressed_size="' // trim(adjustl(bytes_string)) // '"' if (present(compressed_digest)) then if (compressed_digest /= 0) then write(digest_string,"(1z0.16)") compressed_digest string = trim(string) // ' compressed_digest="' // & trim(adjustl(digest_string)) // '"' end if end if else l = index(fname, '.bin') + 3 fname = filename(1:l) end if end if string = trim(string) // '>' // trim(adjustl(fname)) // '' write(lun,'(a)') trim(adjustl(string)) !------------------------------------------------------------------------------- ! end subroutine write_attribute_xml_file ! !=============================================================================== ! ! subroutine WRITE_BINARY_XML: ! --------------------------- ! ! Subroutine writes the input array of bytes in a binary file with ! the provided path and name, and returns the digest of written data. ! ! Arguments: ! ! path, name - the path and name where the array should be written to; ! array - the array of bytes to be written; ! hash_type - the type of digest to hash the data; ! array_bytes - the size of the array in bytes; ! array_digest - the digest of the input array; ! compressed_bytes - the size of the compressed array in bytes; ! compressed_digest - the digest of the compressed array; ! !=============================================================================== ! subroutine write_binary_xml(path, name, array, hash_type, & array_bytes, array_digest, & compressed_bytes, compressed_digest) use compression, only : get_compression, compress use hash , only : get_hash implicit none ! input and output arguments ! character(len=*) , intent(in) :: path, name integer(kind=1), dimension(:), intent(in) :: array integer , intent(in) :: hash_type integer(kind=8), optional , intent(out) :: array_bytes, compressed_bytes integer(kind=8), optional , intent(out) :: array_digest, compressed_digest ! local variables ! character(len=512) :: fname integer :: lun = 123 logical :: written integer :: l, status ! compression buffer ! integer(kind=1), dimension(:), allocatable :: buffer ! !------------------------------------------------------------------------------- ! status = 0 written = .false. array_bytes = size(array, kind=8) if (present(array_digest)) array_digest = get_hash(array, hash_type) write(fname,"(a,'/',a)") trim(path), trim(name) ! try to compress array and write it ! if (present(compressed_bytes) .and. get_compression() > 0) then allocate(buffer(array_bytes), stat = status) if (status == 0) then call compress(array, buffer, compressed_bytes) if (compressed_bytes > 0) then open(newunit = lun, file = fname, form = 'unformatted', & access = 'stream', status = 'replace') write(lun) buffer(1:compressed_bytes) close(lun) written = .true. if (present(compressed_digest)) & compressed_digest = get_hash(buffer(1:compressed_bytes), hash_type) end if deallocate(buffer) end if end if ! compression failed of no compression is used, so writhe the uncompressed array ! if (.not. written) then l = index(fname, '.bin') + 3 open(newunit = lun, file = fname(:l), form = 'unformatted', & access = 'stream', status = 'replace') write(lun) array close(lun) end if !------------------------------------------------------------------------------- ! end subroutine write_binary_xml #ifdef HDF5 ! !=============================================================================== ! ! subroutine READ_SNAPSHOT_PARAMETER_STRING_H5: ! -------------------------------------------- ! ! Subroutine reads a string parameter from the restart snapshot. ! ! Arguments: ! ! pname - the parameter name; ! pvalue - the parameter value; ! iret - the success flag (the success is 0, failure otherwise); ! !=============================================================================== ! subroutine read_snapshot_parameter_string_h5(pname, pvalue, iret) use hdf5 , only : hid_t use hdf5 , only : H5F_ACC_RDONLY_F use hdf5 , only : H5T_NATIVE_CHARACTER use hdf5 , only : hid_t, hsize_t, size_t use hdf5 , only : h5fopen_f, h5fclose_f use hdf5 , only : h5gopen_f, h5gclose_f use hdf5 , only : h5aexists_by_name_f use hdf5 , only : h5tcopy_f, h5tset_size_f use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f use helpers , only : print_message use mpitools , only : nproc use parameters, only : get_parameter implicit none character(len=*), intent(in) :: pname character(len=*), intent(inout) :: pvalue integer , intent(inout) :: iret logical :: info character(len=255) :: rname integer :: np integer(hid_t) :: fid, gid, tid, aid integer(size_t) :: aln integer(hsize_t) :: am(1) = 1 character(len=*), parameter :: loc = & 'IO::read_snapshot_parameter_string_h5()' !------------------------------------------------------------------------------- ! iret = 0 ! generate the filename of the restart snapshot ! info = .false. np = nproc + 1 do while (.not. info .and. np >= 0) np = np - 1 write (rname, "(a,'r',i6.6,'_',i5.5,'.h5')") trim(respath), nrest, np inquire(file = rname, exist = info) end do ! procees if file exists ! if (info) then call h5fopen_f(rname, H5F_ACC_RDONLY_F, fid, iret) if (iret >= 0) then call h5gopen_f(fid, 'attributes', gid, iret) if (iret >= 0) then call h5aexists_by_name_f(gid, '.', trim(pname), info, iret) if (info .and. iret >= 0) then call h5aopen_by_name_f(gid, '.', trim(pname), aid, iret) if (iret >= 0) then aln = len(pvalue) call h5tcopy_f(H5T_NATIVE_CHARACTER, tid, iret) call h5tset_size_f(tid, aln, iret) call h5aread_f(aid, tid, pvalue, am(:), iret) call h5aclose_f(aid, iret) end if end if call h5gclose_f(gid, iret) end if call h5fclose_f(fid, iret) end if else call print_message(loc, "Snapshot " // trim(rname) // & " file does not exist!") iret = 1 end if !------------------------------------------------------------------------------- ! end subroutine read_snapshot_parameter_string_h5 ! !=============================================================================== ! ! subroutine READ_SNAPSHOT_PARAMETER_INTEGER_H5: ! --------------------------------------------- ! ! Subroutine reads an integer parameter from the restart snapshot. ! ! Arguments: ! ! pname - the parameter name; ! pvalue - the parameter value; ! iret - the success flag (the success is 0, failure otherwise); ! !=============================================================================== ! subroutine read_snapshot_parameter_integer_h5(pname, pvalue, iret) use hdf5 , only : hid_t use hdf5 , only : H5F_ACC_RDONLY_F use hdf5 , only : H5T_NATIVE_INTEGER use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5fopen_f, h5fclose_f use hdf5 , only : h5gopen_f, h5gclose_f use hdf5 , only : h5aexists_by_name_f use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f use helpers , only : print_message use mpitools , only : nproc use parameters, only : get_parameter implicit none character(len=*), intent(in) :: pname integer , intent(inout) :: pvalue integer , intent(inout) :: iret logical :: info character(len=255) :: rname integer :: np integer(hid_t) :: fid, gid, aid integer(hsize_t) :: am(1) = 1 character(len=*), parameter :: loc = & 'IO::read_snapshot_parameter_integer_h5()' !------------------------------------------------------------------------------- ! ! reset the success flag ! iret = 0 ! generate the filename of the restart snapshot ! info = .false. np = nproc + 1 do while (.not. info .and. np >= 0) np = np - 1 write (rname, "(a,'r',i6.6,'_',i5.5,'.h5')") trim(respath), nrest, np inquire(file = rname, exist = info) end do ! procees if file exists ! if (info) then call h5fopen_f(rname, H5F_ACC_RDONLY_F, fid, iret) if (iret >= 0) then call h5gopen_f(fid, 'attributes', gid, iret) if (iret >= 0) then call h5aexists_by_name_f(gid, '.', trim(pname), info, iret) if (info .and. iret >= 0) then call h5aopen_by_name_f(gid, '.', trim(pname), aid, iret) if (iret >= 0) then call h5aread_f(aid, H5T_NATIVE_INTEGER, pvalue, am(:), iret) call h5aclose_f(aid, iret) end if end if call h5gclose_f(gid, iret) end if call h5fclose_f(fid, iret) end if else call print_message(loc, "Snapshot " // trim(rname) // & " file does not exist!") iret = 1 end if !------------------------------------------------------------------------------- ! end subroutine read_snapshot_parameter_integer_h5 ! !=============================================================================== ! ! subroutine READ_SNAPSHOT_PARAMETER_DOUBLE_H5: ! -------------------------------------------- ! ! Subroutine reads a double precision real parameter from the restart ! snapshot. ! ! Arguments: ! ! pname - the parameter name; ! pvalue - the parameter value; ! iret - the success flag (the success is 0, failure otherwise); ! !=============================================================================== ! subroutine read_snapshot_parameter_double_h5(pname, pvalue, iret) ! import external procedures ! use hdf5 , only : hid_t use hdf5 , only : H5F_ACC_RDONLY_F use hdf5 , only : H5T_NATIVE_DOUBLE use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5fopen_f, h5fclose_f use hdf5 , only : h5gopen_f, h5gclose_f use hdf5 , only : h5aexists_by_name_f use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f use iso_fortran_env, only : error_unit use mpitools , only : nproc use parameters , only : get_parameter ! local variables are not implicit by default ! implicit none ! subroutine arguments ! character(len=*), intent(in) :: pname real(kind=8) , intent(inout) :: pvalue integer , intent(inout) :: iret ! local variables ! logical :: info character(len=255) :: rname integer :: np integer(hid_t) :: fid, gid, aid integer(hsize_t) :: am(1) = 1 ! local parameters ! character(len=*), parameter :: loc = & 'IO::read_snapshot_parameter_double_h5()' ! !------------------------------------------------------------------------------- ! ! reset the success flag ! iret = 0 ! generate the filename of the restart snapshot ! info = .false. np = nproc + 1 do while (.not. info .and. np >= 0) np = np - 1 write (rname, "(a,'r',i6.6,'_',i5.5,'.h5')") trim(respath), nrest, np inquire(file = rname, exist = info) end do ! procees if file exists ! if (info) then call h5fopen_f(rname, H5F_ACC_RDONLY_F, fid, iret) if (iret >= 0) then call h5gopen_f(fid, 'attributes', gid, iret) if (iret >= 0) then call h5aexists_by_name_f(gid, '.', trim(pname), info, iret) if (info .and. iret >= 0) then call h5aopen_by_name_f(gid, '.', trim(pname), aid, iret) if (iret >= 0) then call h5aread_f(aid, H5T_NATIVE_DOUBLE, pvalue, am(:), iret) call h5aclose_f(aid, iret) end if end if call h5gclose_f(gid, iret) end if call h5fclose_f(fid, iret) end if else write(error_unit,"('[', a, ']: ', a)") trim(loc) & , "Snapshot " // trim(rname) // " file does not exist!" iret = 1 end if !------------------------------------------------------------------------------- ! end subroutine read_snapshot_parameter_double_h5 ! !=============================================================================== ! ! subroutine READ_RESTART_SNAPSHOT_H5: ! ----------------------------------- ! ! Subroutine reads restart snapshot, i.e. parameters, meta and data blocks ! stored in the HDF5 format restart files and reconstructs the data structure ! in order to resume a terminated job. ! ! Arguments: ! ! iret - the return flag to inform if subroutine succeeded or failed; ! !=============================================================================== ! subroutine read_restart_snapshot_h5(iret) ! import external procedures and variables ! use blocks , only : change_blocks_process use forcing , only : einj use hdf5 , only : hid_t use hdf5 , only : H5F_ACC_RDONLY_F use hdf5 , only : h5fis_hdf5_f, h5fopen_f, h5fclose_f use hdf5 , only : h5gopen_f, h5gclose_f use iso_fortran_env, only : error_unit #ifdef MPI use mesh , only : redistribute_blocks #endif /* MPI */ use mpitools , only : nprocs, npmax, nproc ! local variables are not implicit by default ! implicit none ! input and output arguments ! integer, intent(out) :: iret ! local variables ! character(len=255) :: fl, msg integer(hid_t) :: fid, gid integer :: err, lfile logical :: info real(kind=8) :: deinj ! local parameters ! character(len=*), parameter :: loc = 'IO::read_restart_snapshot_h5()' ! !------------------------------------------------------------------------------- ! ! initialize success flag ! iret = 0 !! 1. RESTORE PARAMETERS AND META BLOCKS FROM THE FIRST FILE !! ! prepare the filename using the current process number; in case the file does ! not exist decrease it until the file corresponding to lower process number ! is found; ! write (fl, "(a,'r',i6.6,'_',i5.5,'.h5')") trim(respath), nrest, 0 inquire(file = fl, exist = info) ! quit, if file does not exist ! if (.not. info) then iret = 121 msg = "File " // trim(fl) // " does not exist!" else ! file does exist ! check if this file is in the HDF5 format ! call h5fis_hdf5_f(fl, info, err) ! quit, if the format verification failed or file is not in HDF5 format ! if (err < 0 .or. .not. info) then iret = 122 if (err < 0) msg = "Cannot check the file format!" if (.not. info) msg = "File " // trim(fl) // " is not an HDF5 file!" else ! file is HDF5 ! open the HDF5 file ! call h5fopen_f(fl, H5F_ACC_RDONLY_F, fid, err) ! quit, if file could not be opened ! if (err < 0) then iret = 123 msg = "Cannot open file: " // trim(fl) else ! file is opened ! read global attributes ! call read_attributes_h5(fid) ! read meta blocks and recreate the meta block hierarchy ! call read_metablocks_h5(fid) ! close the file ! call h5fclose_f(fid, err) ! quit, if file could not be closed ! if (err > 0) then iret = 124 msg = "Cannot close file: " // trim(fl) end if end if ! file is opened end if ! file is HDF5 end if ! file does exist !! 1. RESTORE DATA BLOCKS !! ! separate data blocks reading into two cases, when the number of processors is ! larger or equal to the number of files, and when we have less processors than ! files ! ! first, read data blocks to processes which have corresponding restart file ! if (nproc < nfiles) then ! prepare the filename ! write (fl, "(a,'r',i6.6,'_',i5.5,'.h5')") trim(respath), nrest, nproc ! check if the HDF5 file exists ! inquire(file = fl, exist = info) ! quit, if file does not exist ! if (.not. info) then iret = 121 msg = "File " // trim(fl) // " does not exist!" else ! file does exist ! check if this file is in the HDF5 format ! call h5fis_hdf5_f(fl, info, err) ! quit, if the format verification failed or file is not in HDF5 format ! if (err < 0 .or. .not. info) then iret = 122 if (err < 0) msg = "Cannot check the file format!" if (.not. info) msg = "File " // trim(fl) // " is not an HDF5 file!" else ! file is HDF5 ! open the HDF5 file ! call h5fopen_f(fl, H5F_ACC_RDONLY_F, fid, err) ! quit, if file could not be opened ! if (err < 0) then iret = 123 msg = "Cannot open file: " // trim(fl) else ! file is opened ! restore injected energy ! call h5gopen_f(fid, 'attributes', gid, err) if (err >= 0) then call read_attribute(gid, 'einj', einj) call h5gclose_f(gid, err) if (err /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close the group!" end if else write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open the group!" end if ! read data blocks ! call read_datablocks_h5(fid) ! close the file ! call h5fclose_f(fid, err) ! quit, if file could not be closed ! if (err > 0) then iret = 124 msg = "Cannot close file: " // trim(fl) end if end if ! file is opened end if ! file is HDF5 end if ! file exists end if ! nproc < nfiles ! if there are more files than processes, read the remaining files by ! the last process and redistribute blocks after each processed file, ! otherwise only redistribute blocks ! if (nprocs < nfiles) then ! iterate over remaining files and read one by one to the last block ! do lfile = nprocs, nfiles - 1 ! switch meta blocks from the read file to belong to the reading process ! call change_blocks_process(lfile, npmax) ! read the remaining files by the last process only ! if (nproc == npmax) then ! prepare the filename ! write (fl, "(a,'r',i6.6,'_',i5.5,'.h5')") trim(respath), nrest, lfile ! check if the HDF5 file exists ! inquire(file = fl, exist = info) ! quit, if file does not exist ! if (.not. info) then iret = 121 msg = "File " // trim(fl) // " does not exist!" else ! file does exist ! check if this file is in the HDF5 format ! call h5fis_hdf5_f(fl, info, err) ! quit, if the format verification failed or file is not in HDF5 format ! if (err < 0 .or. .not. info) then iret = 122 if (err < 0) msg = "Cannot check the file format!" if (.not. info) msg = "File " // trim(fl) // & " is not an HDF5 file!" else ! file is HDF5 ! open the HDF5 file ! call h5fopen_f(fl, H5F_ACC_RDONLY_F, fid, err) ! quit, if file could not be opened ! if (err < 0) then iret = 123 msg = "Cannot open file: " // trim(fl) else ! file is opened ! restore injected energy ! call h5gopen_f(fid, 'attributes', gid, err) if (err >= 0) then call read_attribute(gid, 'einj', deinj) einj = einj + deinj call h5gclose_f(gid, err) if (err /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close the group!" end if else write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open the group!" end if ! read data blocks ! call read_datablocks_h5(fid) ! close the file ! call h5fclose_f(fid, err) ! quit, if file could not be closed ! if (err > 0) then iret = 124 msg = "Cannot close file: " // trim(fl) end if end if ! file is opened end if ! file is HDF5 end if ! file exists end if ! nproc == npmax #ifdef MPI ! redistribute blocks between processors ! call redistribute_blocks() #endif /* MPI */ end do ! lfile = nprocs, nfiles - 1 else ! nprocs < nfiles #ifdef MPI ! redistribute blocks between processors ! call redistribute_blocks() #endif /* MPI */ end if ! nprocs < nfiles ! deallocate the array of block pointers ! if (allocated(block_array)) deallocate(block_array) ! if there was any problem, print the message ! if (iret > 0) write(error_unit,"('[',a,']: ',a)") trim(loc), trim(msg) !------------------------------------------------------------------------------- ! end subroutine read_restart_snapshot_h5 ! !=============================================================================== ! ! subroutine STORE_RESTART_SNAPSHOT_H5: ! ------------------------------------ ! ! Subroutine stores restart snapshots in the HDF5 format. ! ! Arguments: ! ! problem - the problem's name; ! nrun - the snapshot number; ! status - the subroutine call status; ! !=============================================================================== ! subroutine store_restart_snapshot_h5(problem, nrun, status) use hdf5 use helpers , only : print_message use mpitools, only : nproc implicit none character(len=*), intent(in) :: problem integer , intent(in) :: nrun integer , intent(out) :: status character(len=255) :: fname integer(hid_t) :: file_id character(len=*), parameter :: loc = 'IO::store_restart_snapshot_h5()' !------------------------------------------------------------------------------- ! write(fname, "('r',i6.6,'_',i5.5,'.h5')") nrun, nproc call h5fcreate_f(fname, H5F_ACC_TRUNC_F, file_id, status) if (status < 0) then call print_message(loc, "Could not create file " // trim(fname)) return end if call store_attributes_h5(file_id, problem, .true., status) call store_metablocks_h5(file_id, status) call store_datablocks_h5(file_id, status) call h5fclose_f(file_id, status) if (status < 0) & call print_message(loc, "Could not close file " // trim(fname)) !------------------------------------------------------------------------------- ! end subroutine store_restart_snapshot_h5 ! !=============================================================================== ! ! subroutine STORE_SNAPSHOT_H5: ! ---------------------------- ! ! Subroutine stores the current simulation snapshot, i.e. parameters, ! coordinates and variables as a HDF5 file. ! ! Arguments: ! ! problem - the problem's name; ! status - the subroutine call status; ! !=============================================================================== ! subroutine store_snapshot_h5(problem, status) use hdf5 use helpers , only : print_message use mpitools, only : nproc implicit none character(len=*), intent(in) :: problem integer , intent(out) :: status character(len=255) :: fname integer(hid_t) :: file_id character(len=*), parameter :: loc = 'IO::store_snapshot_h5()' !------------------------------------------------------------------------------- ! write(fname,"('p',i6.6,'_',i5.5,'.h5')") isnap, nproc call h5fcreate_f(fname, H5F_ACC_TRUNC_F, file_id, status) if (status < 0) then call print_message(loc, "Could not create file " // trim(fname)) return end if call store_attributes_h5(file_id, problem, .false., status) call store_coordinates_h5(file_id, status) call store_variables_h5(file_id, status) call h5fclose_f(file_id, status) if (status < 0) & call print_message(loc, "Could not close file " // trim(fname)) !------------------------------------------------------------------------------- ! end subroutine store_snapshot_h5 ! !=============================================================================== ! ! subroutine STORE_ATTRIBUTES_H5: ! ------------------------------ ! ! Subroutine stores the global simulation attributes in the specific location. ! ! Arguments: ! ! loc_id - the HDF5 file identifier; ! problem - the problem's name; ! restart - the flag indicating to store attributes for restart snapshot; ! status - the subroutine call status; ! !=============================================================================== ! subroutine store_attributes_h5(loc_id, problem, restart, status) use hdf5 use blocks , only : get_mblocks, get_dblocks, get_nleafs use blocks , only : get_last_id, nchildren use coordinates, only : minlev, maxlev use coordinates, only : bcells, ncells, nghosts use coordinates, only : bdims => domain_base_dims use coordinates, only : xmin, xmax, ymin, ymax, zmin, zmax use coordinates, only : periodic use equations , only : eqsys, eos, adiabatic_index, csnd, cmax, nv use evolution , only : step, time, dt, dth, dte, cfl, glm_alpha, errs use evolution , only : atol, rtol, mrej, niterations, nrejections use forcing , only : nmodes, einj, fcoefs use helpers , only : print_message use mpitools , only : nprocs, nproc use random , only : gentype, nseeds, get_seeds use sources , only : viscosity, resistivity implicit none integer(hid_t) , intent(in) :: loc_id character(len=*), intent(in) :: problem logical , intent(in) :: restart integer , intent(out) :: status integer(hid_t) :: grp_id integer(hsize_t), dimension(2) :: dims = 1 integer(kind=8), dimension(:,:), allocatable :: seeds real(kind=8) , dimension(:,:), allocatable :: array character(len=*), parameter :: loc = 'IO::store_attributes_h5()' !------------------------------------------------------------------------------- ! call store_attribute_h5(loc_id, 'code', 'AMUN', status) call store_attribute_h5(loc_id, 'version', H5T_NATIVE_REAL, 1, 1.0, status) ! store generic modules ! call h5gcreate_f(loc_id, 'attributes', grp_id, status) if (status >= 0) then call store_attribute_h5(grp_id, "problem", trim(problem), status) call store_attribute_h5(grp_id, "eqsys" , trim(eqsys) , status) call store_attribute_h5(grp_id, 'eos' , trim(eos) , status) call store_attribute_h5(grp_id, 'nprocs', & H5T_NATIVE_INTEGER, 1, nprocs, status) call store_attribute_h5(grp_id, 'nproc', & H5T_NATIVE_INTEGER, 1, nproc, status) call store_attribute_h5(grp_id, 'nvars', & H5T_NATIVE_INTEGER, 1, nv, status) call store_attribute_h5(grp_id, 'ndims', & H5T_NATIVE_INTEGER, 1, NDIMS, status) call store_attribute_h5(grp_id, 'bdims', & H5T_NATIVE_INTEGER, 3, bdims, status) call store_attribute_h5(grp_id, 'minlev', & H5T_NATIVE_INTEGER, 1, minlev, status) call store_attribute_h5(grp_id, 'maxlev', & H5T_NATIVE_INTEGER, 1, maxlev, status) call store_attribute_h5(grp_id, 'ncells', & H5T_NATIVE_INTEGER, 1, ncells, status) call store_attribute_h5(grp_id, 'nghosts', & H5T_NATIVE_INTEGER, 1, nghosts, status) call store_attribute_h5(grp_id, 'bcells', & H5T_NATIVE_INTEGER, 1, bcells, status) call store_attribute_h5(grp_id, 'dblocks', & H5T_NATIVE_INTEGER, 1, get_dblocks(), status) call store_attribute_h5(grp_id, 'nleafs', & H5T_NATIVE_INTEGER, 1, get_nleafs(), status) call store_attribute_h5(grp_id, 'step', & H5T_NATIVE_INTEGER, 1, step, status) call store_attribute_h5(grp_id, 'isnap', & H5T_NATIVE_INTEGER, 1, isnap, status) call store_attribute_h5(grp_id, 'periodic', & H5T_NATIVE_INTEGER, 3, periodic, status) call store_attribute_h5(grp_id, 'adiabatic_index', & H5T_NATIVE_DOUBLE, 1, adiabatic_index, status) call store_attribute_h5(grp_id, 'sound_speed', & H5T_NATIVE_DOUBLE, 1, csnd, status) call store_attribute_h5(grp_id, 'viscosity', & H5T_NATIVE_DOUBLE, 1, viscosity, status) call store_attribute_h5(grp_id, 'resistivity', & H5T_NATIVE_DOUBLE, 1, resistivity, status) call store_attribute_h5(grp_id, 'xmin', & H5T_NATIVE_DOUBLE, 1, xmin, status) call store_attribute_h5(grp_id, 'xmax', & H5T_NATIVE_DOUBLE, 1, xmax, status) call store_attribute_h5(grp_id, 'ymin', & H5T_NATIVE_DOUBLE, 1, ymin, status) call store_attribute_h5(grp_id, 'ymax', & H5T_NATIVE_DOUBLE, 1, ymax, status) call store_attribute_h5(grp_id, 'zmin', & H5T_NATIVE_DOUBLE, 1, zmin, status) call store_attribute_h5(grp_id, 'zmax', & H5T_NATIVE_DOUBLE, 1, zmax, status) call store_attribute_h5(grp_id, 'time', & H5T_NATIVE_DOUBLE, 1, time, status) call store_attribute_h5(grp_id, 'dt' , & H5T_NATIVE_DOUBLE, 1, dt, status) call store_attribute_h5(grp_id, 'cfl' , & H5T_NATIVE_DOUBLE, 1, cfl, status) call store_attribute_h5(grp_id, 'glm_alpha', & H5T_NATIVE_DOUBLE, 1, glm_alpha, status) if (restart) then call store_attribute_h5(grp_id, 'rngtype', trim(gentype), status) call store_attribute_h5(grp_id, 'nchildren', & H5T_NATIVE_INTEGER, 1, nchildren, status) call store_attribute_h5(grp_id, 'mblocks', & H5T_NATIVE_INTEGER, 1, get_mblocks(), status) call store_attribute_h5(grp_id, 'last_id', & H5T_NATIVE_INTEGER, 1, get_last_id(), status) call store_attribute_h5(grp_id, 'maximum_rejections', & H5T_NATIVE_INTEGER, 1, mrej, status) call store_attribute_h5(grp_id, 'niterations', & H5T_NATIVE_INTEGER, 1, niterations, status) call store_attribute_h5(grp_id, 'nrejections', & H5T_NATIVE_INTEGER, 1, nrejections, status) call store_attribute_h5(grp_id, 'nmodes', & H5T_NATIVE_INTEGER, 1, nmodes, status) call store_attribute_h5(grp_id, 'nseeds', & H5T_NATIVE_INTEGER, 1, nseeds, status) call store_attribute_h5(grp_id, 'dth' , & H5T_NATIVE_DOUBLE, 1, dth, status) call store_attribute_h5(grp_id, 'dte' , & H5T_NATIVE_DOUBLE, 1, dte, status) call store_attribute_h5(grp_id, 'cmax', & H5T_NATIVE_DOUBLE, 1, cmax, status) call store_attribute_h5(grp_id, 'absolute_tolerance', & H5T_NATIVE_DOUBLE, 1, atol, status) call store_attribute_h5(grp_id, 'relative_tolerance', & H5T_NATIVE_DOUBLE, 1, rtol, status) call store_attribute_h5(grp_id, 'errs', & H5T_NATIVE_DOUBLE, 3, errs, status) call store_attribute_h5(grp_id, 'einj', & H5T_NATIVE_DOUBLE, 1, einj, status) end if call h5gclose_f(grp_id, status) if (status < 0) & call print_message(loc, "Could not close group 'attributes'!") else call print_message(loc, "Could not create group 'attributes'!") end if ! store forcing coefficients in a different group ! if (restart) then call h5gcreate_f(loc_id, 'forcing', grp_id, status) if (status >= 0) then if (nmodes > 0) then dims(1) = nmodes dims(2) = NDIMS allocate(array(nmodes,NDIMS), stat=status) if (status == 0) then array = real(fcoefs) call store_dataset_h5(grp_id, 'fcoefs_real', H5T_NATIVE_DOUBLE, & dims, array, status) if (status < 0) & call print_message(loc, "Could not store real(fcoefs)!") array = aimag(fcoefs) call store_dataset_h5(grp_id, 'fcoefs_imag', H5T_NATIVE_DOUBLE, & dims, array, status) if (status < 0) & call print_message(loc, "Could not store imag(fcoefs)!") deallocate(array, stat=status) if (status /= 0) & call print_message(loc, "Could not deallocate space for fcoefs!") else call print_message(loc, "Could not allocate space for fcoefs!") end if end if call h5gclose_f(grp_id, status) if (status < 0) & call print_message(loc, "Could not close group 'forcing'!") else call print_message(loc, "Could not create group 'forcing'!") end if ! store random seeds in a different group ! call h5gcreate_f(loc_id, 'random', grp_id, status) if (status >= 0) then if (nseeds > 0) then dims(1) = 4 dims(2) = nseeds allocate(seeds(4,nseeds), stat=status) if (status == 0) then call get_seeds(seeds) call store_dataset_h5(grp_id, 'seeds', H5T_STD_I64LE, & dims, seeds, status) if (status < 0) & call print_message(loc, "Could not store seeds!") deallocate(seeds, stat=status) if (status /= 0) & call print_message(loc, "Could not deallocate space for seeds!") else call print_message(loc, "Could not allocate space for seeds!") end if end if call h5gclose_f(grp_id, status) if (status < 0) & call print_message(loc, "Could not close group 'random'!") else call print_message(loc, "Could not create group 'random'!") end if end if !------------------------------------------------------------------------------- ! end subroutine store_attributes_h5 ! !=============================================================================== ! ! subroutine READ_ATTRIBUTES_H5: ! ----------------------------- ! ! Subroutine restores global attributes from an HDF5 file provided by its ! identifier. ! ! Arguments: ! ! fid - the HDF5 file identifier; ! !=============================================================================== ! subroutine read_attributes_h5(fid) ! import external procedures and variables ! use blocks , only : block_meta use blocks , only : append_metablock use blocks , only : set_last_id, get_last_id use blocks , only : get_mblocks, get_dblocks, get_nleafs use coordinates , only : ncells use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax use equations , only : cmax, cmax2 use evolution , only : step, time, dt, dth, dte use evolution , only : niterations, nrejections, errs use forcing , only : nmodes, fcoefs use hdf5 , only : hid_t use hdf5 , only : h5gopen_f, h5gclose_f use iso_fortran_env, only : error_unit use random , only : set_seeds, gentype ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer(hid_t), intent(in) :: fid ! local variables ! integer(hid_t) :: gid integer :: ierr, l integer :: lndims, lmblocks, lnleafs, llast_id integer :: lncells, lnproc, lnseeds, lnmodes integer :: status ! local pointers ! type(block_meta), pointer :: pmeta ! allocatable arrays ! integer(kind=8), dimension(:,:), allocatable :: seeds ! local parameters ! character(len=*), parameter :: loc = 'IO::read_attributes_h5()' ! !------------------------------------------------------------------------------- ! ! open the global attributes group ! call h5gopen_f(fid, 'attributes', gid, ierr) ! check if the group has been opened successfuly ! if (ierr < 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open the group!" return end if ! restore integer attributes ! call read_attribute(gid, 'ndims' , lndims ) call read_attribute(gid, 'nprocs' , nfiles ) call read_attribute(gid, 'nproc' , lnproc ) call read_attribute(gid, 'mblocks', lmblocks) call read_attribute(gid, 'nleafs' , lnleafs ) call read_attribute(gid, 'last_id', llast_id) call read_attribute(gid, 'ncells' , lncells ) call read_attribute(gid, 'nseeds' , lnseeds ) call read_attribute(gid, 'step' , step ) call read_attribute(gid, 'isnap' , isnap ) call read_attribute(gid, 'niterations', niterations) call read_attribute(gid, 'nrejections', nrejections) ! restore double precision attributes ! call read_attribute(gid, 'xmin', xmin) call read_attribute(gid, 'xmax', xmax) call read_attribute(gid, 'ymin', ymin) call read_attribute(gid, 'ymax', ymax) call read_attribute(gid, 'zmin', zmin) call read_attribute(gid, 'zmax', zmax) call read_attribute(gid, 'time', time) call read_attribute(gid, 'dt' , dt ) call read_attribute(gid, 'dth' , dth ) call read_attribute(gid, 'dte' , dte ) call read_attribute(gid, 'cmax', cmax) cmax2 = cmax * cmax call read_attribute(gid, 'errs(1)', errs(1)) call read_attribute(gid, 'errs(2)', errs(2)) call read_attribute(gid, 'errs(3)', errs(3)) ! check the number of dimensions ! if (lndims /= NDIMS) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "The number of dimensions does not match!" return end if ! check the block dimensions ! if (lncells /= ncells) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "The block dimensions do not match!" end if ! restore forcing coefficients ! call read_attribute(gid, 'nmodes', lnmodes) if (lnmodes == nmodes) then if (lnmodes > 0) call read_attribute(gid, 'fcoefs', fcoefs) else write(error_unit,"('[',a,']: ',a)") trim(loc) & , "The number of driving modes does not match!" end if ! restore seeds ! if (trim(gentype) == "same") then allocate(seeds(4,0:lnseeds-1)) call read_attribute(gid, 'seeds', seeds(:,:)) call set_seeds(lnseeds, seeds(:,:), .false.) deallocate(seeds) end if ! allocate all metablocks ! do l = 1, lmblocks call append_metablock(pmeta, status) end do ! check if the number of created metablocks is equal to lbmcloks ! if (lmblocks /= get_mblocks()) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Number of metablocks does not match!" end if ! allocate an array of pointers with the size llast_id ! allocate(block_array(llast_id)) ! set the last_id ! call set_last_id(llast_id) ! close the group ! call h5gclose_f(gid, ierr) ! check if the group has been closed successfuly ! if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close the group!" end if !------------------------------------------------------------------------------- ! end subroutine read_attributes_h5 ! !=============================================================================== ! ! subroutine STORE_METABLOCKS_H5: ! ------------------------------ ! ! Subroutine stores all meta blocks' data in the group 'metablocks'. ! ! Arguments: ! ! loc_id - the location in which store the datablocks; ! status - the subroutine call status; ! !=============================================================================== ! subroutine store_metablocks_h5(loc_id, status) use hdf5 use blocks , only : block_meta, list_meta use blocks , only : ns => nsides, nc => nchildren use blocks , only : get_last_id, get_mblocks use helpers, only : print_message implicit none integer(hid_t), intent(in) :: loc_id integer , intent(out) :: status type(block_meta), pointer :: pmeta integer(hid_t) :: grp_id integer :: nm, l, n, p, i, j #if NDIMS == 3 integer :: k #endif /* NDIMS == 3 */ integer(hsize_t), dimension(5) :: dims integer(kind=4), dimension(:,:) , allocatable :: fields integer(kind=4), dimension(:,:) , allocatable :: children #if NDIMS == 2 integer(kind=4), dimension(:,:,:,:) , allocatable :: edges integer(kind=4), dimension(:,:,:) , allocatable :: corners #endif /* NDIMS == 2 */ #if NDIMS == 3 integer(kind=4), dimension(:,:,:,:,:), allocatable :: faces integer(kind=4), dimension(:,:,:,:,:), allocatable :: edges integer(kind=4), dimension(:,:,:,:) , allocatable :: corners #endif /* NDIMS == 3 */ real(kind=8) , dimension(:,:,:) , allocatable :: bounds character(len=*), parameter :: loc = 'IO::store_metablocks_h5()' !------------------------------------------------------------------------------- ! call h5gcreate_f(loc_id, 'metablocks', grp_id, status) if (status >= 0) then nm = get_mblocks() if (nm > 0) then allocate(fields(nm,16), children(nm,nc), bounds(nm,3,2), & #if NDIMS == 3 faces(nm,NDIMS,ns,ns,ns), & edges(nm,NDIMS,ns,ns,ns), corners(nm,ns,ns,ns), & #else /* NDIMS == 3 */ edges(nm,NDIMS,ns,ns), corners(nm,ns,ns), & #endif /* NDIMS == 3 */ stat = status) if (status == 0) then fields(:,:) = -1 children(:,:) = -1 #if NDIMS == 3 faces(:,:,:,:,:) = -1 edges(:,:,:,:,:) = -1 corners(:,:,:,:) = -1 #else /* NDIMS == 3 */ edges(:,:,:,:) = -1 corners(:,:,:) = -1 #endif /* NDIMS == 3 */ bounds(:,:,:) = 0.0d+00 l = 0 pmeta => list_meta do while(associated(pmeta)) l = l + 1 fields(l, 1) = pmeta%id fields(l, 2) = pmeta%process fields(l, 3) = pmeta%level fields(l, 4) = pmeta%conf fields(l, 5) = pmeta%refine fields(l, 6) = pmeta%pos(1) fields(l, 7) = pmeta%pos(2) #if NDIMS == 3 fields(l, 8) = pmeta%pos(3) #endif /* NDIMS == 3 */ fields(l, 9) = pmeta%coords(1) fields(l,10) = pmeta%coords(2) #if NDIMS == 3 fields(l,11) = pmeta%coords(3) #endif /* NDIMS == 3 */ if (pmeta%leaf) fields(l,12) = 1 if (associated(pmeta%data) ) fields(l,13) = 1 if (associated(pmeta%parent)) fields(l,14) = pmeta%parent%id do p = 1, nc if (associated(pmeta%child(p)%ptr)) & children(l,p) = pmeta%child(p)%ptr%id end do #if NDIMS == 2 do j = 1, ns do i = 1, ns do n = 1, NDIMS if (associated(pmeta%edges(i,j,n)%ptr)) & edges(l,n,i,j) = pmeta%edges(i,j,n)%ptr%id end do ! NDIMS if (associated(pmeta%corners(i,j)%ptr)) & corners(l,i,j) = pmeta%corners(i,j)%ptr%id end do ! i = 1, ns end do ! j = 1, ns #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, ns do j = 1, ns do i = 1, ns do n = 1, NDIMS if (associated(pmeta%faces(i,j,k,n)%ptr)) & faces(l,n,i,j,k) = pmeta%faces(i,j,k,n)%ptr%id if (associated(pmeta%edges(i,j,k,n)%ptr)) & edges(l,n,i,j,k) = pmeta%edges(i,j,k,n)%ptr%id end do ! NDIMS if (associated(pmeta%corners(i,j,k)%ptr)) & corners(l,i,j,k) = pmeta%corners(i,j,k)%ptr%id end do ! i = 1, ns end do ! j = 1, ns end do ! k = 1, ns #endif /* NDIMS == 3 */ bounds(l,1,1) = pmeta%xmin bounds(l,1,2) = pmeta%xmax bounds(l,2,1) = pmeta%ymin bounds(l,2,2) = pmeta%ymax #if NDIMS == 3 bounds(l,3,1) = pmeta%zmin bounds(l,3,2) = pmeta%zmax #endif /* NDIMS == 3 */ pmeta => pmeta%next end do ! over all meta blocks dims(1:2) = [ nm, 16 ] call store_dataset_h5(grp_id, 'fields', & H5T_NATIVE_INTEGER, dims(1:2), fields, status) dims(1:2) = [ nm, nc ] call store_dataset_h5(grp_id, 'children', & H5T_NATIVE_INTEGER, dims(1:2), children, status) #if NDIMS == 3 dims(1:5) = [ nm, NDIMS, ns, ns, ns ] call store_dataset_h5(grp_id, 'faces', & H5T_NATIVE_INTEGER, dims(1:5), faces, status) dims(1:5) = [ nm, NDIMS, ns, ns, ns ] call store_dataset_h5(grp_id, 'edges', & H5T_NATIVE_INTEGER, dims(1:5), edges, status) dims(1:4) = [ nm, ns, ns, ns ] call store_dataset_h5(grp_id, 'corners', & H5T_NATIVE_INTEGER, dims(1:4), corners, status) #else /* NDIMS == 3 */ dims(1:4) = [ nm, NDIMS, ns, ns ] call store_dataset_h5(grp_id, 'edges', & H5T_NATIVE_INTEGER, dims(1:4), edges, status) dims(1:3) = [ nm, ns, ns ] call store_dataset_h5(grp_id, 'corners', & H5T_NATIVE_INTEGER, dims(1:3), corners, status) #endif /* NDIMS == 3 */ dims(1:3) = [ nm, 3, 2 ] call store_dataset_h5(grp_id, 'bounds', & H5T_NATIVE_DOUBLE, dims(1:3), bounds, status) else call print_message(loc, "Could not allocate space for metablocks!") end if #if NDIMS == 3 deallocate(fields, children, bounds, faces, edges, corners, stat=status) #else /* NDIMS == 3 */ deallocate(fields, children, bounds, edges, corners, stat=status) #endif /* NDIMS == 3 */ if (status /= 0) & call print_message(loc, "Could not deallocate space of metablocks!") end if ! meta blocks > 0 call h5gclose_f(grp_id, status) if (status < 0) & call print_message(loc, "Could not close group 'metablocks'!") else call print_message(loc, "Could not create group 'metablocks'!") end if !------------------------------------------------------------------------------- ! end subroutine store_metablocks_h5 ! !=============================================================================== ! ! subroutine READ_METABLOCKS_H5: ! ----------------------------- ! ! Subroutine restores all meta blocks with their complete fields from ! 'metablock' group in a provided restart file identifier. ! ! Arguments: ! ! fid - the HDF5 file identifier; ! !=============================================================================== ! subroutine read_metablocks_h5(fid) ! import procedures and variables from other modules ! use blocks , only : block_meta, list_meta use blocks , only : ndims, nchildren, nsides use blocks , only : get_mblocks use blocks , only : metablock_set_id, metablock_set_process use blocks , only : metablock_set_refinement use blocks , only : metablock_set_configuration use blocks , only : metablock_set_level, metablock_set_position use blocks , only : metablock_set_coordinates, metablock_set_bounds use blocks , only : metablock_set_leaf use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5gopen_f, h5gclose_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer(hid_t), intent(in) :: fid ! local variables ! integer(hid_t) :: gid integer(kind=4) :: i, j, l, p, n, ip #if NDIMS == 3 integer(kind=4) :: k #endif /* NDIMS == 3 */ integer :: err integer(hsize_t), dimension(1) :: am integer(hsize_t), dimension(2) :: dm, pm #if NDIMS == 2 integer(hsize_t), dimension(4) :: nm #endif /* NDIMS == 2 */ #if NDIMS == 3 integer(hsize_t), dimension(5) :: nm #endif /* NDIMS == 3 */ ! local allocatable arrays ! integer(kind=4), dimension(:) , allocatable :: par, dat integer(kind=4), dimension(:) , allocatable :: id, cpu, lev, cfg, ref, lea real (kind=8), dimension(:) , allocatable :: xmn, xmx, ymn, ymx, zmn, zmx integer(kind=4), dimension(:,:), allocatable :: chl, pos, cor #if NDIMS == 2 integer(kind=4), dimension(:,:,:,:) , allocatable :: edges integer(kind=4), dimension(:,:,:) , allocatable :: corners #endif /* NDIMS == 2 */ #if NDIMS == 3 integer(kind=4), dimension(:,:,:,:,:), allocatable :: faces integer(kind=4), dimension(:,:,:,:,:), allocatable :: edges integer(kind=4), dimension(:,:,:,:) , allocatable :: corners #endif /* NDIMS == 3 */ ! local pointers ! type(block_meta), pointer :: pmeta ! local parameters ! character(len=*), parameter :: loc = 'IO::read_metablocks_h5()' ! !------------------------------------------------------------------------------- ! ! open metablock group ! call h5gopen_f(fid, 'metablocks', gid, err) ! check if the group has been opened successfuly ! if (err >= 0) then ! prepate dimensions ! am(1) = get_mblocks() dm(1) = get_mblocks() dm(2) = nchildren pm(1) = get_mblocks() pm(2) = NDIMS nm(1) = get_mblocks() nm(2) = nsides nm(3) = nsides #if NDIMS == 2 nm(4) = ndims #endif /* NDIMS == 2 */ #if NDIMS == 3 nm(4) = nsides nm(5) = ndims #endif /* NDIMS == 3 */ ! allocate arrays to restore metablocks data ! allocate(id (am(1))) allocate(cpu(am(1))) allocate(lev(am(1))) allocate(par(am(1))) allocate(dat(am(1))) allocate(cfg(am(1))) allocate(ref(am(1))) allocate(lea(am(1))) allocate(xmn(am(1))) allocate(xmx(am(1))) allocate(ymn(am(1))) allocate(ymx(am(1))) allocate(zmn(am(1))) allocate(zmx(am(1))) allocate(chl(dm(1),dm(2))) allocate(pos(pm(1),pm(2))) allocate(cor(pm(1),pm(2))) #if NDIMS == 2 allocate(edges (nm(1),nm(2),nm(3),nm(4))) allocate(corners(nm(1),nm(2),nm(3))) #endif /* NDIMS == 2 */ #if NDIMS == 3 allocate(faces (nm(1),nm(2),nm(3),nm(4),nm(5))) allocate(edges (nm(1),nm(2),nm(3),nm(4),nm(5))) allocate(corners(nm(1),nm(2),nm(3),nm(4))) #endif /* NDIMS == 3 */ ! reset vectors ! par(:) = -1 dat(:) = -1 lea(:) = -1 chl(:,:) = -1 #if NDIMS == 2 edges(:,:,:,:) = -1 corners(:,:,:) = -1 #endif /* NDIMS == 2 */ #if NDIMS == 3 faces(:,:,:,:,:) = -1 edges(:,:,:,:,:) = -1 corners(:,:,:,:) = -1 #endif /* NDIMS == 3 */ ! read metablock fields from the HDF5 file ! call read_array(gid, 'id' , am(:), id (:)) call read_array(gid, 'cpu' , am(:), cpu(:)) call read_array(gid, 'level' , am(:), lev(:)) call read_array(gid, 'config' , am(:), cfg(:)) call read_array(gid, 'refine' , am(:), ref(:)) call read_array(gid, 'leaf' , am(:), lea(:)) call read_array(gid, 'parent' , am(:), par(:)) call read_array(gid, 'xmin' , am(:), xmn(:)) call read_array(gid, 'xmax' , am(:), xmx(:)) call read_array(gid, 'ymin' , am(:), ymn(:)) call read_array(gid, 'ymax' , am(:), ymx(:)) call read_array(gid, 'zmin' , am(:), zmn(:)) call read_array(gid, 'zmax' , am(:), zmx(:)) call read_array(gid, 'pos' , pm(:), pos(:,:)) call read_array(gid, 'coord' , pm(:), cor(:,:)) call read_array(gid, 'child' , dm(:), chl(:,:)) #if NDIMS == 2 call read_array(gid, 'edges' , nm(1:4), edges(:,:,:,:)) call read_array(gid, 'corners', nm(1:3), corners(:,:,:)) #endif /* NDIMS == 2 */ #if NDIMS == 3 call read_array(gid, 'faces' , nm(1:5), faces(:,:,:,:,:)) call read_array(gid, 'edges' , nm(1:5), edges(:,:,:,:,:)) call read_array(gid, 'corners', nm(1:4), corners(:,:,:,:)) #endif /* NDIMS == 3 */ ! reset the block counter ! l = 0 ! associate pmeta with the first block on the meta block list ! pmeta => list_meta ! iterate over all meta blocks and restore their fields ! do while(associated(pmeta)) ! increase the block counter ! l = l + 1 ! restore meta block fields ! block_array(id(l))%ptr => pmeta call metablock_set_id (pmeta, id (l)) call metablock_set_process (pmeta, cpu(l)) call metablock_set_refinement (pmeta, ref(l)) call metablock_set_configuration(pmeta, cfg(l)) call metablock_set_level (pmeta, lev(l)) call metablock_set_position (pmeta, pos(l,:)) call metablock_set_coordinates (pmeta, cor(l,:)) call metablock_set_bounds (pmeta, xmn(l), xmx(l), ymn(l), ymx(l) & , zmn(l), zmx(l)) if (lea(l) == 1) call metablock_set_leaf(pmeta) ! associate pmeta with the next block on the list ! pmeta => pmeta%next end do ! over all meta blocks ! reset the block counter ! l = 0 ! associate pmeta with the first block on the meta block list ! pmeta => list_meta ! iterate over all meta blocks and restore their pointers ! do while(associated(pmeta)) ! increase the block counter ! l = l + 1 ! restore %parent pointer ! if (par(l) > 0) pmeta%parent => block_array(par(l))%ptr ! restore %child pointers ! do p = 1, nchildren if (chl(l,p) > 0) then pmeta%child(p)%ptr => block_array(chl(l,p))%ptr end if end do ! p = 1, nchildren ! restore %faces, %edges and %corners neighbor pointers ! #if NDIMS == 2 do i = 1, nsides do j = 1, nsides do n = 1, ndims ip = edges(l,i,j,n) if (ip > 0) pmeta%edges(i,j,n)%ptr => block_array(ip)%ptr end do ! n = 1, ndims ip = corners(l,i,j) if (ip > 0) pmeta%corners(i,j)%ptr => block_array(ip)%ptr end do ! i = 1, nsides end do ! j = 1, nsides #endif /* NDIMS == 2 */ #if NDIMS == 3 do i = 1, nsides do j = 1, nsides do k = 1, nsides do n = 1, ndims ip = faces(l,i,j,k,n) if (ip > 0) pmeta%faces(i,j,k,n)%ptr => block_array(ip)%ptr ip = edges(l,i,j,k,n) if (ip > 0) pmeta%edges(i,j,k,n)%ptr => block_array(ip)%ptr end do ! n = 1, ndims ip = corners(l,i,j,k) if (ip > 0) pmeta%corners(i,j,k)%ptr => block_array(ip)%ptr end do ! i = 1, nsides end do ! j = 1, nsides end do ! k = 1, nsides #endif /* NDIMS == 3 */ ! associate pmeta with the next block on the list ! pmeta => pmeta%next end do ! over all meta blocks ! deallocate allocatable arrays ! if (allocated(id) ) deallocate(id ) if (allocated(par)) deallocate(par) if (allocated(dat)) deallocate(dat) if (allocated(cpu)) deallocate(cpu) if (allocated(lev)) deallocate(lev) if (allocated(cfg)) deallocate(cfg) if (allocated(ref)) deallocate(ref) if (allocated(lea)) deallocate(lea) if (allocated(xmn)) deallocate(xmn) if (allocated(xmx)) deallocate(xmx) if (allocated(ymn)) deallocate(ymn) if (allocated(ymx)) deallocate(ymx) if (allocated(zmn)) deallocate(zmn) if (allocated(zmx)) deallocate(zmx) if (allocated(chl)) deallocate(chl) if (allocated(cor)) deallocate(cor) #if NDIMS == 3 if (allocated(faces)) deallocate(faces) #endif /* NDIMS == 3 */ if (allocated(edges)) deallocate(edges) if (allocated(corners)) deallocate(corners) ! close the group ! call h5gclose_f(gid, err) ! check if the group has been closed successfuly ! if (err > 0) then ! print error about the problem with closing the group ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close the metablock group!" end if else ! print error about the problem with opening the group ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open the metablock group!" end if !------------------------------------------------------------------------------- ! end subroutine read_metablocks_h5 ! !=============================================================================== ! ! subroutine RESTORE_DATABLOCKS_H5: ! -------------------------------- ! ! Subroutine reads all data blocks stored in the group 'datablocks' of ! the provided handler to the HDF5 restart file. ! ! Arguments: ! ! loc_id - the location in which store the datablocks; ! status - the subroutine call status; ! !=============================================================================== ! subroutine restore_datablocks_h5(loc_id, status) use hdf5 use blocks , only : block_meta, block_data use blocks , only : append_datablock, link_blocks use helpers , only : print_message implicit none integer(hid_t), intent(in) :: loc_id integer , intent(out) :: status type(block_data), pointer :: pdata integer(hid_t) :: grp_id, blk_id character(len=64) :: blk_name integer(kind=4) :: dblocks, l, id integer(hsize_t), dimension(4) :: pdims = 1 integer(hsize_t), dimension(5) :: cdims = 1 character(len=*), parameter :: loc = 'IO::restore_datablocks_h5()' !------------------------------------------------------------------------------- ! call h5gopen_f(loc_id, 'attributes', grp_id, status) if (status /= 0) then call print_message(loc, "Could not open group 'attributes'!") return end if call restore_attribute_h5(grp_id, 'dblocks', & H5T_NATIVE_INTEGER, 1, dblocks, status) call h5gclose_f(grp_id, status) if (status /= 0) then call print_message(loc, "Could not close group 'attributes'!") return end if if (dblocks == 0) return call h5gopen_f(loc_id, 'datablocks', grp_id, status) if (status /= 0) then call print_message(loc, "Could not open group 'datablocks'!") return end if do l = 1, dblocks write(blk_name, "('datablock_', i0)") l call append_datablock(pdata, status) if (status /= 0) then call print_message(loc, "Could not append new datablock!") go to 1000 end if call h5gopen_f(grp_id, blk_name, blk_id, status) if (status == 0) then call restore_attribute_h5(blk_id, 'meta', & H5T_NATIVE_INTEGER, 1, id, status) call link_blocks(block_array(id)%ptr, pdata) pdims = shape(pdata%q) cdims = shape(pdata%uu) call read_dataset_h5(blk_id, 'primitive_variables', & H5T_NATIVE_DOUBLE, pdims, pdata%q , status) call read_dataset_h5(blk_id, 'conservative_variables', & H5T_NATIVE_DOUBLE, cdims, pdata%uu, status) call h5gclose_f(blk_id, status) if (status /= 0) & call print_message(loc, & "Could not close group '" // trim(blk_name) // "'!") else call print_message(loc, & "Could not open group '" // trim(blk_name) // "'!") end if end do 1000 continue call h5gclose_f(grp_id, status) if (status /= 0) & call print_message(loc, "Could not close group 'datablocks'!") !------------------------------------------------------------------------------- ! end subroutine restore_datablocks_h5 ! !=============================================================================== ! ! subroutine STORE_DATABLOCKS_H5: ! ------------------------------ ! ! Subroutine stores all data blocks in the group 'datablocks'. ! ! Arguments: ! ! loc_id - the location in which store the datablocks; ! status - the subroutine call status; ! !=============================================================================== ! subroutine store_datablocks_h5(loc_id, status) use hdf5 use blocks , only : block_meta, block_data, list_data, get_dblocks use helpers, only : print_message implicit none integer(hid_t), intent(in) :: loc_id integer , intent(out) :: status type(block_data), pointer :: pdata character(len=64) :: blk_name integer(hid_t) :: grp_id, blk_id integer(kind=4) :: l integer(hsize_t), dimension(4) :: pdims = 1 integer(hsize_t), dimension(5) :: cdims = 1 character(len=*), parameter :: loc = 'IO::store_datablocks_h5()' !------------------------------------------------------------------------------- ! status = 0 call h5gcreate_f(loc_id, 'datablocks', grp_id, status) if (status /= 0) then call print_message(loc, "Could not create group 'datablocks'!") end if if (get_dblocks() > 0) then l = 0 pdata => list_data do while(associated(pdata)) l = l + 1 write(blk_name, "('datablock_', i0)") l call h5gcreate_f(grp_id, blk_name, blk_id, status) if (status == 0) then call store_attribute_h5(blk_id, 'meta', & H5T_NATIVE_INTEGER, 1, pdata%meta%id, status) pdims = shape(pdata%q) cdims = shape(pdata%uu) call store_dataset_h5(blk_id, 'primitive_variables', & H5T_NATIVE_DOUBLE, pdims, pdata%q, status) if (status /= 0) & call print_message(loc, & "Could not store the primitive variables in " // & trim(blk_name) // "!") call store_dataset_h5(blk_id, 'conservative_variables', & H5T_NATIVE_DOUBLE, cdims, pdata%uu, status) if (status /= 0) & call print_message(loc, & "Could not store the conservative variables in " // & trim(blk_name) // "!") call h5gclose_f(blk_id, status) if (status /= 0) & call print_message(loc, & "Could not close group for " // trim(blk_name) // "!") else call print_message(loc, & "Could not create group for " // trim(blk_name) // "!") end if pdata => pdata%next end do end if call h5gclose_f(grp_id, status) if (status /= 0) & call print_message(loc, "Could not close group 'datablocks'!") !------------------------------------------------------------------------------- ! end subroutine store_datablocks_h5 ! !=============================================================================== ! ! subroutine READ_DATABLOCKS_H5: ! ----------------------------- ! ! Subroutine reads all data blocks stored in the group 'datablocks' of ! the provided handler to the HDF5 restart file. ! ! Arguments: ! ! fid - the HDF5 file identifier; ! !=============================================================================== ! subroutine read_datablocks_h5(fid) ! import external procedures and variables ! use blocks , only : block_meta, block_data use blocks , only : append_datablock, link_blocks, nregs use coordinates , only : nn => bcells, ng => nghosts use equations , only : nv use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5gopen_f, h5gclose_f, h5lexists_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! subroutine variables ! integer(hid_t), intent(in) :: fid ! local pointers ! type(block_data), pointer :: pdata ! local variables ! character(len=16) :: bname integer(hid_t) :: gid, bid integer(hsize_t) :: l integer(kind=4) :: i integer :: dblocks, ncells, nghosts, nc, nb, ne, status ! local arrays ! integer(hsize_t), dimension(5) :: dm = 1 ! local allocatable arrays ! real(kind=8), dimension(:,:,:,:,:), allocatable :: uu ! local parameters ! character(len=*), parameter :: loc = 'IO::read_datablocks_h5()' ! !------------------------------------------------------------------------------- ! ! read the number of data blocks ! call h5gopen_f(fid, 'attributes', gid, status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open the attribute group!" return end if call read_attribute(gid, 'dblocks', dblocks) call read_attribute(gid, 'ncells' , ncells ) call read_attribute(gid, 'nghosts', nghosts) call h5gclose_f(gid, status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close the attribute group!" return end if ! restore data blocks only if there are any ! if (dblocks > 0) then ! open the group 'datablocks' ! call h5gopen_f(fid, 'datablocks', gid, status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open the data block group!" else ! prepare the number of cells in the stored blocks and corresponding indices ! nc = ncells + 2 * nghosts if (nghosts >= ng) then nb = 1 + (nghosts - ng) ne = nc - (nghosts - ng) else nb = 1 + (ng - nghosts) ne = nn - (ng - nghosts) end if ! fill out dimensions ! dm(1) = dblocks dm(2) = nv dm(3) = nc dm(4) = nc #if NDIMS == 3 dm(5) = nc #endif /* NDIMS == 3 */ ! allocate arrays for input data ! allocate(uu(3,dm(2),dm(3),dm(4),dm(5)), stat = status) ! check if allocation was successful ! if (status == 0) then ! iterate over data blocks ! do l = 1, dm(1) ! allocate and append to the end of the list a new datablock ! call append_datablock(pdata, status) ! create name for the current block ! write(bname, "('dblk_', i11.11)") l ! open the group for the current block fields ! call h5gopen_f(gid, bname, bid, status) if (status == 0) then ! get the id of the linked meta block ! call read_attribute(bid, 'meta', i) ! associate the corresponding meta block with the current data block ! call link_blocks(block_array(i)%ptr, pdata) ! read the data ! call read_array(bid, 'pvar' , dm(2:5), uu(1,:,:,:,:)) call read_array(bid, 'cvar0', dm(2:5), uu(2,:,:,:,:)) call read_array(bid, 'cvar1', dm(2:5), uu(3,:,:,:,:)) ! fill out the block arrays taking into account the change of nghosts ! if (nghosts >= ng) then #if NDIMS == 3 pdata%q (:,:,:,:) = uu(1,:,nb:ne,nb:ne,nb:ne) pdata%uu(:,:,:,:,1) = uu(2,:,nb:ne,nb:ne,nb:ne) if (nregs > 1) pdata%uu(:,:,:,:,2) = uu(3,:,nb:ne,nb:ne,nb:ne) #else /* NDIMS == 3 */ pdata%q (:,:,:,:) = uu(1,:,nb:ne,nb:ne, : ) pdata%uu(:,:,:,:,1) = uu(2,:,nb:ne,nb:ne, : ) if (nregs > 1) pdata%uu(:,:,:,:,2) = uu(3,:,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ else #if NDIMS == 3 pdata%q (:,nb:ne,nb:ne,nb:ne) = uu(1,:,:,:,:) pdata%uu(:,nb:ne,nb:ne,nb:ne,1) = uu(2,:,:,:,:) if (nregs > 1) pdata%uu(:,nb:ne,nb:ne,nb:ne,2) = uu(3,:,:,:,:) #else /* NDIMS == 3 */ pdata%q (:,nb:ne,nb:ne, : ) = uu(1,:,:,:,:) pdata%uu(:,nb:ne,nb:ne, : ,1) = uu(2,:,:,:,:) if (nregs > 1) pdata%uu(:,nb:ne,nb:ne, : ,2) = uu(3,:,:,:,:) #endif /* NDIMS == 3 */ end if ! close the current data block group ! call h5gclose_f(bid, status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close group '", trim(bname), "'!" end if else write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open group '", trim(bname), "'!" end if end do ! l = 1, dm(1) else ! allocate write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot allocate temporary array!" end if ! allocate ! deallocate allocatable arrays ! if (allocated(uu)) deallocate(uu) ! close the data block group ! call h5gclose_f(gid, status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close the data block group!" end if end if end if ! dblocks > 0 !------------------------------------------------------------------------------- ! end subroutine read_datablocks_h5 ! !=============================================================================== ! ! subroutine STORE_COORDINATES_H5: ! ------------------------------- ! ! Subroutine stores blocks' data such as their IDs, levels, coordinates, etc. ! in a specific location. ! ! Arguments: ! ! loc_id - the location in which store the coordinates; ! status - the subroutine call status; ! !=============================================================================== ! subroutine store_coordinates_h5(loc_id, status) use hdf5 use blocks , only : block_meta, block_data, list_data use blocks , only : get_dblocks use coordinates, only : toplev use coordinates, only : adx, ady #if NDIMS == 3 use coordinates, only : adz #endif /* NDIMS == 3 */ use helpers , only : print_message implicit none integer(hid_t), intent(in) :: loc_id integer , intent(out) :: status integer(hid_t) :: grp_id integer :: n type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata integer(hsize_t), dimension(1) :: am = 1, im = 1 integer(hsize_t), dimension(2) :: cm = 1 integer(hsize_t), dimension(3) :: bm = 1 integer(kind=4), dimension(:) , allocatable :: ids, levels integer(kind=4), dimension(:,:) , allocatable :: coords real (kind=8), dimension(:,:,:), allocatable :: bounds character(len=*), parameter :: loc = 'IO::store_coordinates_h5()' !------------------------------------------------------------------------------- ! status = 0 call h5gcreate_f(loc_id, 'coordinates', grp_id, status) if (status /= 0) then call print_message(loc, "Could not create group 'coordinates'!") return end if am(1) = toplev call store_dataset_h5(grp_id, 'dx', H5T_NATIVE_DOUBLE, am, adx, status) call store_dataset_h5(grp_id, 'dy', H5T_NATIVE_DOUBLE, am, ady, status) #if NDIMS == 3 call store_dataset_h5(grp_id, 'dz', H5T_NATIVE_DOUBLE, am, adz, status) #endif /* NDIMS == 3 */ if (get_dblocks() > 0) then n = get_dblocks() im(1) = n cm(1) = NDIMS cm(2) = n bm(1) = NDIMS bm(2) = 2 bm(3) = n allocate(ids(n), levels(n), coords(NDIMS,n), & bounds(NDIMS,2,n), stat=status) if (status /= 0) then call print_message(loc, "Could not allocate space for coordinates!") else n = 0 pdata => list_data do while(associated(pdata)) pmeta => pdata%meta n = n + 1 ids(n) = pmeta%id levels(n) = pmeta%level coords(:,n) = pmeta%coords(:) bounds(1,:,n) = [ pmeta%xmin, pmeta%xmax ] bounds(2,:,n) = [ pmeta%ymin, pmeta%ymax ] #if NDIMS == 3 bounds(3,:,n) = [ pmeta%zmin, pmeta%zmax ] #endif /* NDIMS == 3 */ pdata => pdata%next end do call store_dataset_h5(grp_id, 'ids', & H5T_NATIVE_INTEGER, im, ids, status) call store_dataset_h5(grp_id, 'levels', & H5T_NATIVE_INTEGER, im, levels, status) call store_dataset_h5(grp_id, 'coords', & H5T_NATIVE_INTEGER, cm, coords, status) call store_dataset_h5(grp_id, 'bounds', & H5T_NATIVE_DOUBLE, bm, bounds, status) deallocate(ids, levels, coords, bounds, stat=status) if (status > 0) & call print_message(loc, "Could not deallocate the coordinate space!") end if end if call h5gclose_f(grp_id, status) if (status < 0) & call print_message(loc, "Could not close group 'coordinates'!") !------------------------------------------------------------------------------- ! end subroutine store_coordinates_h5 ! !=============================================================================== ! ! subroutine STORE_VARIABLES_H5: ! ----------------------------- ! ! Subroutine stores primitive variables in a specific location. ! ! Arguments: ! ! loc_id - the location in which store the variables; ! status - the subroutine call status; ! !=============================================================================== ! subroutine store_variables_h5(loc_id, status) use hdf5 use blocks , only : block_data, list_data use blocks , only : get_dblocks use coordinates , only : bcells use equations , only : nv, pvars use helpers , only : print_message use iso_c_binding, only : c_loc implicit none integer(hid_t), intent(in) :: loc_id integer , intent(out) :: status integer(hid_t) :: grp_id integer :: n, p integer(hsize_t), dimension(4) :: dims = 1 type(block_data), pointer :: pdata real(kind=8), dimension(:,:,:,:), allocatable, target :: array character(len=*), parameter :: loc = 'IO::store_variables_h5()' !------------------------------------------------------------------------------- ! status = 0 call h5gcreate_f(loc_id, 'variables', grp_id, status) if (status /= 0) then call print_message(loc, "Could not create group 'variables'!") return end if if (get_dblocks() > 0) then dims(1:NDIMS) = bcells dims(4) = get_dblocks() allocate(array(dims(1),dims(2),dims(3),dims(4)), stat=status) if (status /= 0) then call print_message(loc, "Could not allocate space for variables!") else do p = 1, nv n = 0 pdata => list_data do while(associated(pdata)) n = n + 1 array(:,:,:,n) = pdata%q(p,:,:,:) pdata => pdata%next end do call store_dataset_h5(grp_id, trim(pvars(p)), & H5T_NATIVE_DOUBLE, dims, array, status) end do deallocate(array, stat=status) if (status /= 0) & call print_message(loc, "Could not deallocate the variable space!") end if end if call h5gclose_f(grp_id, status) if (status /= 0) & call print_message(loc, "Could not close group 'variables'!") !------------------------------------------------------------------------------- ! end subroutine store_variables_h5 ! !=============================================================================== ! ! subroutine RESTORE_ATTRIBUTE_STRING_H5: ! -------------------------------------- ! ! Subroutine restores a string attribute from a given location. ! ! Arguments: ! ! loc_id - the location in which the dataset is stored; ! name - the attribute's name; ! string - the attribute's buffer; ! status - the subroutine call status; ! !=============================================================================== ! subroutine restore_attribute_string_h5(loc_id, name, string, status) use hdf5 use helpers , only : print_message use iso_c_binding, only : c_loc implicit none integer(hid_t) , intent(in) :: loc_id character(len=*) , intent(in) :: name character(len=*), target, intent(inout) :: string integer , intent(out) :: status integer(hid_t) :: attr_id, type_id, mem_id integer(hsize_t) :: length, attr_size logical :: flag type(c_ptr) :: str_ptr character(len=*), parameter :: loc = 'IO::restore_attribute_string_h5' !------------------------------------------------------------------------------- ! call h5aexists_by_name_f(loc_id, '.', trim(name), flag, status) if (status /= 0) then call print_message(loc, & "Could not check if attribute '" // trim(name) // "' exists!") return end if if (flag) then call h5aopen_by_name_f(loc_id, '.', trim(name), attr_id, status) if (status /= 0) then call print_message(loc, & "Could not open attribute '" // trim(name) // "'!") return end if call h5tcopy_f(H5T_NATIVE_CHARACTER, type_id, status) if (status /= 0) then call print_message(loc, & "Could not copy the datatype for attribute '" // trim(name) // "'!") go to 1000 end if call h5aget_storage_size_f(attr_id, attr_size, status) if (status /= 0) then call print_message(loc, & "Could not get the datatype size of attribute '" // & trim(name) // "'!") go to 1000 end if call h5tset_size_f(type_id, attr_size, status) if (status /= 0) then call print_message(loc, & "Could not set the datatype size for attribute '" // & trim(name) // "'!") go to 1000 end if call h5aget_type_f(attr_id, mem_id, status) if (status /= 0) then call print_message(loc, & "Could not get the datatype of attribute '" // trim(name) // "'!") go to 1000 end if call h5tequal_f(type_id, mem_id, flag, status) if (status /= 0) then call print_message(loc, & "The datatypes of the input string and attribute '" // & trim(name) // "' do not match!") go to 1000 end if length = len(string) if (length < attr_size) then call print_message(loc, & "The string is too small for storing attribute '" // & trim(name) // "'!") go to 1000 end if if (flag) then string = "" str_ptr = c_loc(string) call h5aread_f(attr_id, type_id, str_ptr, status) if (status /= 0) then call print_message(loc, & "Could not read attribute '" // trim(name) // "'!") end if end if 1000 continue call h5aclose_f(attr_id, status) if (status /= 0) & call print_message(loc, & "Could not close attribute '" // trim(name) // "'!") else call print_message(loc, "Attribute '" // trim(name) // "' not found!") end if !------------------------------------------------------------------------------- ! end subroutine restore_attribute_string_h5 ! !=============================================================================== ! ! subroutine RESTORE_ATTRIBUTE_NUMBER_H5: ! -------------------------------------- ! ! Subroutine restores a number (integer or real) attribute from a given ! location. Both scalar and vectors are supported. ! ! Arguments: ! ! loc_id - the location in which the dataset is stored; ! name - the attribute's name; ! type_id - the HDF5 type of attribute; ! length - the number of attribute's elements; ! buffer - the attribute's data; ! status - the subroutine call status; ! !=============================================================================== ! subroutine restore_attribute_number_h5(loc_id, name, type_id, & length, buffer, status) use hdf5 use helpers , only : print_message use iso_c_binding, only : c_loc implicit none integer(hid_t) , intent(in) :: loc_id, type_id character(len=*) , intent(in) :: name integer , intent(in) :: length type(*), target , dimension(..), intent(inout) :: buffer integer , intent(out) :: status integer(hid_t) :: attr_id, spc_id, mem_id integer(hsize_t), dimension(1) :: dims, mdims logical :: flag integer :: rank, cls_type type(c_ptr) :: buffer_ptr character(len=*), parameter :: loc = 'IO::restore_attribute_h5' !------------------------------------------------------------------------------- ! call h5aexists_by_name_f(loc_id, '.', trim(name), flag, status) if (status /= 0) then call print_message(loc, & "Could not check if attribute '" // trim(name) // "' exists!") return end if if (flag) then call h5aopen_by_name_f(loc_id, '.', trim(name), attr_id, status) if (status /= 0) then call print_message(loc, & "Could not open attribute '" // trim(name) // "'!") return end if call h5aget_type_f(attr_id, mem_id, status) if (status /= 0) then call print_message(loc, & "Could not get the datatype of attribute '" // trim(name) // "'!") go to 1000 end if call h5tequal_f(type_id, mem_id, flag, status) if (status /= 0) then call print_message(loc, & "The datatypes of the input buffer and attribute '" // & trim(name) // "' do not match!") go to 1000 end if if (flag) then call h5aget_space_f(attr_id, spc_id, status) if (status /= 0) then call print_message(loc, & "Could not get the dataspace of attribute '" // trim(name) // "'!") go to 1000 end if call h5sget_simple_extent_type_f(spc_id, cls_type, status) if (status /= 0) then call print_message(loc, & "Could not get the dataspace type for attribute '" // & trim(name) // "'!") go to 900 end if if (cls_type == H5S_SCALAR_F) then buffer_ptr = c_loc(buffer) call h5aread_f(attr_id, type_id, buffer_ptr, status) if (status /= 0) then call print_message(loc, & "Could not read attribute '" // trim(name) // "'!") end if else if (cls_type == H5S_SIMPLE_F) then call h5sget_simple_extent_dims_f(spc_id, dims, mdims, rank) if (rank /= 1) then call print_message(loc, "Only rank equal 1 is supported!") go to 800 end if mdims(1) = length if (dims(1) /= mdims(1)) then call print_message(loc, & "The dataspace dimensions the input buffer and attribute '" // & trim(name) // "' do not match!") go to 800 end if buffer_ptr = c_loc(buffer) call h5aread_f(attr_id, type_id, buffer_ptr, status) if (status /= 0) then call print_message(loc, & "Could not read attribute '" // trim(name) // "'!") end if 800 continue else call print_message(loc, & "The dataspace type of attribute '" // trim(name) // & "' is not supported!") end if 900 continue call h5sclose_f(spc_id, status) if (status /= 0) & call print_message(loc, & "Could not close the dataspace of attribute '" // trim(name) // "'!") end if 1000 continue call h5aclose_f(attr_id, status) if (status /= 0) & call print_message(loc, & "Could not close attribute '" // trim(name) // "'!") else call print_message(loc, "Attribute '" // trim(name) // "' not found!") end if !------------------------------------------------------------------------------- ! end subroutine restore_attribute_number_h5 ! !=============================================================================== ! ! subroutine STORE_ATTRIBUTE_STRING_H5: ! ------------------------------------ ! ! Subroutine stores a string attribute in a given location. ! ! Arguments: ! ! loc_id - the location in which the dataset is stored; ! name - the attribute's name; ! string - the attribute's text; ! status - the subroutine call status; ! !=============================================================================== ! subroutine store_attribute_string_h5(loc_id, name, string, status) use hdf5 use helpers , only : print_message use iso_c_binding, only : c_loc implicit none integer(hid_t) , intent(in) :: loc_id character(len=*) , intent(in) :: name character(len=*), target, intent(in) :: string integer , intent(out) :: status integer(hid_t) :: mem_id, spc_id, attr_id integer(hsize_t) :: length character(len=*), parameter :: loc = 'IO::store_attribute_string_h5()' !------------------------------------------------------------------------------- ! call h5tcopy_f(H5T_NATIVE_CHARACTER, mem_id, status) if (status /= 0) then call print_message(loc, & "Could not copy the datatype for attribute '" // trim(name) // "'!") return end if length = len(trim(string)) call h5tset_size_f(mem_id, length, status) if (status /= 0) then call print_message(loc, & "Could not set the datatype size for attribute '" // trim(name) // "'!") return end if call h5screate_f(H5S_SCALAR_F, spc_id, status) if (status /= 0) then call print_message(loc, & "Could not create the dataspace for attribute '" // trim(name) // "'!") return end if call h5acreate_f(loc_id, trim(name), mem_id, spc_id, attr_id, status) if (status /= 0) then call print_message(loc, & "Could not create attribute '" // trim(name) // "'!") go to 1000 end if call h5awrite_f(attr_id, mem_id, c_loc(string), status) if (status /= 0) & call print_message(loc, & "Could not write attribute " // trim(name) // "'!") call h5aclose_f(attr_id, status) if (status /= 0) & call print_message(loc, & "Could not close attribute '" // trim(name) // "'!") 1000 continue call h5sclose_f(spc_id, status) if (status /= 0) & call print_message(loc, & "Could not close the dataspace for attribute '" // trim(name) // "'!") !------------------------------------------------------------------------------- ! end subroutine store_attribute_string_h5 ! !=============================================================================== ! ! subroutine STORE_ATTRIBUTE_NUMBER_H5: ! ------------------------------------ ! ! Subroutine stores a number (integer or real) attribute in a given location. ! Both scalar and vectors are supported. ! ! Arguments: ! ! loc_id - the location in which the dataset is stored; ! name - the attribute's name; ! type_id - the HDF5 type of attribute; ! length - the number of attribute elements ! buffer - the attribute's data; ! status - the subroutine call status; ! !=============================================================================== ! subroutine store_attribute_number_h5(loc_id, name, type_id, & length, buffer, status) use hdf5 use helpers , only : print_message use iso_c_binding, only : c_loc implicit none integer(hid_t) , intent(in) :: loc_id, type_id character(len=*) , intent(in) :: name integer , intent(in) :: length type(*), target , dimension(..), intent(in) :: buffer integer , intent(out) :: status integer(hid_t) :: spc_id, attr_id integer(hsize_t), dimension(1) :: dims = 1 character(len=*), parameter :: loc = 'IO::store_attribute_number_h5()' !------------------------------------------------------------------------------- ! dims(1) = length if (length > 1) then call h5screate_simple_f(1, dims, spc_id, status) else call h5screate_f(H5S_SCALAR_F, spc_id, status) end if if (status /= 0) then call print_message(loc, & "Could not create the dataspace for attribute '" // trim(name) // "'!") return end if call h5acreate_f(loc_id, trim(name), type_id, spc_id, attr_id, status) if (status /= 0) then call print_message(loc, & "Could not create attribute '" // trim(name) // "'!") go to 1000 end if call h5awrite_f(attr_id, type_id, c_loc(buffer), status) if (status /= 0) & call print_message(loc, & "Could not write attribute " // trim(name) // "'!") call h5aclose_f(attr_id, status) if (status /= 0) & call print_message(loc, & "Could not close attribute '" // trim(name) // "'!") 1000 continue call h5sclose_f(spc_id, status) if (status /= 0) & call print_message(loc, & "Could not close the dataspace for attribute '" // trim(name) // "'!") !------------------------------------------------------------------------------- ! end subroutine store_attribute_number_h5 ! !=============================================================================== ! ! READ_ATTRIBUTE SUBROUTINES ! !=============================================================================== ! ! subroutine READ_SCALAR_ATTRIBUTE_INTEGER_H5: ! ------------------------------------------- ! ! Subroutine reads a value of the integer attribute provided by the group ! identifier to which it is linked and its name. ! ! Arguments: ! ! gid - the group identifier to which the attribute is linked; ! aname - the attribute name; ! avalue - the attribute value; ! !=============================================================================== ! subroutine read_scalar_attribute_integer_h5(gid, aname, avalue) ! import procedures and variables from other modules ! use hdf5 , only : H5T_NATIVE_INTEGER use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5aexists_by_name_f use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! attribute arguments ! integer(hid_t) , intent(in) :: gid character(len=*), intent(in) :: aname integer(kind=4) , intent(inout) :: avalue ! local variables ! logical :: exists = .false. integer(hid_t) :: aid integer(hsize_t), dimension(1) :: am = (/ 1 /) integer :: ierr ! local parameters ! character(len=*), parameter :: loc = 'IO::read_scalar_attribute_integer_h5()' ! !------------------------------------------------------------------------------- ! ! check if the attribute exists in the group provided by gid ! call h5aexists_by_name_f(gid, '.', aname, exists, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot check if attribute exists :" // trim(aname) return end if if (.not. exists) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Attribute does not exist :" // trim(aname) return end if ! open the attribute ! call h5aopen_by_name_f(gid, '.', aname, aid, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open attribute :" // trim(aname) return end if ! read attribute value ! call h5aread_f(aid, H5T_NATIVE_INTEGER, avalue, am(:), ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot read attribute :" // trim(aname) end if ! close the attribute ! call h5aclose_f(aid, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close attribute :" // trim(aname) return end if !------------------------------------------------------------------------------- ! end subroutine read_scalar_attribute_integer_h5 ! !=============================================================================== ! ! subroutine READ_SCALAR_ATTRIBUTE_DOUBLE_H5: ! ------------------------------------------ ! ! Subroutine reads a value of the double precision attribute provided by ! the group identifier to which it is linked and its name. ! ! Arguments: ! ! gid - the group identifier to which the attribute is linked; ! aname - the attribute name; ! avalue - the attribute value; ! !=============================================================================== ! subroutine read_scalar_attribute_double_h5(gid, aname, avalue) ! import procedures and variables from other modules ! use hdf5 , only : H5T_NATIVE_DOUBLE use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5aexists_by_name_f use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! attribute arguments ! integer(hid_t) , intent(in) :: gid character(len=*), intent(in) :: aname real(kind=8) , intent(inout) :: avalue ! local variables ! logical :: exists = .false. integer(hid_t) :: aid integer(hsize_t), dimension(1) :: am = (/ 1 /) integer :: ierr ! local parameters ! character(len=*), parameter :: loc = 'IO::read_scalar_attribute_double_h5()' ! !------------------------------------------------------------------------------- ! ! check if the attribute exists in the group provided by gid ! call h5aexists_by_name_f(gid, '.', aname, exists, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot check if attribute exists :" // trim(aname) return end if if (.not. exists) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Attribute does not exist :" // trim(aname) return end if ! open the attribute ! call h5aopen_by_name_f(gid, '.', aname, aid, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open attribute :" // trim(aname) return end if ! read attribute value ! call h5aread_f(aid, H5T_NATIVE_DOUBLE, avalue, am(:), ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot read attribute :" // trim(aname) end if ! close the attribute ! call h5aclose_f(aid, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close attribute :" // trim(aname) return end if !------------------------------------------------------------------------------- ! end subroutine read_scalar_attribute_double_h5 ! !=============================================================================== ! ! subroutine READ_ARRAY_ATTRIBUTE_LONG_H5: ! --------------------------------------- ! ! Subroutine reads a 2D array of the integer attribute provided by the group ! identifier to which it is linked and its name. ! ! Arguments: ! ! gid - the group identifier to which the attribute is linked; ! aname - the attribute name; ! avalue - the attribute value; ! !=============================================================================== ! subroutine read_array_attribute_long_h5(gid, aname, avalue) ! import procedures and variables from other modules ! use hdf5 , only : H5T_STD_I64LE use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5aexists_by_name_f, h5aget_space_f use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f use hdf5 , only : h5sclose_f, h5sget_simple_extent_dims_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! attribute arguments ! integer(hid_t) , intent(in) :: gid character(len=*) , intent(in) :: aname integer(kind=8) , dimension(:,:), intent(inout) :: avalue ! local variables ! logical :: exists = .false. integer(hid_t) :: aid, sid integer(hsize_t), dimension(2) :: am, bm integer :: ierr ! local parameters ! character(len=*), parameter :: loc = 'IO::read_array_attribute_long_h5()' ! !------------------------------------------------------------------------------- ! ! check if the attribute exists in the group provided by gid ! call h5aexists_by_name_f(gid, '.', aname, exists, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot check if attribute exists :" // trim(aname) return end if if (.not. exists) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Attribute does not exist :" // trim(aname) return end if ! open the attribute ! call h5aopen_by_name_f(gid, '.', aname, aid, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open attribute :" // trim(aname) return end if ! get the attribute space ! call h5aget_space_f(aid, sid, ierr) if (ierr == 0) then call h5sget_simple_extent_dims_f(sid, am, bm, ierr) if (ierr /= 2) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot get attribute dimensions :" // trim(aname) end if call h5sclose_f(sid, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close the attribute space :" // trim(aname) end if else write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot get the attribute space :" // trim(aname) return end if ! check if the output array is large enough ! if (am(1) * am(2) > size(avalue)) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Attribute too large for output argument :" // trim(aname) return end if ! read attribute value ! call h5aread_f(aid, H5T_STD_I64LE, avalue, am(:), ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot read attribute :" // trim(aname) end if ! close the attribute ! call h5aclose_f(aid, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close attribute :" // trim(aname) return end if !------------------------------------------------------------------------------- ! end subroutine read_array_attribute_long_h5 ! !=============================================================================== ! ! subroutine READ_ARRAY_ATTRIBUTE_COMPLEX_H5: ! ------------------------------------------ ! ! Subroutine reads a 2D array of the double precision complex attribute ! provided by the group identifier to which it is linked and its name. ! ! Arguments: ! ! gid - the group identifier to which the attribute is linked; ! aname - the attribute name; ! avalue - the attribute value; ! !=============================================================================== ! subroutine read_array_attribute_complex_h5(gid, aname, avalue) ! import procedures and variables from other modules ! use hdf5 , only : H5T_NATIVE_DOUBLE use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5aexists_by_name_f, h5aget_space_f use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f use hdf5 , only : h5sclose_f, h5sget_simple_extent_dims_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! attribute arguments ! integer(hid_t) , intent(in) :: gid character(len=*) , intent(in) :: aname complex(kind=8) , dimension(:,:), intent(inout) :: avalue ! local variables ! logical :: exists = .false. integer(hid_t) :: aid, sid integer(hsize_t), dimension(3) :: am, bm integer :: ierr ! local allocatable arrays ! real(kind=8), dimension(:,:,:), allocatable :: tvalue ! local parameters ! character(len=*), parameter :: loc = 'IO::read_array_attribute_complex_h5()' ! !------------------------------------------------------------------------------- ! ! check if the attribute exists in the group provided by gid ! call h5aexists_by_name_f(gid, '.', aname, exists, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot check if attribute exists :" // trim(aname) return end if if (.not. exists) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Attribute does not exist :" // trim(aname) return end if ! open the attribute ! call h5aopen_by_name_f(gid, '.', aname, aid, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open attribute :" // trim(aname) return end if ! get the attribute space ! call h5aget_space_f(aid, sid, ierr) if (ierr == 0) then call h5sget_simple_extent_dims_f(sid, am, bm, ierr) if (ierr /= 3) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot get attribute dimensions :" // trim(aname) end if call h5sclose_f(sid, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close the attribute space :" // trim(aname) end if else write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot get the attribute space :" // trim(aname) return end if ! allocate temporary array for attribute ! allocate(tvalue(am(1),am(2),am(3))) ! check if the output array is large enough ! if (am(1) * am(2) * am(3) > size(tvalue)) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Attribute too large for output argument :" // trim(aname) return end if ! read attribute value ! call h5aread_f(aid, H5T_NATIVE_DOUBLE, tvalue, am(:), ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot read attribute :" // trim(aname) end if ! copy array to complex one ! avalue(:,:) = cmplx(tvalue(1,:,:), tvalue(2,:,:), kind=8) ! deallocate temporary array ! deallocate(tvalue) ! close the attribute ! call h5aclose_f(aid, ierr) if (ierr /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close attribute :" // trim(aname) return end if !------------------------------------------------------------------------------- ! end subroutine read_array_attribute_complex_h5 ! !=============================================================================== ! ! subroutine READ_DATASET_H5: ! --------------------------- ! ! Subroutine reads a generic dataset. ! ! Arguments: ! ! loc_id - the location in which the dataset is stored; ! name - the dataset name; ! type_id - the dataset type; ! dims - the dataset dimensions; ! buffer - the dataset buffer; ! status - the subroutine call status; ! !=============================================================================== ! subroutine read_dataset_h5(loc_id, name, type_id, dims, buffer, status) use hdf5 use helpers , only : print_message use iso_c_binding, only : c_loc implicit none integer(hid_t) , intent(in) :: loc_id, type_id character(len=*) , intent(in) :: name integer(hsize_t), dimension(:), intent(in) :: dims type(*), target, dimension(..), intent(inout) :: buffer integer , intent(out) :: status type(c_ptr) :: buffer_ptr logical :: flag integer :: rank integer(hid_t) :: space_id, dset_id, mem_id integer(hsize_t), dimension(size(dims)) :: ddims, mdims character(len=*), parameter :: loc = 'IO::read_dataset_h5()' !------------------------------------------------------------------------------- ! status = 0 call h5dopen_f(loc_id, trim(name), dset_id, status) if (status /= 0) then call print_message(loc, "Could not open dataset '" // trim(name) // "'!") status = 1 return end if call h5dget_type_f(dset_id, mem_id, status) if (status /= 0) then call print_message(loc, & "Could not get the datatype for dataset '" // trim(name) // "'!") go to 1000 end if call h5tequal_f(type_id, mem_id, flag, status) if (status < 0) then call print_message(loc, & "Could not compare the buffer and dataset '" // trim(name) // & "' types!") go to 1000 end if if (.not. flag) then call print_message(loc, & "Datatypes of the buffer and dataset '" // trim(name) // & "' do not match!") go to 1000 end if call h5dget_space_f(dset_id, space_id, status) if (status /= 0) then call print_message(loc, & "Could not get the dataspace for dataset '" // trim(name) // "'!") go to 1000 end if call h5sget_simple_extent_dims_f(space_id, ddims, mdims, rank) if (rank /= size(dims)) then call print_message(loc, "Wrong rank of dataset '" // trim(name) // "'!") status = 1 go to 900 end if if (any(ddims /= dims)) then call print_message(loc, & "Wrong dimensions of dataset '" // trim(name) // "'!") status = 1 go to 900 end if buffer_ptr = c_loc(buffer) call h5dread_f(dset_id, type_id, buffer_ptr, status) if (status /= 0) & call print_message(loc, "Could not read dataset '" // trim(name) // "'!") 900 continue call h5sclose_f(space_id, status) if (status /= 0) & call print_message(loc, & "Could not close the dataspace for dataset '" // trim(name) // "'!") 1000 continue call h5dclose_f(dset_id, status) if (status /= 0) & call print_message(loc, "Could not close dataset '" // trim(name) // "'!") !------------------------------------------------------------------------------- ! end subroutine read_dataset_h5 ! !=============================================================================== ! ! subroutine STORE_DATASET_H5: ! --------------------------- ! ! Subroutine stores a generic dataset. ! ! Arguments: ! ! loc_id - the location in which the dataset is stored; ! name - the dataset name; ! type_id - the dataset type; ! dims - the dataset dimensions; ! buffer - the dataset buffer to store; ! status - the subroutine call status; ! !=============================================================================== ! subroutine store_dataset_h5(loc_id, name, type_id, dims, buffer, status) use hdf5 use helpers , only : print_message use iso_c_binding, only : c_loc implicit none integer(hid_t) , intent(in) :: loc_id, type_id character(len=*) , intent(in) :: name integer(hsize_t), dimension(:), intent(in) :: dims type(*), target, dimension(..), intent(in) :: buffer integer , intent(out) :: status integer :: rank integer(hid_t) :: space_id, dset_id character(len=*), parameter :: loc = 'IO::store_dataset_h5()' !------------------------------------------------------------------------------- ! rank = size(dims) call h5screate_simple_f(rank, dims, space_id, status) if (status /= 0) then call print_message(loc, & "Could not create the dataspace for dataset '" // trim(name) // "'!") return end if call h5pset_chunk_f(prp_id, rank, dims, status) if (status /= 0) then call print_message(loc, & "Could not set the chunk size for dataset '" // trim(name) // "'!") go to 1000 end if call h5dcreate_f(loc_id, name, type_id, space_id, dset_id, status) if (status /= 0) then call print_message(loc, & "Could not create dataset '" // trim(name) // "'!") go to 1000 end if call h5dwrite_f(dset_id, type_id, c_loc(buffer), status) if (status /= 0) then call print_message(loc, "Could not store dataset '" // trim(name) // "'!") end if call h5dclose_f(dset_id, status) if (status /= 0) & call print_message(loc, "Could not close dataset '" // trim(name) // "'!") 1000 continue call h5sclose_f(space_id, status) if (status /= 0) & call print_message(loc, & "Could not close the dataspace for dataset '" // trim(name) // "'!") !------------------------------------------------------------------------------- ! end subroutine store_dataset_h5 ! !=============================================================================== ! ! READ_ARRAY SUBROUTINES ! !=============================================================================== ! ! subroutine READ_1D_ARRAY_INTEGER_H5: ! ----------------------------------- ! ! Subroutine restores a one-dimensional integer array from a group specified ! by identifier. ! ! Arguments: ! ! gid - the HDF5 group identifier ! name - the string name describing the array ! dm - the array dimensions ! value - the array values ! !=============================================================================== ! subroutine read_1d_array_integer_h5(gid, name, dm, var) ! import procedures and variables from other modules ! use hdf5 , only : H5T_NATIVE_INTEGER use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5dopen_f, h5dread_f, h5dclose_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer(hid_t) , intent(in) :: gid character(len=*) , intent(in) :: name integer(hsize_t), dimension(1), intent(inout) :: dm integer(kind=4) , dimension(:), intent(inout) :: var ! local variables ! integer(hid_t) :: did integer :: iret ! local parameters ! character(len=*), parameter :: loc = 'IO::read_1d_array_integer_h5()' ! !------------------------------------------------------------------------------- ! ! open the dataset ! call h5dopen_f(gid, name, did, iret) ! check if the dataset has been opened successfuly ! if (iret < 0) then ! print error about the problem with opening the data space ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open dataset: " // trim(name) ! quit the subroutine ! return end if ! read dataset data ! call h5dread_f(did, H5T_NATIVE_INTEGER, var(:), dm(1:1), iret) ! check if the dataset has been read successfuly ! if (iret > 0) then ! print error about the problem with reading the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot read dataset: " // trim(name) end if ! close the dataset ! call h5dclose_f(did, iret) ! check if the dataset has been closed successfuly ! if (iret > 0) then ! print error about the problem with closing the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close dataset: " // trim(name) end if !------------------------------------------------------------------------------- ! end subroutine read_1d_array_integer_h5 ! !=============================================================================== ! ! subroutine READ_2D_ARRAY_INTEGER_H5: ! ----------------------------------- ! ! Subroutine restores a two-dimensional integer array from a group specified ! by identifier. ! ! Arguments: ! ! gid - the HDF5 group identifier ! name - the string name describing the array ! dm - the array dimensions ! value - the array values ! !=============================================================================== ! subroutine read_2d_array_integer_h5(gid, name, dm, var) ! import procedures and variables from other modules ! use hdf5 , only : H5T_NATIVE_INTEGER use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5dopen_f, h5dread_f, h5dclose_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer(hid_t) , intent(in) :: gid character(len=*) , intent(in) :: name integer(hsize_t), dimension(2) , intent(inout) :: dm integer(kind=4) , dimension(:,:), intent(inout) :: var ! local variables ! integer(hid_t) :: did integer :: iret ! local parameters ! character(len=*), parameter :: loc = 'IO::read_2d_array_integer_h5()' ! !------------------------------------------------------------------------------- ! ! open the dataset ! call h5dopen_f(gid, name, did, iret) ! check if the dataset has been opened successfuly ! if (iret < 0) then ! print error about the problem with opening the data space ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open dataset: " // trim(name) ! quit the subroutine ! return end if ! read dataset data ! call h5dread_f(did, H5T_NATIVE_INTEGER, var(:,:), dm(1:2), iret) ! check if the dataset has been read successfuly ! if (iret > 0) then ! print error about the problem with reading the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot read dataset: " // trim(name) end if ! close the dataset ! call h5dclose_f(did, iret) ! check if the dataset has been closed successfuly ! if (iret > 0) then ! print error about the problem with closing the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close dataset: " // trim(name) end if !------------------------------------------------------------------------------- ! end subroutine read_2d_array_integer_h5 #if NDIMS == 2 ! !=============================================================================== ! ! subroutine READ_3D_ARRAY_INTEGER_H5: ! ----------------------------------- ! ! Subroutine restores a three-dimensional integer array from a group specified ! by identifier. ! ! Arguments: ! ! gid - the HDF5 group identifier ! name - the string name describing the array ! dm - the array dimensions ! value - the array values ! !=============================================================================== ! subroutine read_3d_array_integer_h5(gid, name, dm, var) ! import procedures and variables from other modules ! use hdf5 , only : H5T_NATIVE_INTEGER use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5dopen_f, h5dread_f, h5dclose_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer(hid_t) , intent(in) :: gid character(len=*) , intent(in) :: name integer(hsize_t), dimension(3) , intent(inout) :: dm integer(kind=4) , dimension(:,:,:), intent(inout) :: var ! local variables ! integer(hid_t) :: did integer :: iret ! local parameters ! character(len=*), parameter :: loc = 'IO::read_3d_array_integer_h5()' ! !------------------------------------------------------------------------------- ! ! open the dataset ! call h5dopen_f(gid, name, did, iret) ! check if the dataset has been opened successfuly ! if (iret < 0) then ! print error about the problem with opening the data space ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open dataset: " // trim(name) ! quit the subroutine ! return end if ! read dataset data ! call h5dread_f(did, H5T_NATIVE_INTEGER, var(:,:,:), dm(1:3), iret) ! check if the dataset has been read successfuly ! if (iret > 0) then ! print error about the problem with reading the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot read dataset: " // trim(name) end if ! close the dataset ! call h5dclose_f(did, iret) ! check if the dataset has been closed successfuly ! if (iret > 0) then ! print error about the problem with closing the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close dataset: " // trim(name) end if !------------------------------------------------------------------------------- ! end subroutine read_3d_array_integer_h5 #endif /* NDIMS == 2 */ ! !=============================================================================== ! ! subroutine READ_4D_ARRAY_INTEGER_H5: ! ----------------------------------- ! ! Subroutine restores a four-dimensional integer array from a group specified ! by identifier. ! ! Arguments: ! ! gid - the HDF5 group identifier ! name - the string name describing the array ! dm - the array dimensions ! value - the array values ! !=============================================================================== ! subroutine read_4d_array_integer_h5(gid, name, dm, var) ! import procedures and variables from other modules ! use hdf5 , only : H5T_NATIVE_INTEGER use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5dopen_f, h5dread_f, h5dclose_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer(hid_t) , intent(in) :: gid character(len=*) , intent(in) :: name integer(hsize_t), dimension(4) , intent(inout) :: dm integer(kind=4) , dimension(:,:,:,:), intent(inout) :: var ! local variables ! integer(hid_t) :: did integer :: iret ! local parameters ! character(len=*), parameter :: loc = 'IO::read_4d_array_integer_h5()' ! !------------------------------------------------------------------------------- ! ! open the dataset ! call h5dopen_f(gid, name, did, iret) ! check if the dataset has been opened successfuly ! if (iret < 0) then ! print error about the problem with opening the data space ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open dataset: " // trim(name) ! quit the subroutine ! return end if ! read dataset data ! call h5dread_f(did, H5T_NATIVE_INTEGER, var(:,:,:,:), dm(1:4), iret) ! check if the dataset has been read successfuly ! if (iret > 0) then ! print error about the problem with reading the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot read dataset: " // trim(name) end if ! close the dataset ! call h5dclose_f(did, iret) ! check if the dataset has been closed successfuly ! if (iret > 0) then ! print error about the problem with closing the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close dataset: " // trim(name) end if !------------------------------------------------------------------------------- ! end subroutine read_4d_array_integer_h5 #if NDIMS == 3 ! !=============================================================================== ! ! subroutine READ_5D_ARRAY_INTEGER_H5: ! ----------------------------------- ! ! Subroutine restores a five-dimensional integer array from a group specified ! by identifier. ! ! Arguments: ! ! gid - the HDF5 group identifier ! name - the string name describing the array ! dm - the array dimensions ! value - the array values ! !=============================================================================== ! subroutine read_5d_array_integer_h5(gid, name, dm, var) ! import procedures and variables from other modules ! use hdf5 , only : H5T_NATIVE_INTEGER use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5dopen_f, h5dread_f, h5dclose_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer(hid_t) , intent(in) :: gid character(len=*) , intent(in) :: name integer(hsize_t), dimension(5) , intent(inout) :: dm integer(kind=4) , dimension(:,:,:,:,:), intent(inout) :: var ! local variables ! integer(hid_t) :: did integer :: iret ! local parameters ! character(len=*), parameter :: loc = 'IO::read_5d_array_integer_h5()' ! !------------------------------------------------------------------------------- ! ! open the dataset ! call h5dopen_f(gid, name, did, iret) ! check if the dataset has been opened successfuly ! if (iret < 0) then ! print error about the problem with opening the data space ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open dataset: " // trim(name) ! quit the subroutine ! return end if ! read dataset data ! call h5dread_f(did, H5T_NATIVE_INTEGER, var(:,:,:,:,:), dm(1:5), iret) ! check if the dataset has been read successfuly ! if (iret > 0) then ! print error about the problem with reading the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot read dataset: " // trim(name) end if ! close the dataset ! call h5dclose_f(did, iret) ! check if the dataset has been closed successfuly ! if (iret > 0) then ! print error about the problem with closing the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close dataset: " // trim(name) end if !------------------------------------------------------------------------------- ! end subroutine read_5d_array_integer_h5 #endif /* NDIMS == 3 */ ! !=============================================================================== ! ! subroutine READ_1D_ARRAY_DOUBLE_H5: ! ---------------------------------- ! ! Subroutine restores a one-dimensional double precision array from a group ! specified by identifier. ! ! Arguments: ! ! gid - the HDF5 group identifier ! name - the string name describing the array ! dm - the array dimensions ! value - the array values ! !=============================================================================== ! subroutine read_1d_array_double_h5(gid, name, dm, var) ! import procedures and variables from other modules ! use hdf5 , only : H5T_NATIVE_DOUBLE use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5dopen_f, h5dread_f, h5dclose_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer(hid_t) , intent(in) :: gid character(len=*) , intent(in) :: name integer(hsize_t), dimension(1), intent(inout) :: dm real(kind=8) , dimension(:), intent(inout) :: var ! local variables ! integer(hid_t) :: did integer :: iret ! local parameters ! character(len=*), parameter :: loc = 'IO::read_1d_array_double_h5()' ! !------------------------------------------------------------------------------- ! ! open the dataset ! call h5dopen_f(gid, name, did, iret) ! check if the dataset has been opened successfuly ! if (iret < 0) then ! print error about the problem with opening the data space ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open dataset: " // trim(name) ! quit the subroutine ! return end if ! read dataset data ! call h5dread_f(did, H5T_NATIVE_DOUBLE, var(:), dm(1:1), iret) ! check if the dataset has been read successfuly ! if (iret > 0) then ! print error about the problem with reading the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot read dataset: " // trim(name) end if ! close the dataset ! call h5dclose_f(did, iret) ! check if the dataset has been closed successfuly ! if (iret > 0) then ! print error about the problem with closing the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close dataset: " // trim(name) end if !------------------------------------------------------------------------------- ! end subroutine read_1d_array_double_h5 ! !=============================================================================== ! ! subroutine READ_4D_ARRAY_DOUBLE_H5: ! ----------------------------------- ! ! Subroutine restores a four-dimensional double precision array from a group ! specified by identifier. ! ! Arguments: ! ! gid - the HDF5 group identifier ! name - the string name describing the array ! dm - the array dimensions ! value - the array values ! !=============================================================================== ! subroutine read_4d_array_double_h5(gid, name, dm, var) ! import procedures and variables from other modules ! use hdf5 , only : H5T_NATIVE_DOUBLE use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5dopen_f, h5dread_f, h5dclose_f use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer(hid_t) , intent(in) :: gid character(len=*) , intent(in) :: name integer(hsize_t), dimension(4) , intent(inout) :: dm real(kind=8) , dimension(:,:,:,:), intent(inout) :: var ! local variables ! integer(hid_t) :: did integer :: iret ! local parameters ! character(len=*), parameter :: loc = 'IO::read_4d_array_double_h5()' ! !------------------------------------------------------------------------------- ! ! open the dataset ! call h5dopen_f(gid, name, did, iret) ! check if the dataset has been opened successfuly ! if (iret < 0) then ! print error about the problem with opening the data space ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot open dataset: " // trim(name) ! quit the subroutine ! return end if ! read dataset data ! call h5dread_f(did, H5T_NATIVE_DOUBLE, var(:,:,:,:), dm(1:4), iret) ! check if the dataset has been read successfuly ! if (iret > 0) then ! print error about the problem with reading the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot read dataset: " // trim(name) end if ! close the dataset ! call h5dclose_f(did, iret) ! check if the dataset has been closed successfuly ! if (iret > 0) then ! print error about the problem with closing the dataset ! write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot close dataset: " // trim(name) end if !------------------------------------------------------------------------------- ! end subroutine read_4d_array_double_h5 ! !=============================================================================== ! ! subroutine WRITE_SNAPSHOT_XDMF: ! ------------------------------ ! ! Subroutine writes an XDMF file per snapshot per MPI process. ! XDMF file is just a wrapper which helps to load data in a visualization ! tools like Paraview or Visit. ! ! Based on the subroutine by Pierre Kestener ! (see https://bitbucket.org/pkestene/amun-code). ! ! !=============================================================================== ! subroutine write_snapshot_xdmf() ! import procedures and variables from other modules ! use blocks , only : block_data, list_data use blocks , only : get_dblocks use equations , only : nv, pvars use mpitools , only : nproc use coordinates , only : ni => ncells, ng => nghosts use coordinates , only : adx, ady #if NDIMS == 3 use coordinates , only : adz #endif /* NDIMS == 3 */ use evolution , only : time ! local variables are not implicit by default ! implicit none ! local pointers ! type(block_data), pointer :: pdata ! local variables ! character(len=64) :: fname, hname character(len=128) :: stmp, ttmp, sdim, bdim, pdim integer(kind=4) :: l, p integer(kind=4) :: ip, jp #if NDIMS == 3 integer(kind=4) :: kp #endif /* NDIMS == 3 */ ! local arrays ! integer, dimension(12) :: slab ! local parameters ! integer, parameter :: xdmf = 101 ! !------------------------------------------------------------------------------- ! ! prepare the XDMF and HDF5 file names write(fname, "('p',i6.6,'_',i5.5,'.xdmf')") isnap, nproc write(hname, "('p',i6.6,'_',i5.5,'.h5' )") isnap, nproc ! open the XDMF file ! open (unit = xdmf, file = fname, status = 'replace') ! write the header ! write(xdmf, "(a)") '' write(xdmf, "(a)") '' write(xdmf, "(a)") ' ' write(stmp, "(1i16)") nproc write(xdmf, "(a)") ' ' write(stmp, "(1g15.8)") time write(xdmf, "(a)") ' ' write(xdmf, "(a)") ' ' write(xdmf, "(a)") '' ! close the XDMF file ! close(xdmf) !------------------------------------------------------------------------------- ! end subroutine write_snapshot_xdmf ! !=============================================================================== ! ! subroutine WRITE_SNAPSHOT_XDMF_MASTER: ! ------------------------------------- ! ! Subroutine writes one XDMF file per snapshot in root MPI process to gather ! all MPI subdomains. ! ! Based on the subroutine by Pierre Kestener ! (see https://bitbucket.org/pkestene/amun-code). ! ! !=============================================================================== ! subroutine write_snapshot_xdmf_master() ! import procedures and variables from other modules ! use mpitools , only : npmax ! local variables are not implicit by default ! implicit none ! local variables ! character(len=64) :: fname, pname integer(kind=4) :: np ! local parameters ! integer, parameter :: xdmf = 102 ! !------------------------------------------------------------------------------- ! ! prepare the XDMF and HDF5 file names write(fname, "('p',i6.6,'.xdmf')") isnap ! open the XDMF file ! open (unit = xdmf, file = fname, status = 'replace') ! write the header ! write(xdmf, "(a)") '' write(xdmf, "(a)") '' write(xdmf, "(a)") ' ' write(xdmf, "(a)") ' ' ! write references to MPI subdomain files ! do np = 0, npmax write(pname, "('p',i6.6,'_',i5.5,'.xdmf')") isnap, np write(xdmf, "(a)") ' ' end do ! close the XDMF structures ! write(xdmf, "(a)") ' ' write(xdmf, "(a)") ' ' write(xdmf, "(a)") '' ! close the XDMF file ! close(xdmf) !------------------------------------------------------------------------------- ! end subroutine write_snapshot_xdmf_master #endif /* HDF5 */ !=============================================================================== ! end module