diff --git a/unstructured/M3Dmodules.f90 b/unstructured/M3Dmodules.f90 index c87e5770c..6ab5eff17 100644 --- a/unstructured/M3Dmodules.f90 +++ b/unstructured/M3Dmodules.f90 @@ -424,6 +424,8 @@ module basic integer :: iread_neo ! 1 = read velocity profiles from NEO output integer :: ineo_subtract_diamag ! 1 = subtract v* from input v profile + integer :: islice_serial=0 ! (=0; slice output using parallel hdf5), (=1; using serial hdf5) + ! adaptation options integer :: iadapt ! 1,2 = adapts mesh after initialization real :: adapt_factor diff --git a/unstructured/hdf5_output.f90 b/unstructured/hdf5_output.f90 index 535e848fc..eef90ed4a 100644 --- a/unstructured/hdf5_output.f90 +++ b/unstructured/hdf5_output.f90 @@ -360,6 +360,9 @@ end subroutine write_vec_attr ! ================= subroutine output_field(parent_id, name, values, ndofs, nelms, error) use hdf5 + use h5lt + use mpi + use basic implicit none @@ -374,100 +377,140 @@ subroutine output_field(parent_id, name, values, ndofs, nelms, error) integer(HSIZE_T), dimension(rank) :: local_dims, global_dims integer(HSSIZE_T), dimension(rank) :: off + integer, parameter :: sp = selected_real_kind(6,37) + real, dimension(ndofs,global_elms) :: recv_Gvalues + real(sp), dimension(:,:), allocatable :: Gvalues_sp ! single precision data + integer(HSIZE_T) dtype_id + integer recv_counts_Lsize(maxrank), recv_displs_Lsize(maxrank) + integer i + #ifdef USETAU integer :: dummy ! this is necessary to prevent TAU from dummy = 0 ! breaking formatting requirements #endif - - local_dims(1) = ndofs - local_dims(2) = nelms - global_dims(1) = ndofs - global_dims(2) = global_elms - off(1) = 0 - off(2) = offset - - ! Create global dataset - call h5screate_simple_f(rank, global_dims, filespace, error) - if(error.ne.0) then - write(*,*) error,rank," after h5screate_simple_f" - call safestop(101) - endif - if(idouble_out.eq.1) then - call h5dcreate_f(parent_id, name, H5T_NATIVE_DOUBLE, filespace, & - dset_id, error) - else - call h5dcreate_f(parent_id, name, H5T_NATIVE_REAL, filespace, & - dset_id, error) - end if - if(error.ne.0) then - write(*,*) error,rank," after h5dcreate_f" - call safestop(101) - endif - call h5sclose_f(filespace, error) - if(error.ne.0) then - write(*,*) error,rank," after hsclose_f" - call safestop(101) - endif - - ! Select local hyperslab within dataset - call h5screate_simple_f(rank, local_dims, memspace, error) - if(error.ne.0) then - write(*,*) error,rank," after h5screate_simple_f" - call safestop(102) - endif - call h5dget_space_f(dset_id, filespace, error) - if(error.ne.0) then - write(*,*) error,rank," after h5dget_space_f" - call safestop(102) - endif - call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, off, local_dims, & - error) - if(error.ne.0) then - write(*,*) error,rank," after h5sselect_hyperslab_f" - call safestop(102) - endif - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error) - if(error.ne.0) then - write(*,*) error,rank," after h5pcreate_f" - call safestop(102) - endif - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, error) - if(error.ne.0) then - write(*,*) error,rank," after h5pset_dxpl_mpio_f" - call safestop(102) - endif - - ! Write the dataset - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, values, global_dims, error, & - file_space_id=filespace, mem_space_id=memspace, xfer_prp=plist_id) - if(error.ne.0) then - write(*,*) error,rank," h5dwrite_f" - call safestop(103) + + if(islice_serial==0) then + local_dims(1) = ndofs + local_dims(2) = nelms + global_dims(1) = ndofs + global_dims(2) = global_elms + off(1) = 0 + off(2) = offset + + ! Create global dataset + call h5screate_simple_f(rank, global_dims, filespace, error) + if(error.ne.0) then + write(*,*) error,rank," after h5screate_simple_f" + call safestop(101) + endif + if(idouble_out.eq.1) then + call h5dcreate_f(parent_id, name, H5T_NATIVE_DOUBLE, filespace, & + dset_id, error) + else + call h5dcreate_f(parent_id, name, H5T_NATIVE_REAL, filespace, & + dset_id, error) + end if + if(error.ne.0) then + write(*,*) error,rank," after h5dcreate_f" + call safestop(101) + endif + call h5sclose_f(filespace, error) + if(error.ne.0) then + write(*,*) error,rank," after hsclose_f" + call safestop(101) + endif + + ! Select local hyperslab within dataset + call h5screate_simple_f(rank, local_dims, memspace, error) + if(error.ne.0) then + write(*,*) error,rank," after h5screate_simple_f" + call safestop(102) + endif + call h5dget_space_f(dset_id, filespace, error) + if(error.ne.0) then + write(*,*) error,rank," after h5dget_space_f" + call safestop(102) + endif + call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, off, local_dims, & + error) + if(error.ne.0) then + write(*,*) error,rank," after h5sselect_hyperslab_f" + call safestop(102) + endif + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error) + if(error.ne.0) then + write(*,*) error,rank," after h5pcreate_f" + call safestop(102) + endif + call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, error) + if(error.ne.0) then + write(*,*) error,rank," after h5pset_dxpl_mpio_f" + call safestop(102) + endif + + ! Write the dataset + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, values, global_dims, error, & + file_space_id=filespace, mem_space_id=memspace, xfer_prp=plist_id) + if(error.ne.0) then + write(*,*) error,rank," h5dwrite_f" + call safestop(103) + endif + + ! Close HDF5 handles + call h5sclose_f(filespace, error) + if(error.ne.0) then + write(*,*) error,rank," h5sclose_f" + call safestop(104) + endif + call h5sclose_f(memspace, error) + if(error.ne.0) then + write(*,*) error,rank," h5sclose_f" + call safestop(105) + endif + + call h5dclose_f(dset_id, error) + if(error.ne.0) then + write(*,*) error,rank," h5dclose_f" + call safestop(104) + endif + call h5pclose_f(plist_id, error) + if(error.ne.0) then + write(*,*) error,rank," h5pclose_f" + call safestop(105) + endif + else ! islice_serial==1 + global_dims(1) = ndofs + global_dims(2) = global_elms + + ! Get counts and displacements of Lsize in MPI_COMM_WORLD + call MPI_ALLGATHER(ndofs*nelms,1,MPI_INTEGER,recv_counts_Lsize,1,MPI_INTEGER,MPI_COMM_WORLD,error) + recv_displs_Lsize(1)=0 + do i=1,maxrank-1 + recv_displs_Lsize(i+1)=recv_displs_Lsize(i)+recv_counts_Lsize(i) + enddo + + call MPI_Gatherv(values, ndofs*nelms, MPI_DOUBLE_PRECISION, recv_Gvalues, recv_counts_Lsize, recv_displs_Lsize, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, error) + + if(myrank==0) then + if(idouble_out.eq.1) then + call h5ltmake_dataset_f(parent_id,name,rank,global_dims,H5T_NATIVE_DOUBLE, recv_Gvalues, error) + else + allocate(Gvalues_sp(ndofs,global_elms)) + Gvalues_sp = real(recv_Gvalues, kind=sp) + dtype_id = h5kind_to_type(sp,H5_REAL_KIND) + call h5screate_simple_f(rank, global_dims, filespace, error) + call h5dcreate_f(parent_id, name, dtype_id, filespace, dset_id, error) + call h5dwrite_f(dset_id, dtype_id, Gvalues_sp, global_dims, error) + call h5dclose_f(dset_id,error) + call h5sclose_f(filespace,error) + deallocate(Gvalues_sp) + end if + if(error.ne.0) then + write(*,*) error,rank," Error. h5ltmake_dataset_f" + call safestop(106) + endif + endif endif - - ! Close HDF5 handles - call h5sclose_f(filespace, error) - if(error.ne.0) then - write(*,*) error,rank," h5sclose_f" - call safestop(104) - endif - call h5sclose_f(memspace, error) - if(error.ne.0) then - write(*,*) error,rank," h5sclose_f" - call safestop(105) - endif - - call h5dclose_f(dset_id, error) - if(error.ne.0) then - write(*,*) error,rank," h5dclose_f" - call safestop(104) - endif - call h5pclose_f(plist_id, error) - if(error.ne.0) then - write(*,*) error,rank," h5pclose_f" - call safestop(105) - endif - end subroutine output_field diff --git a/unstructured/input.f90 b/unstructured/input.f90 index 2fa7317a1..6234260e4 100644 --- a/unstructured/input.f90 +++ b/unstructured/input.f90 @@ -1037,6 +1037,8 @@ subroutine set_defaults "1: Use double-precision floating points in output", output_grp) call add_var_int("irestart_slice", irestart_slice, -1, & "Field output slice from which to restart", output_grp) + call add_var_int("islice_serial", islice_serial, 0, & + "Field ouput slice using parallel (0) or serial (1) hdf5", output_grp) call add_var_int("iveldif", iveldif, 0, & "ne.0: veldif plot contains only partial results ", output_grp) diff --git a/unstructured/output.f90 b/unstructured/output.f90 index 5aeee5ace..8c79813c2 100644 --- a/unstructured/output.f90 +++ b/unstructured/output.f90 @@ -620,6 +620,10 @@ subroutine hdf5_write_time_slice(equilibrium, error) integer(HID_T) :: time_file_id, time_root_id, plist_id integer :: info logical :: link_exists + logical :: flag_call=.true. + + ! In the case of serial run, non-root processors do not call some functions + if(islice_serial==1 .and. myrank/=0) flag_call=.false. call hdf5_get_local_elms(nelms, error) @@ -643,43 +647,46 @@ subroutine hdf5_write_time_slice(equilibrium, error) ! Create new file for timeslice ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - ! Set up the file access property list with parallel I/O - call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, error) - info = MPI_INFO_NULL - call h5pset_fapl_mpio_f(plist_id, MPI_COMM_WORLD, info, error) - - ! Open the new file - call h5fcreate_f(time_file_name, H5F_ACC_TRUNC_F, time_file_id, error, & - access_prp = plist_id) + if(islice_serial==0) then + ! Set up the file access property list with parallel I/O + call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, error) + info = MPI_INFO_NULL + call h5pset_fapl_mpio_f(plist_id, MPI_COMM_WORLD, info, error) + ! Open the new file + call h5fcreate_f(time_file_name, H5F_ACC_TRUNC_F, time_file_id, error, & + access_prp = plist_id) + else + if(myrank==0) call h5fcreate_f(time_file_name, H5F_ACC_TRUNC_F, time_file_id, error) + endif if(error.lt.0) then print *, "Error: could not open ", time_file_name, & " for HDF5 output. error = ", error return endif - + ! open the root group - call h5gopen_f(time_file_id, "/", time_root_id, error) + if(flag_call) call h5gopen_f(time_file_id, "/", time_root_id, error) if(myrank.eq.0 .and. iprint.ge.1) & print *, ' Writing time slice file ', time_file_name ! Write attributes if(myrank.eq.0 .and. iprint.ge.1) print *, ' Writing attr ' - call write_real_attr(time_root_id, "time", time, error) + if(flag_call) call write_real_attr(time_root_id, "time", time, error) #ifdef USE3D - call write_int_attr(time_root_id, "nspace", 3, error) + if(flag_call) call write_int_attr(time_root_id, "nspace", 3, error) #else - call write_int_attr(time_root_id, "nspace", 2, error) + if(flag_call) call write_int_attr(time_root_id, "nspace", 2, error) #endif ! Time step associated with this time slice if(equilibrium.eq.1) then - call write_int_attr(time_root_id, "ntimestep", 0, error) + if(flag_call) call write_int_attr(time_root_id, "ntimestep", 0, error) else - call write_int_attr(time_root_id, "ntimestep", ntime, error) + if(flag_call) call write_int_attr(time_root_id, "ntimestep", ntime, error) end if ! Write version number - call write_int_attr(time_root_id, "version", version, error) + if(flag_call) call write_int_attr(time_root_id, "version", version, error) ! Output the mesh data if(myrank.eq.0 .and. iprint.ge.1) print *, ' Writing mesh ' @@ -693,13 +700,13 @@ subroutine hdf5_write_time_slice(equilibrium, error) ! output wall regions #ifndef USE3D if(myrank.eq.0 .and. iprint.ge.1) print *, ' Writing wall regions' - call output_regions(time_root_id, error) + if(flag_call) call output_regions(time_root_id, error) #endif ! Close the file - call h5gclose_f(time_root_id, error) - call h5fclose_f(time_file_id, error) - call h5pclose_f(plist_id, error) + if(flag_call) call h5gclose_f(time_root_id, error) + if(flag_call) call h5fclose_f(time_file_id, error) + if(islice_serial==0) call h5pclose_f(plist_id, error) ! Add timeslice link in main file @@ -731,6 +738,7 @@ subroutine output_mesh(time_group_id, nelms, error) use mesh_mod use basic use boundary_conditions + use mpi implicit none @@ -756,24 +764,30 @@ subroutine output_mesh(time_group_id, nelms, error) integer :: idim(3) real :: bound integer :: izone + logical :: flag_call=.true. - ! Create the group - call h5gcreate_f(time_group_id, "mesh", mesh_group_id, error) + ! In the case of serial run, non-root processors do not call some functions + if(islice_serial==1 .and. myrank/=0) flag_call=.false. - ! Write attributes - call write_int_attr(mesh_group_id, "nelms", global_elms, error) - call get_bounding_box_size(alx, alz) - call write_real_attr(mesh_group_id, "width", alx, error) - call write_real_attr(mesh_group_id, "height", alz, error) + if(flag_call) then + ! Create the group + call h5gcreate_f(time_group_id, "mesh", mesh_group_id, error) + + ! Write attributes + call write_int_attr(mesh_group_id, "nelms", global_elms, error) + call get_bounding_box_size(alx, alz) + call write_real_attr(mesh_group_id, "width", alx, error) + call write_real_attr(mesh_group_id, "height", alz, error) #ifdef USE3D - call write_int_attr(mesh_group_id, "3D", 1, error) + call write_int_attr(mesh_group_id, "3D", 1, error) #else - call write_int_attr(mesh_group_id, "3D", 0, error) + call write_int_attr(mesh_group_id, "3D", 0, error) #endif - call write_int_attr(mesh_group_id, "nplanes", nplanes, error) - call write_int_attr(mesh_group_id, "nperiods", nperiods, error) - call write_int_attr(mesh_group_id, "ifull_torus", ifull_torus, error) - call write_real_attr(mesh_group_id, "period", toroidal_period, error) + call write_int_attr(mesh_group_id, "nplanes", nplanes, error) + call write_int_attr(mesh_group_id, "nperiods", nperiods, error) + call write_int_attr(mesh_group_id, "ifull_torus", ifull_torus, error) + call write_real_attr(mesh_group_id, "period", toroidal_period, error) + endif ! Output the mesh data do i=1, nelms @@ -809,18 +823,20 @@ subroutine output_mesh(time_group_id, nelms, error) call output_field(mesh_group_id, "elements", elm_data, vals_per_elm, & nelms, error) + if(flag_call) then #ifdef USE3D - allocate(phi(nplanes)) - do i=1, nplanes - call m3dc1_plane_getphi(i-1, phi(i)) - end do - call write_vec_attr(mesh_group_id, "phi", phi, nplanes, error) - deallocate(phi) + allocate(phi(nplanes)) + do i=1, nplanes + call m3dc1_plane_getphi(i-1, phi(i)) + end do + call write_vec_attr(mesh_group_id, "phi", phi, nplanes, error) + deallocate(phi) #endif + ! Close the group + call h5gclose_f(mesh_group_id, error) + endif - ! Close the group - call h5gclose_f(mesh_group_id, error) end subroutine output_mesh subroutine output_regions(group_id, error) @@ -929,6 +945,10 @@ subroutine output_fields(time_group_id, equilibrium, error) integer :: i, nelms, ilin vectype, allocatable :: dum(:,:), dum2(:,:) character(len=64) :: field_name + logical :: flag_call = .true. + + ! In the case of serial run, non-root processors do not call some functions + if(islice_serial==1 .and. myrank/=0) flag_call=.false. ilin = 1 - equilibrium @@ -940,7 +960,7 @@ subroutine output_fields(time_group_id, equilibrium, error) allocate(dum(coeffs_per_element,nelms)) ! Create the fields group - call h5gcreate_f(time_group_id, "fields", group_id, error) + if(flag_call) call h5gcreate_f(time_group_id, "fields", group_id, error) ! Output the fields ! ~~~~~~~~~~~~~~~~~ @@ -1257,7 +1277,7 @@ subroutine output_fields(time_group_id, equilibrium, error) !!$ end if ! Close the mesh group - call h5gclose_f(group_id, error) + if(flag_call) call h5gclose_f(group_id, error) deallocate(dum)