HDF5_utilities.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------!
!------------------------------------------------------------------------------!
module hdf5_utilities
#include <PB3D_macros.h>
    use num_vars, only: max_str_ln, dp, plot_dir, data_dir, script_dir
    use messages
    use str_utilities, only: i2str, r2str, r2strt
    use hdf5_vars
    use hdf5
    
    implicit none
    private
    public probe_hdf5_group, set_1d_vars
#if ldebug
    public list_all_vars_in_group, &
        &debug_set_1d_vars
#endif
    
#if ldebug
    ! global variables
    logical :: debug_set_1d_vars = .false.                                      
#endif
    
contains
    integer function set_1d_vars(lim_tot,lim_loc,space_id,c_plist_id) &
        &result(ierr)
#if ldebug
        use num_vars, only: rank
#endif
        
        character(*), parameter :: rout_name = 'set_1D_vars'
        
        ! input / output
        integer, intent(in) :: lim_tot(:,:)
        integer, intent(in) :: lim_loc(:,:)
        integer(HID_T), intent(in), optional :: space_id
        integer(HID_T), intent(inout), optional :: c_plist_id
        
        ! local variables
        integer :: id, kd                                                       ! counters
        integer :: n_dims                                                       ! nr. of dimensions
        integer :: n_prod(2)                                                    ! n_prod(1) = n_1*n_2*..*n_(i-1), n_prod(2) = n_prod(1)*n_i
        integer :: div_dim                                                      ! first divided dimension
        integer :: chunk_size                                                   ! size of chunk
        integer(HSIZE_T) :: block(1)                                            ! block size in memory
        integer(HSIZE_T) :: stride(1)                                           ! stride in memory
        integer(HSIZE_T) :: count(1)                                            ! nr. of repetitions of block in memory
        integer(HSIZE_T) :: offset(1)                                           ! offset in memory
        integer(HSIZE_T) :: chunk_dims(1)                                       ! chunk dimensions
        real(dp) :: max_chunk_size                                              ! maximum 4 GB (from manual)
        real(dp) :: min_chunk_size                                              ! minimum 10 MB (from experience)
        character(len=max_str_ln) :: err_msg                                    ! error message
#if ldebug
        integer :: istat                                                        ! status
#endif
        
        ! initialize ierr
        ierr = 0
        
        ! set n_dims
        n_dims = size(lim_tot,1)
 
#if ldebug
        if (size(lim_loc,1).ne.n_dims) then
            ierr = 1
            err_msg = 'lim_tot and lim_loc are not compatible'
            chckerr(err_msg)
        end if
        if (size(lim_tot,2).ne.2 .or. size(lim_loc,2).ne.2) then
            ierr = 1
            err_msg = 'lim_tot and_lim_loc need to contain 2 columns'
            chckerr(err_msg)
        end if
        do id = 1,n_dims
            if (lim_tot(id,1).gt.lim_loc(id,1) .or. &
                &lim_tot(id,2).lt.lim_loc(id,2)) then
                write(*,*) rank, 'lim_tot = ', lim_tot
                write(*,*) rank, 'lim_loc = ', lim_loc
                ierr = 1
                err_msg = 'Total limits must comprise local ones'
                chckerr(err_msg)
            end if
        end do
        
        if (debug_set_1d_vars) then
            write(*,*,iostat=istat) rank, 'in set_1D_vars:'
            write(*,*,iostat=istat) rank, 'lim_tot_lo = ', lim_tot(:,1)
            write(*,*,iostat=istat) rank, 'lim_tot_hi = ', lim_tot(:,2)
            write(*,*,iostat=istat) rank, 'lim_loc_lo = ', lim_loc(:,1)
            write(*,*,iostat=istat) rank, 'lim_loc_hi = ', lim_loc(:,2)
        end if
#endif
        
        ! initialize divided dimension
        div_dim = n_dims+1
        
        ! set hyperslab to be everything if requested
        if (present(space_id)) then
            call h5sselect_all_f(space_id,ierr) 
            chckerr('Failed to select hyperslab')
        end if
        
        ! loop over dimensions to possibly restrict hyperslab and set div_dim
        do id = 1,n_dims
            ! only restrict hyperslab if local range differ from total one
            if (lim_loc(id,1).gt.lim_tot(id,1) .or. &
                &lim_loc(id,2).lt.lim_tot(id,2)) then
                ! set divided dimension if first time
                if (div_dim.eq.n_dims+1) div_dim = id
                
                ! calculations if hyperslab restriction requested
                if (present(space_id)) then
                    ! set auxiliary variables
                    n_prod(1) = product(lim_tot(1:id-1,2)-lim_tot(1:id-1,1)+1)
                    n_prod(2) = n_prod(1)*(lim_tot(id,2)-lim_tot(id,1)+1)
                    
#if ldebug
                    if (debug_set_1d_vars) then
                        write(*,*,iostat=istat) rank, 'dimension', id, 'of', &
                            &n_dims
                        write(*,*,iostat=istat) rank, 'n_prod = ', n_prod
                    end if
#endif
                    
                    block = (lim_loc(id,2)-lim_loc(id,1)+1) * n_prod(1)
                    stride = n_prod(2)
                    offset = (lim_loc(id,1)-lim_tot(id,1)) * n_prod(1)
                    count = 1
                    do kd = id+1,n_dims
                        count = count*(lim_tot(kd,2)-lim_tot(kd,1)+1)
                    end do
                    
#if ldebug
                    if (debug_set_1d_vars) then
                        write(*,*,iostat=istat) rank, 'block, offset, stride, &
                            &count = ', block, offset, stride, count
                    end if
#endif
                    
                    call h5sselect_hyperslab_f(space_id,h5s_select_and_f,&
                        &offset,count,ierr,stride=stride,block=block)
                    chckerr('Failed to select hyperslab')
                end if
#if ldebug
                if (debug_set_1d_vars) then
                    write(*,*,iostat=istat) rank, 'total range DIFFERS from &
                        &local'
                end if
            else
                if (debug_set_1d_vars) then
                    write(*,*,iostat=istat) rank, 'total range EQUAL to local'
                end if
#endif
            end if
        end do
        
        ! set chunk property list if requested
        max_chunk_size = 4.e9_dp / sizeof(1._dp)                                ! maximum 4 GB (from manual)
        min_chunk_size = min(10.e6_dp/sizeof(1._dp),&
            &product(lim_tot(:,2)-lim_tot(:,1)+1._dp))                          ! maximum 10 MB (from experience), or limit to total variable size
        if (present(c_plist_id)) then
            ! set up variables
            chunk_size = 1
            id = 1
            do while (chunk_size.le.min_chunk_size .and. id.le.size(lim_tot,1))
                chunk_size = chunk_size * (lim_tot(id,2)-lim_tot(id,1)+1)
                id = id + 1
            end do
#if ldebug
            if (debug_set_1d_vars) then
                write(*,*,iostat=istat) rank, 'suggested chunk size', chunk_size
            end if
#endif
            chunk_size = chunk_size / (1 + floor(chunk_size/max_chunk_size))    ! limit to max chunk size
#if ldebug
            if (debug_set_1d_vars) then
                write(*,*,iostat=istat) rank, 'def. chunk size', chunk_size
            end if
#endif
            
            ! set creation property list if provided
            if (present(c_plist_id)) then
                chunk_dims = chunk_size
                call h5pcreate_f(h5p_dataset_create_f,c_plist_id,ierr) 
                chckerr('Failed to create property list')
                call h5pset_chunk_f(c_plist_id,1,chunk_dims,ierr)               ! one dimension for 1D chunk
                chckerr('Failed to create chunk')
            end if
        end if
    end function set_1d_vars
    
    integer function probe_hdf5_group(HDF5_name,group_name,group_exists) &
        &result(ierr)
        use mpi_vars, only: hdf5_lock
        use mpi_utilities, only: lock_req_acc, lock_return_acc
        
        character(*), parameter :: rout_name = 'probe_HDF5_group'
        
        ! input / output
        character(len=*), intent(in) :: hdf5_name
        character(len=*), intent(in) :: group_name
        logical, intent(inout) :: group_exists
        
        ! local variables
        integer :: istat                                                        ! status
        integer(HID_T) :: hdf5_i                                                ! file identifier 
        integer(HID_T) :: group_id                                              ! group identifier
        
        ! initialize ierr
        ierr = 0
        
        ! initialize FORTRAN predefined datatypes
        call h5open_f(ierr) 
        chckerr('Failed to initialize HDF5')
        
        ! wait for file access in a non-blocking way
        ierr = lock_req_acc(hdf5_lock,blocking=.false.)
        chckerr('')
        
        ! open the file
        call h5fopen_f(hdf5_name,h5f_acc_rdonly_f,hdf5_i,ierr)
        chckerr('Failed to open file. Does it exist?')
        
        ! disable error messages
        call h5eset_auto_f(0,ierr)
        chckerr('Failed to disable error printing')
        
        ! try to open group
        call h5gopen_f(hdf5_i,group_name,group_id,istat)
        group_exists = istat.eq.0
        
        ! reenable error messages
        call h5eset_auto_f(1,ierr)
        chckerr('Failed to enable error printing')
        
        ! close group
        if (group_exists) then
            call h5gclose_f(group_id,ierr)
            chckerr('Failed to close head group')
        end if
        
        ! close the HDF5 file
        call h5fclose_f(hdf5_i,ierr)
        chckerr('failed to close HDF5 file')
        
        ! return lock
        ierr = lock_return_acc(hdf5_lock)
        chckerr('')
        
        ! close FORTRAN interfaces and HDF5 library.
        call h5close_f(ierr)
        chckerr('Failed to close FORTRAN HDF5 interface')
    end function probe_hdf5_group
 
#if ldebug
 
    integer function list_all_vars_in_group(group_id) result(ierr)
        character(*), parameter :: rout_name = 'list_all_vars_in_group'
        
        ! input / output
        integer(HID_T), intent(in) :: group_id
        
        ! local variables
        integer :: storage_type                                                 ! type of storage used in HDF5 file
        integer :: nr_lnks                                                      ! number of links in group
        integer :: max_corder                                                   ! current maximum creation order value for group
        integer(HSIZE_T) :: id                                                  ! counter
        integer :: id_loc                                                       ! local id
        integer(SIZE_T) :: name_len                                             ! length of name of group
        character(len=max_str_ln) :: var_name                                   ! name of variable
        
        ! get number of objects in group
        call h5gget_info_f(group_id,storage_type,nr_lnks,max_corder,&
            &ierr)
        chckerr('Failed to get group info')
        
        ! user output
        call writo('The group has '//trim(i2str(nr_lnks))//' elements')
        
        ! loop over all elements
        do id = 1,nr_lnks
            ! get variable name
            call h5lget_name_by_idx_f(group_id,'.',h5_index_name_f,&
                &h5_iter_native_f,id-1,var_name,ierr,size=name_len)
            chckerr('Failed to get name')
            
            ! user output
            id_loc = int(id,1)
            call writo('Element '//trim(i2str(id_loc))//'/'//&
                &trim(i2str(nr_lnks))//': '//trim(var_name))
        end do
    end function list_all_vars_in_group
#endif
end module hdf5_utilities