Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- PROGRAM phdf5write
- USE MPI
- USE H5LT
- USE timing
- IMPLICIT NONE
- INTEGER, PARAMETER :: minglobalsize = 100, maxglobalsize = 300, stepsize = 100 ! Size of the whole 3D domain (per side)
- ! Will do tests for all sizes from minglobalsize to maxglobalsize
- ! in steps of stepsize
- INTEGER, PARAMETER :: nblocks = 4 ! Number of blocks per dimension. with nblocks = 4 we end up
- ! with 64 blocks (4^3), so we will need to run the code with 64 processors
- ! !! Make sure that all the possible step are divisible by nblocks !!
- INTEGER, PARAMETER :: layersize = 24, nvar = 12 ! Variables controlling the size of the PML layer to be written.
- ! nvar is always the last dimension in the array datapml
- ! layersize is the first dimension for mode = 1, the second for mode = 2
- ! and the third for mode = 3
- !
- ! For example, when layersize=24 and nvar=24 and 3D size is 100, then
- ! the dimensions for PMLX will be 24x100x100x24
- !===========================================================================================
- DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: datapml
- INTEGER :: error, proc_number, my_rank
- INTEGER :: size, mysize, temp,x,y,z,timerx,timery,timerz
- CALL timer_create("WRITINGPMLX",timerx)
- CALL timer_create("WRITINGPMLY",timery)
- CALL timer_create("WRITINGPMLZ",timerz)
- CALL MPI_INIT ( error )
- CALL MPI_COMM_SIZE ( MPI_COMM_WORLD, proc_number, error )
- CALL MPI_COMM_RANK ( MPI_COMM_WORLD, my_rank, error )
- DO size=minglobalsize,maxglobalsize,stepsize
- IF (my_rank .EQ. 0) PRINT*, ""
- IF (my_rank .EQ. 0) PRINT*, "===================", size, "====================="
- IF (my_rank .EQ. 0) PRINT*, ""
- CALL timer_reset(timerx)
- CALL timer_reset(timery)
- CALL timer_reset(timerz)
- mysize = size / nblocks
- ALLOCATE(datapml(mysize,mysize,layersize,nvar))
- datapml = my_rank
- CALL timer_start(timerx)
- CALL write_data_pml(1)
- CALL timer_stop(timerx)
- CALL timer_start(timery)
- CALL write_data_pml(2)
- CALL timer_stop(timery)
- CALL timer_start(timerz)
- CALL write_data_pml(3)
- CALL timer_stop(timerz)
- DEALLOCATE(datapml)
- IF (my_rank .EQ. 0) CALL timer_report(0)
- END DO
- CALL MPI_Finalize (error)
- CONTAINS
- SUBROUTINE write_data_pml(mode)
- INTEGER :: mode
- CHARACTER(LEN=11) :: filename
- CHARACTER(LEN=4) :: modetxt
- INTEGER(HID_T) :: plist_id ! Property list identifier
- INTEGER(HID_T) :: file_id ! File identifier
- INTEGER(HID_T) :: ddata ! Dataset identifier
- INTEGER(HID_T) :: filespace ! Dataspace identifier in file
- INTEGER(HID_T) :: memspace ! Dataspace identifier in memory
- INTEGER(HSIZE_T), DIMENSION(4) :: dims,count
- INTEGER(HSIZE_T), DIMENSION(4) :: offset
- INTEGER :: error ! Error flag
- INTEGER :: info
- SELECT CASE (mode)
- CASE (1)
- modetxt = "PMLX"
- CASE (2)
- modetxt = "PMLY"
- CASE(3)
- modetxt = "PMLZ"
- END SELECT
- IF (my_rank .EQ. 0) PRINT*, "===================", modetxt, "====================="
- CALL h5open_f (error)
- dims = (/ size, size, layersize,nvar /)
- WRITE(filename,'(a4,i3,a3)'), modetxt, size, '.h5'
- ! Setup file access property list with parallel I/O access.
- !==========================================================
- CALL h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, error)
- !!$ call MPI_Info_create(info, error)
- !!$ call MPI_Info_set(info,"IBM_largeblock_io","true", error)
- CALL h5pset_fapl_mpio_f(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL, error)
- ! Create the file collectively.
- ! =============================
- CALL h5fcreate_f(filename, H5F_ACC_TRUNC_F, file_id, error, access_prp = plist_id)
- ! Close propertly list
- ! ====================
- CALL h5pclose_f(plist_id, error)
- CALL h5screate_simple_f(4, dims, filespace, error)
- CALL h5dcreate_f(file_id, "datapml", H5T_NATIVE_DOUBLE, filespace, ddata, error)
- CALL h5sclose_f(filespace, error)
- ! Each process defines dataset in memory and writes it to the hyperslab in the file.
- ! ==================================================================================
- count = (/ 0,0,0,0 /)
- SELECT CASE(mode)
- CASE(1)
- IF (my_rank < nblocks*nblocks) count = (/ mysize,mysize,layersize,nvar /) ! X PML
- CASE(2)
- IF (MOD(my_rank,nblocks*nblocks) < nblocks) count = (/ mysize,mysize,layersize,nvar /) ! Y PML
- CASE(3)
- IF (MOD(my_rank,nblocks) .EQ. 0) count = (/ mysize,mysize,layersize,nvar /) ! Z PML
- END SELECT
- offset = (/ 0,0,0,0 /)
- DO temp=0,nblocks*nblocks - 1
- x = temp / nblocks
- y = MOD(temp,nblocks)
- SELECT CASE(mode)
- CASE(1)
- IF (my_rank .EQ. temp) THEN ! Like a X PML layer
- offset(1) = x * count(1)
- offset(2) = y * count(2)
- PRINT*, "rank, offs", my_rank, offset(1), offset(2)
- END IF
- CASE(2)
- IF (my_rank .EQ. temp/nblocks * (nblocks*nblocks) + MOD(temp,nblocks)) THEN ! Like a Y PML layer
- offset(1) = x * count(1)
- offset(2) = y * count(2)
- PRINT*, "rank, offs", my_rank, offset(1), offset(2)
- END IF
- CASE(3)
- IF (my_rank .EQ. temp * nblocks) THEN ! Like a Z PML layer
- offset(1) = x * count(1)
- offset(2) = y * count(2)
- PRINT*, "rank, offs", my_rank, offset(1), offset(2)
- END IF
- END SELECT
- END DO
- CALL h5screate_simple_f(4, count, memspace, error)
- !
- ! Select hyperslab in the file.
- !
- CALL h5dget_space_f(ddata, filespace, error)
- CALL h5sselect_hyperslab_f (filespace, H5S_SELECT_SET_F, offset, count, error)
- !
- ! Create property list for collective dataset write
- !
- CALL h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error)
- CALL h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, error)
- ! Do the actual writing
- ! =====================
- CALL h5dwrite_f(ddata,H5T_NATIVE_DOUBLE,datapml,dims,error,file_space_id=filespace,&
- mem_space_id = memspace, xfer_prp = plist_id)
- ! Close dataspaces.
- ! ================
- CALL h5sclose_f(filespace, error)
- CALL h5sclose_f(memspace, error)
- CALL h5dclose_f(ddata, error)
- CALL h5pclose_f(plist_id, error)
- CALL h5fclose_f(file_id, error)
- END SUBROUTINE write_data_pml
- END PROGRAM phdf5write
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
