! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! Maintainers : support@fluidnumerics.com ! Official Repository : https://github.com/FluidNumerics/self/ ! ! Copyright © 2024 Fluid Numerics LLC ! ! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in ! the documentation and/or other materials provided with the distribution. ! ! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from ! this software without specific prior written permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT ! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT ! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY ! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF ! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! module SELF_MultiSetIO_3D !! HDF5 helpers for writing and reading per-set element data in a multi-set model file, !! and for checking a dataset's extent before reading it. !! !! A set's data on rank r covers set-global elements elemOffset+1 .. elemOffset+nLocal of !! a dataset whose last dimension is the set's global element count. A rank may own no !! elements of a set (nLocal = 0); it still takes part in every collective call and !! selects nothing in both the file and memory spaces. SELF_HDF5's parallel array !! routines always select the local shape, which is why these exist. !! !! Collectives. With parallel = .true., dataset creation and the transfers are collective !! over the file's communicator: every rank must call each routine for every dataset, in !! the same order. With parallel = .false. the file is serial and the whole array is !! written or read. use SELF_Constants use HDF5 implicit none contains subroutine WriteSetArray_HDF5(fileId,arrayName,hfArray,elemOffset,nGlobalElem,parallel) !! Create arrayName with shape (Np,Np,Np,nGlobalElem) and write this rank's block !! hfArray(1:Np,1:Np,1:Np,1:nLocal) at element offset elemOffset. !! !! Input !! - fileId : open HDF5 file !! - arrayName : dataset path; its group must exist !! - hfArray(Np,Np,Np,nLocal) : this rank's elements; nLocal may be 0 !! - elemOffset : set-global elements on lower ranks !! - nGlobalElem : set-global element count, >= 1 !! - parallel : .true. for a file opened with the MPI-IO driver implicit none integer(HID_T),intent(in) :: fileId character(*),intent(in) :: arrayName real(prec),intent(in) :: hfArray(:,:,:,:) integer,intent(in) :: elemOffset integer,intent(in) :: nGlobalElem logical,intent(in) :: parallel ! Local integer(HID_T) :: dsetId,filespace,memspace,plistId,dtype integer(HSIZE_T) :: dims(1:4),globalDims(1:4),offset(1:4) integer :: error dtype = h5kind_to_type(prec,H5_REAL_KIND) dims = shape(hfArray) globalDims = (/dims(1),dims(2),dims(3),int(nGlobalElem,HSIZE_T)/) offset = (/0_hsize_t,0_hsize_t,0_hsize_t,int(elemOffset,HSIZE_T)/) call h5screate_simple_f(4,globalDims,filespace,error) call h5screate_simple_f(4,dims,memspace,error) call h5dcreate_f(fileId,trim(arrayName),dtype,filespace,dsetId,error) if(error /= 0) then print*,__FILE__," : failed to create ",trim(arrayName) error stop "WriteSetArray_HDF5 : dataset creation failed" endif if(dims(4) > 0) then call h5sselect_hyperslab_f(filespace,H5S_SELECT_SET_F,offset,dims,error) else call h5sselect_none_f(filespace,error) call h5sselect_none_f(memspace,error) endif call h5pcreate_f(H5P_DATASET_XFER_F,plistId,error) if(parallel) call h5pset_dxpl_mpio_f(plistId,H5FD_MPIO_COLLECTIVE_F,error) call h5dwrite_f(dsetId,dtype,hfArray,dims,error, & mem_space_id=memspace,file_space_id=filespace,xfer_prp=plistId) if(error /= 0) then print*,__FILE__," : failed to write ",trim(arrayName) error stop "WriteSetArray_HDF5 : write failed" endif call h5pclose_f(plistId,error) call h5sclose_f(memspace,error) call h5sclose_f(filespace,error) call h5dclose_f(dsetId,error) endsubroutine WriteSetArray_HDF5 subroutine ReadSetArray_HDF5(fileId,arrayName,hfArray,elemOffset,nGlobalElem,parallel) !! Read this rank's block of arrayName, elements elemOffset+1 .. elemOffset+nLocal, !! into hfArray(1:Np,1:Np,1:Np,1:nLocal). nLocal may be 0. !! !! The dataset must have shape (Np,Np,Np,nGlobalElem), matching hfArray's nodal !! extents: a dataset written at another polynomial degree holds values at other !! nodes, so this stops rather than reading a sub-block of it. The check uses only !! the dataset's metadata, so every rank, including one with nLocal = 0, reaches the !! same decision. implicit none integer(HID_T),intent(in) :: fileId character(*),intent(in) :: arrayName real(prec),intent(inout) :: hfArray(:,:,:,:) integer,intent(in) :: elemOffset integer,intent(in) :: nGlobalElem logical,intent(in) :: parallel ! Local integer(HID_T) :: dsetId,filespace,memspace,plistId,dtype integer(HSIZE_T) :: dims(1:4),offset(1:4),fileDims(1:4),maxDims(1:4) integer :: error,ndims dtype = h5kind_to_type(prec,H5_REAL_KIND) dims = shape(hfArray) offset = (/0_hsize_t,0_hsize_t,0_hsize_t,int(elemOffset,HSIZE_T)/) call h5dopen_f(fileId,trim(arrayName),dsetId,error) if(error /= 0) then print*,__FILE__," : failed to open ",trim(arrayName) error stop "ReadSetArray_HDF5 : dataset open failed" endif call h5dget_space_f(dsetId,filespace,error) call h5sget_simple_extent_ndims_f(filespace,ndims,error) if(ndims /= 4) then print*,__FILE__," : ",trim(arrayName)," has rank ",ndims,", expected 4" error stop "ReadSetArray_HDF5 : dataset shape mismatch" endif call h5sget_simple_extent_dims_f(filespace,fileDims,maxDims,error) if(any(fileDims(1:3) /= dims(1:3)) .or. fileDims(4) /= int(nGlobalElem,HSIZE_T)) then print*,__FILE__," : ",trim(arrayName)," has shape ",fileDims, & ", expected ",dims(1:3),nGlobalElem error stop "ReadSetArray_HDF5 : dataset shape mismatch" endif call h5screate_simple_f(4,dims,memspace,error) if(dims(4) > 0) then call h5sselect_hyperslab_f(filespace,H5S_SELECT_SET_F,offset,dims,error) else call h5sselect_none_f(filespace,error) call h5sselect_none_f(memspace,error) endif call h5pcreate_f(H5P_DATASET_XFER_F,plistId,error) if(parallel) call h5pset_dxpl_mpio_f(plistId,H5FD_MPIO_COLLECTIVE_F,error) call h5dread_f(dsetId,dtype,hfArray,dims,error, & mem_space_id=memspace,file_space_id=filespace,xfer_prp=plistId) if(error /= 0) then print*,__FILE__," : failed to read ",trim(arrayName) error stop "ReadSetArray_HDF5 : read failed" endif call h5pclose_f(plistId,error) call h5sclose_f(memspace,error) call h5sclose_f(filespace,error) call h5dclose_f(dsetId,error) endsubroutine ReadSetArray_HDF5 function DatasetLength1D_HDF5(fileId,arrayName) result(n) !! Number of entries in the one-dimensional dataset arrayName, or -1 if it does not !! exist or has a different rank. implicit none integer(HID_T),intent(in) :: fileId character(*),intent(in) :: arrayName integer :: n ! Local integer(HID_T) :: dsetId,filespace integer(HSIZE_T) :: dims(1:1),maxDims(1:1) integer :: error,ndims logical :: exists n = -1 call h5lexists_f(fileId,trim(arrayName),exists,error) if(error /= 0 .or. .not. exists) return call h5dopen_f(fileId,trim(arrayName),dsetId,error) if(error /= 0) return call h5dget_space_f(dsetId,filespace,error) call h5sget_simple_extent_ndims_f(filespace,ndims,error) if(ndims == 1) then call h5sget_simple_extent_dims_f(filespace,dims,maxDims,error) if(error >= 0) n = int(dims(1)) endif call h5sclose_f(filespace,error) call h5dclose_f(dsetId,error) endfunction DatasetLength1D_HDF5 endmodule SELF_MultiSetIO_3D