diff --git a/installation/mods_Abaqus/abaqus_v6.env b/installation/mods_Abaqus/abaqus_v6.env index 0811d0f65..83cc2ed33 100644 --- a/installation/mods_Abaqus/abaqus_v6.env +++ b/installation/mods_Abaqus/abaqus_v6.env @@ -16,7 +16,7 @@ if False: # use hdf5 compiler wrapper in $PATH fortCmd = os.popen('h5fc -shlib -show').read().replace('\n','') # complicated way needed to pass in DAMASKVERSION string link_sl += fortCmd.split()[1:] - fortCmd +=" -DDAMASKHDF5" + fortCmd +=" -DDAMASK_HDF5" else: # Use the version in $PATH fortCmd = "ifort" diff --git a/installation/mods_Abaqus/abaqus_v6_debug.env b/installation/mods_Abaqus/abaqus_v6_debug.env index 943f40bfa..943d0d10e 100644 --- a/installation/mods_Abaqus/abaqus_v6_debug.env +++ b/installation/mods_Abaqus/abaqus_v6_debug.env @@ -16,7 +16,7 @@ if False: # use hdf5 compiler wrapper in $PATH fortCmd = os.popen('h5fc -shlib -show').read().replace('\n','') # complicated way needed to pass in DAMASKVERSION string link_sl += fortCmd.split()[1:] - fortCmd +=" -DDAMASKHDF5" + fortCmd +=" -DDAMASK_HDF5" else: # Use the version in $PATH fortCmd = "ifort" diff --git a/installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64 b/installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64 index 661d3aaca..538434ad0 100644 --- a/installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64 +++ b/installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64 @@ -102,7 +102,7 @@ fi if test "$DAMASK_HDF5" = "ON";then H5FC="$(h5fc -shlib -show)" HDF5_LIB=${H5FC//ifort/} - FCOMP="$H5FC -DDAMASKHDF5" + FCOMP="$H5FC -DDAMASK_HDF5" echo $FCOMP else FCOMP=ifort diff --git a/installation/mods_MarcMentat/2018/Marc_tools/include_linux64 b/installation/mods_MarcMentat/2018/Marc_tools/include_linux64 index 270184af2..d3151ac6c 100644 --- a/installation/mods_MarcMentat/2018/Marc_tools/include_linux64 +++ b/installation/mods_MarcMentat/2018/Marc_tools/include_linux64 @@ -63,7 +63,6 @@ else INTEGER_PATH=/$MARC_INTEGER_SIZE fi -FCOMP=ifort INTELPATH="/opt/intel/compilers_and_libraries_2017/linux" # find the root directory of the compiler installation: @@ -99,6 +98,16 @@ else FCOMPROOT= fi +# DAMASK uses the HDF5 compiler wrapper around the Intel compiler +if test "$DAMASK_HDF5" = "ON";then + H5FC="$(h5fc -shlib -show)" + HDF5_LIB=${H5FC//ifort/} + FCOMP="$H5FC -DDAMASK_HDF5" + echo $FCOMP +else + FCOMP=ifort +fi + # AEM if test "$MARCDLLOUTDIR" = ""; then DLLOUTDIR="$MARC_LIB" @@ -535,6 +544,7 @@ else DAMASKVERSION="'N/A'" fi + # DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3 DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \ -fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018 -DDAMASKVERSION=$DAMASKVERSION \ @@ -737,7 +747,7 @@ SECLIBS="-L$MARC_LIB -llapi" SOLVERLIBS="${BCSSOLVERLIBS} ${VKISOLVERLIBS} ${CASISOLVERLIBS} ${MF2SOLVERLIBS} \ $MKLLIB -L$MARC_MKL -liomp5 \ - $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/kdtree2.a " + $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/kdtree2.a $HDF5_LIB " SOLVERLIBS_DLL=${SOLVERLIBS} if test "$AEM_DLL" -eq 1 diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index d34a79bf7..d2eaa7979 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -72,6 +72,12 @@ subroutine CPFEM_initAll(el,ip) mesh_init use material, only: & material_init +#ifdef DAMASK_HDF5 + use HDF5_utilities, only: & + HDF5_utilities_init + use results, only: & + results_init +#endif use lattice, only: & lattice_init use constitutive, only: & @@ -100,6 +106,10 @@ subroutine CPFEM_initAll(el,ip) call FE_init call mesh_init(ip, el) call lattice_init +#ifdef DAMASK_HDF5 + call HDF5_utilities_init + call results_init +#endif call material_init call constitutive_init call crystallite_init diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 13d7f06c4..45e423fe3 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -72,9 +72,9 @@ subroutine CPFEM_initAll() call FE_init call mesh_init call lattice_init - call material_init call HDF5_utilities_init call results_init + call material_init call constitutive_init call crystallite_init call homogenization_init @@ -300,6 +300,8 @@ subroutine CPFEM_results(inc,time) use HDF5_utilities use constitutive, only: & constitutive_results + use crystallite, only: & + crystallite_results implicit none integer(pInt), intent(in) :: inc @@ -307,7 +309,8 @@ subroutine CPFEM_results(inc,time) call results_openJobFile call results_addIncrement(inc,time) - call constitutive_results() + call constitutive_results + call crystallite_results call results_removeLink('current') ! ToDo: put this into closeJobFile call results_closeJobFile diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index a81aaee0e..a2593e1cb 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -5,11 +5,11 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !-------------------------------------------------------------------------------------------------- module HDF5_utilities - use prec - use IO - use HDF5 + use prec + use IO + use HDF5 #ifdef PETSc - use PETSC + use PETSC #endif implicit none @@ -19,76 +19,81 @@ module HDF5_utilities !> @brief reads pInt or pReal data of defined shape from file ! ToDo: order of arguments wrong !> @details for parallel IO, all dimension except for the last need to match !-------------------------------------------------------------------------------------------------- - interface HDF5_read - module procedure HDF5_read_pReal1 - module procedure HDF5_read_pReal2 - module procedure HDF5_read_pReal3 - module procedure HDF5_read_pReal4 - module procedure HDF5_read_pReal5 - module procedure HDF5_read_pReal6 - module procedure HDF5_read_pReal7 - - module procedure HDF5_read_pInt1 - module procedure HDF5_read_pInt2 - module procedure HDF5_read_pInt3 - module procedure HDF5_read_pInt4 - module procedure HDF5_read_pInt5 - module procedure HDF5_read_pInt6 - module procedure HDF5_read_pInt7 - - end interface HDF5_read + interface HDF5_read + module procedure HDF5_read_real1 + module procedure HDF5_read_real2 + module procedure HDF5_read_real3 + module procedure HDF5_read_real4 + module procedure HDF5_read_real5 + module procedure HDF5_read_real6 + module procedure HDF5_read_real7 + + module procedure HDF5_read_int1 + module procedure HDF5_read_int2 + module procedure HDF5_read_int3 + module procedure HDF5_read_int4 + module procedure HDF5_read_int5 + module procedure HDF5_read_int6 + module procedure HDF5_read_int7 + + end interface HDF5_read !-------------------------------------------------------------------------------------------------- !> @brief writes pInt or pReal data of defined shape to file ! ToDo: order of arguments wrong !> @details for parallel IO, all dimension except for the last need to match !-------------------------------------------------------------------------------------------------- - interface HDF5_write - module procedure HDF5_write_pReal1 - module procedure HDF5_write_pReal2 - module procedure HDF5_write_pReal3 - module procedure HDF5_write_pReal4 - module procedure HDF5_write_pReal5 - module procedure HDF5_write_pReal6 - module procedure HDF5_write_pReal7 - - module procedure HDF5_write_pInt1 - module procedure HDF5_write_pInt2 - module procedure HDF5_write_pInt3 - module procedure HDF5_write_pInt4 - module procedure HDF5_write_pInt5 - module procedure HDF5_write_pInt6 - module procedure HDF5_write_pInt7 - - end interface HDF5_write + interface HDF5_write + module procedure HDF5_write_real1 + module procedure HDF5_write_real2 + module procedure HDF5_write_real3 + module procedure HDF5_write_real4 + module procedure HDF5_write_real5 + module procedure HDF5_write_real6 + module procedure HDF5_write_real7 + + module procedure HDF5_write_int1 + module procedure HDF5_write_int2 + module procedure HDF5_write_int3 + module procedure HDF5_write_int4 + module procedure HDF5_write_int5 + module procedure HDF5_write_int6 + module procedure HDF5_write_int7 + + module procedure HDF5_write_rotation + + end interface HDF5_write !-------------------------------------------------------------------------------------------------- !> @brief attached attributes of type char,pInt or pReal to a file/dataset/group !-------------------------------------------------------------------------------------------------- - interface HDF5_addAttribute - module procedure HDF5_addAttribute_str - module procedure HDF5_addAttribute_pInt - module procedure HDF5_addAttribute_pReal - end interface HDF5_addAttribute + interface HDF5_addAttribute + module procedure HDF5_addAttribute_str + module procedure HDF5_addAttribute_int + module procedure HDF5_addAttribute_real + end interface HDF5_addAttribute !-------------------------------------------------------------------------------------------------- - public :: & - HDF5_utilities_init, & - HDF5_openFile, & - HDF5_closeFile, & - HDF5_addAttribute, & - HDF5_closeGroup ,& - HDF5_openGroup, & - HDF5_addGroup, & - HDF5_read, & - HDF5_write, & - HDF5_setLink, & - HDF5_objectExists + public :: & + HDF5_utilities_init, & + HDF5_openFile, & + HDF5_closeFile, & + HDF5_addAttribute, & + HDF5_closeGroup ,& + HDF5_openGroup, & + HDF5_addGroup, & + HDF5_read, & + HDF5_write, & + HDF5_setLink, & + HDF5_objectExists contains + +!-------------------------------------------------------------------------------------------------- +!> @brief open libary and do sanity checks +!-------------------------------------------------------------------------------------------------- subroutine HDF5_utilities_init - implicit none integer :: hdferr integer(SIZE_T) :: typeSize @@ -117,46 +122,45 @@ end subroutine HDF5_utilities_init !-------------------------------------------------------------------------------------------------- integer(HID_T) function HDF5_openFile(fileName,mode,parallel) ! ToDo: simply "open" is enough - implicit none - character(len=*), intent(in) :: fileName - character, intent(in), optional :: mode - logical, intent(in), optional :: parallel - - character :: m - integer(HID_T) :: plist_id - integer :: hdferr - - if (present(mode)) then - m = mode - else - m = 'r' - endif - - call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5pcreate_f') + character(len=*), intent(in) :: fileName + character, intent(in), optional :: mode + logical, intent(in), optional :: parallel + + character :: m + integer(HID_T) :: plist_id + integer :: hdferr + + if (present(mode)) then + m = mode + else + m = 'r' + endif + + call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5pcreate_f') #ifdef PETSc - if (present(parallel)) then; if (parallel) then - call h5pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5pset_fapl_mpio_f') - endif; endif + if (present(parallel)) then; if (parallel) then + call h5pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5pset_fapl_mpio_f') + endif; endif #endif - if (m == 'w') then - call h5fcreate_f(fileName,H5F_ACC_TRUNC_F,HDF5_openFile,hdferr,access_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fcreate_f (w)') - elseif(m == 'a') then - call h5fopen_f(fileName,H5F_ACC_RDWR_F,HDF5_openFile,hdferr,access_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fopen_f (a)') - elseif(m == 'r') then - call h5fopen_f(fileName,H5F_ACC_RDONLY_F,HDF5_openFile,hdferr,access_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fopen_f (r)') - else - call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fopen_f unknown access mode: '//trim(m)) - endif - - call h5pclose_f(plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5pclose_f') + if (m == 'w') then + call h5fcreate_f(fileName,H5F_ACC_TRUNC_F,HDF5_openFile,hdferr,access_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fcreate_f (w)') + elseif(m == 'a') then + call h5fopen_f(fileName,H5F_ACC_RDWR_F,HDF5_openFile,hdferr,access_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fopen_f (a)') + elseif(m == 'r') then + call h5fopen_f(fileName,H5F_ACC_RDONLY_F,HDF5_openFile,hdferr,access_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fopen_f (r)') + else + call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fopen_f unknown access mode: '//trim(m)) + endif + + call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5pclose_f') end function HDF5_openFile @@ -166,13 +170,12 @@ end function HDF5_openFile !-------------------------------------------------------------------------------------------------- subroutine HDF5_closeFile(fileHandle) - implicit none - integer(HID_T), intent(in) :: fileHandle - - integer :: hdferr - - call h5fclose_f(fileHandle,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_closeFile: h5fclose_f') + integer(HID_T), intent(in) :: fileHandle + + integer :: hdferr + + call h5fclose_f(fileHandle,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_closeFile: h5fclose_f') end subroutine HDF5_closeFile @@ -182,29 +185,30 @@ end subroutine HDF5_closeFile !-------------------------------------------------------------------------------------------------- integer(HID_T) function HDF5_addGroup(fileHandle,groupName) - implicit none - integer(HID_T), intent(in) :: fileHandle - character(len=*), intent(in) :: groupName + integer(HID_T), intent(in) :: fileHandle + character(len=*), intent(in) :: groupName + + integer :: hdferr + integer(HID_T) :: aplist_id - integer :: hdferr - integer(HID_T) :: aplist_id +!------------------------------------------------------------------------------------------------- +! creating a property list for data access properties + call h5pcreate_f(H5P_GROUP_ACCESS_F, aplist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup: h5pcreate_f ('//trim(groupName)//')') - !------------------------------------------------------------------------------------------------- - ! creating a property list for data access properties - call h5pcreate_f(H5P_GROUP_ACCESS_F, aplist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup: h5pcreate_f ('//trim(groupName)//')') - - !------------------------------------------------------------------------------------------------- - ! setting I/O mode to collective +!------------------------------------------------------------------------------------------------- +! setting I/O mode to collective #ifdef PETSc - call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup: h5pset_all_coll_metadata_ops_f ('//trim(groupName)//')') + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup: h5pset_all_coll_metadata_ops_f ('//trim(groupName)//')') #endif - !------------------------------------------------------------------------------------------------- - ! Create group - call h5gcreate_f(fileHandle, trim(groupName), HDF5_addGroup, hdferr, OBJECT_NAMELEN_DEFAULT_F,gapl_id = aplist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup: h5gcreate_f ('//trim(groupName)//')') +!------------------------------------------------------------------------------------------------- +! Create group + call h5gcreate_f(fileHandle, trim(groupName), HDF5_addGroup, hdferr, OBJECT_NAMELEN_DEFAULT_F,gapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup: h5gcreate_f ('//trim(groupName)//')') + + call h5pclose_f(aplist_id,hdferr) end function HDF5_addGroup @@ -214,32 +218,33 @@ end function HDF5_addGroup !-------------------------------------------------------------------------------------------------- integer(HID_T) function HDF5_openGroup(fileHandle,groupName) - implicit none - integer(HID_T), intent(in) :: fileHandle - character(len=*), intent(in) :: groupName - - - integer :: hdferr - integer(HID_T) :: aplist_id - logical :: is_collective + integer(HID_T), intent(in) :: fileHandle + character(len=*), intent(in) :: groupName + + + integer :: hdferr + integer(HID_T) :: aplist_id + logical :: is_collective !------------------------------------------------------------------------------------------------- ! creating a property list for data access properties - call h5pcreate_f(H5P_GROUP_ACCESS_F, aplist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup: h5pcreate_f ('//trim(groupName)//')') + call h5pcreate_f(H5P_GROUP_ACCESS_F, aplist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup: h5pcreate_f ('//trim(groupName)//')') !------------------------------------------------------------------------------------------------- ! setting I/O mode to collective #ifdef PETSc - call h5pget_all_coll_metadata_ops_f(aplist_id, is_collective, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup: h5pset_all_coll_metadata_ops_f ('//trim(groupName)//')') + call h5pget_all_coll_metadata_ops_f(aplist_id, is_collective, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup: h5pset_all_coll_metadata_ops_f ('//trim(groupName)//')') #endif !------------------------------------------------------------------------------------------------- ! opening the group - call h5gopen_f(fileHandle, trim(groupName), HDF5_openGroup, hdferr, gapl_id = aplist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup: h5gopen_f ('//trim(groupName)//')') + call h5gopen_f(fileHandle, trim(groupName), HDF5_openGroup, hdferr, gapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup: h5gopen_f ('//trim(groupName)//')') + + call h5pclose_f(aplist_id,hdferr) end function HDF5_openGroup @@ -249,12 +254,11 @@ end function HDF5_openGroup !-------------------------------------------------------------------------------------------------- subroutine HDF5_closeGroup(group_id) - implicit none - integer(HID_T), intent(in) :: group_id - integer :: hdferr - - call h5gclose_f(group_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_closeGroup: h5gclose_f (el is ID)', el = int(group_id,pInt)) + integer(HID_T), intent(in) :: group_id + integer :: hdferr + + call h5gclose_f(group_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_closeGroup: h5gclose_f (el is ID)', el = int(group_id,pInt)) end subroutine HDF5_closeGroup @@ -264,25 +268,24 @@ end subroutine HDF5_closeGroup !-------------------------------------------------------------------------------------------------- logical function HDF5_objectExists(loc_id,path) - implicit none - integer(HID_T), intent(in) :: loc_id - character(len=*), intent(in), optional :: path - integer :: hdferr - character(len=256) :: p + integer(HID_T), intent(in) :: loc_id + character(len=*), intent(in), optional :: path + integer :: hdferr + character(len=256) :: p + + if (present(path)) then + p = trim(path) + else + p = '.' + endif - if (present(path)) then - p = trim(path) - else - p = '.' - endif - - call h5lexists_f(loc_id, p, HDF5_objectExists, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_objectExists: h5oexists_by_name_f') - - if(HDF5_objectExists) then - call h5oexists_by_name_f(loc_id, p, HDF5_objectExists, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_objectExists: h5oexists_by_name_f') - endif + call h5lexists_f(loc_id, p, HDF5_objectExists, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_objectExists: h5oexists_by_name_f') + + if(HDF5_objectExists) then + call h5oexists_by_name_f(loc_id, p, HDF5_objectExists, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_objectExists: h5oexists_by_name_f') + endif end function HDF5_objectExists @@ -292,43 +295,42 @@ end function HDF5_objectExists !-------------------------------------------------------------------------------------------------- subroutine HDF5_addAttribute_str(loc_id,attrLabel,attrValue,path) - implicit none - integer(HID_T), intent(in) :: loc_id - character(len=*), intent(in) :: attrLabel, attrValue - character(len=*), intent(in), optional :: path - integer :: hdferr - integer(HID_T) :: attr_id, space_id, type_id - logical :: attrExists - character(len=256) :: p + integer(HID_T), intent(in) :: loc_id + character(len=*), intent(in) :: attrLabel, attrValue + character(len=*), intent(in), optional :: path + integer :: hdferr + integer(HID_T) :: attr_id, space_id, type_id + logical :: attrExists + character(len=256) :: p + + if (present(path)) then + p = trim(path) + else + p = '.' + endif - if (present(path)) then - p = trim(path) - else - p = '.' - endif - - call h5screate_f(H5S_SCALAR_F,space_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5screate_f') - call h5tcopy_f(H5T_NATIVE_CHARACTER, type_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5tcopy_f') - call h5tset_size_f(type_id, int(len(trim(attrValue)),HSIZE_T), hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5tset_size_f') - call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5aexists_by_name_f') - if (attrExists) then - call h5adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5adelete_by_name_f') - endif - call h5acreate_by_name_f(loc_id,trim(p),trim(attrLabel),type_id,space_id,attr_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5acreate_f') - call h5awrite_f(attr_id, type_id, trim(attrValue), int([1],HSIZE_T), hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5awrite_f') - call h5aclose_f(attr_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5aclose_f') - call h5tclose_f(type_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5tclose_f') - call h5sclose_f(space_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5sclose_f') + call h5screate_f(H5S_SCALAR_F,space_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5screate_f') + call h5tcopy_f(H5T_NATIVE_CHARACTER, type_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5tcopy_f') + call h5tset_size_f(type_id, int(len(trim(attrValue)),HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5tset_size_f') + call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5aexists_by_name_f') + if (attrExists) then + call h5adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5adelete_by_name_f') + endif + call h5acreate_by_name_f(loc_id,trim(p),trim(attrLabel),type_id,space_id,attr_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5acreate_f') + call h5awrite_f(attr_id, type_id, trim(attrValue), int([1],HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5awrite_f') + call h5aclose_f(attr_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5aclose_f') + call h5tclose_f(type_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5tclose_f') + call h5sclose_f(space_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5sclose_f') end subroutine HDF5_addAttribute_str @@ -336,117 +338,113 @@ end subroutine HDF5_addAttribute_str !-------------------------------------------------------------------------------------------------- !> @brief adds a integer attribute to the path given relative to the location !-------------------------------------------------------------------------------------------------- -subroutine HDF5_addAttribute_pInt(loc_id,attrLabel,attrValue,path) +subroutine HDF5_addAttribute_int(loc_id,attrLabel,attrValue,path) - implicit none - integer(HID_T), intent(in) :: loc_id - character(len=*), intent(in) :: attrLabel - integer(pInt), intent(in) :: attrValue - character(len=*), intent(in), optional :: path - integer :: hdferr - integer(HID_T) :: attr_id, space_id, type_id - logical :: attrExists - character(len=256) :: p + integer(HID_T), intent(in) :: loc_id + character(len=*), intent(in) :: attrLabel + integer(pInt), intent(in) :: attrValue + character(len=*), intent(in), optional :: path + integer :: hdferr + integer(HID_T) :: attr_id, space_id, type_id + logical :: attrExists + character(len=256) :: p + + if (present(path)) then + p = trim(path) + else + p = '.' + endif - if (present(path)) then - p = trim(path) - else - p = '.' - endif + call h5screate_f(H5S_SCALAR_F,space_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5screate_f') + call h5tcopy_f(H5T_NATIVE_INTEGER, type_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5tcopy_f') + call h5tset_size_f(type_id, 1_HSIZE_T, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5tset_size_f') + call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5aexists_by_name_f') + if (attrExists) then + call h5adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5adelete_by_name_f') + endif + call h5acreate_by_name_f(loc_id,trim(p),trim(attrLabel),type_id,space_id,attr_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5acreate_f') + call h5awrite_f(attr_id, type_id, attrValue, int([1],HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5awrite_f') + call h5aclose_f(attr_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5aclose_f') + call h5tclose_f(type_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5tclose_f') + call h5sclose_f(space_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5sclose_f') - call h5screate_f(H5S_SCALAR_F,space_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5screate_f') - call h5tcopy_f(H5T_NATIVE_INTEGER, type_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5tcopy_f') - call h5tset_size_f(type_id, 1_HSIZE_T, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5tset_size_f') - call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5aexists_by_name_f') - if (attrExists) then - call h5adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5adelete_by_name_f') - endif - call h5acreate_by_name_f(loc_id,trim(p),trim(attrLabel),type_id,space_id,attr_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5acreate_f') - call h5awrite_f(attr_id, type_id, attrValue, int([1],HSIZE_T), hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5awrite_f') - call h5aclose_f(attr_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5aclose_f') - call h5tclose_f(type_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5tclose_f') - call h5sclose_f(space_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pInt: h5sclose_f') - -end subroutine HDF5_addAttribute_pInt +end subroutine HDF5_addAttribute_int !-------------------------------------------------------------------------------------------------- !> @brief adds a integer attribute to the path given relative to the location !-------------------------------------------------------------------------------------------------- -subroutine HDF5_addAttribute_pReal(loc_id,attrLabel,attrValue,path) +subroutine HDF5_addAttribute_real(loc_id,attrLabel,attrValue,path) - implicit none - integer(HID_T), intent(in) :: loc_id - character(len=*), intent(in) :: attrLabel - real(pReal), intent(in) :: attrValue - character(len=*), intent(in), optional :: path - integer :: hdferr - integer(HID_T) :: attr_id, space_id, type_id - logical :: attrExists - character(len=256) :: p + integer(HID_T), intent(in) :: loc_id + character(len=*), intent(in) :: attrLabel + real(pReal), intent(in) :: attrValue + character(len=*), intent(in), optional :: path + integer :: hdferr + integer(HID_T) :: attr_id, space_id, type_id + logical :: attrExists + character(len=256) :: p + + if (present(path)) then + p = trim(path) + else + p = '.' + endif - if (present(path)) then - p = trim(path) - else - p = '.' - endif + call h5screate_f(H5S_SCALAR_F,space_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5screate_f') + call h5tcopy_f(H5T_NATIVE_DOUBLE, type_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5tcopy_f') + call h5tset_size_f(type_id, 8_HSIZE_T, hdferr) ! ToDo + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5tset_size_f') + call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5aexists_by_name_f') + if (attrExists) then + call h5adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5adelete_by_name_f') + endif + call h5acreate_by_name_f(loc_id,trim(p),trim(attrLabel),type_id,space_id,attr_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5acreate_f') + call h5awrite_f(attr_id, type_id, attrValue, int([1],HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5awrite_f') + call h5aclose_f(attr_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5aclose_f') + call h5tclose_f(type_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5tclose_f') + call h5sclose_f(space_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5sclose_f') - call h5screate_f(H5S_SCALAR_F,space_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5screate_f') - call h5tcopy_f(H5T_NATIVE_DOUBLE, type_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5tcopy_f') - call h5tset_size_f(type_id, 8_HSIZE_T, hdferr) ! ToDo - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5tset_size_f') - call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5aexists_by_name_f') - if (attrExists) then - call h5adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5adelete_by_name_f') - endif - call h5acreate_by_name_f(loc_id,trim(p),trim(attrLabel),type_id,space_id,attr_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5acreate_f') - call h5awrite_f(attr_id, type_id, attrValue, int([1],HSIZE_T), hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5awrite_f') - call h5aclose_f(attr_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5aclose_f') - call h5tclose_f(type_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5tclose_f') - call h5sclose_f(space_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_pReal: h5sclose_f') - -end subroutine HDF5_addAttribute_pReal +end subroutine HDF5_addAttribute_real !-------------------------------------------------------------------------------------------------- !> @brief set link to object in results file !-------------------------------------------------------------------------------------------------- subroutine HDF5_setLink(loc_id,target_name,link_name) - use hdf5 - implicit none - character(len=*), intent(in) :: target_name, link_name + character(len=*), intent(in) :: target_name, link_name integer(HID_T), intent(in) :: loc_id - integer :: hdferr - logical :: linkExists - - call h5lexists_f(loc_id, link_name,linkExists, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_setLink: h5lexists_soft_f ('//trim(link_name)//')') - if (linkExists) then - call h5ldelete_f(loc_id,link_name, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_setLink: h5ldelete_soft_f ('//trim(link_name)//')') - endif - call h5lcreate_soft_f(target_name, loc_id, link_name, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_setLink: h5lcreate_soft_f ('//trim(target_name)//' '//trim(link_name)//')') + integer :: hdferr + logical :: linkExists + + call h5lexists_f(loc_id, link_name,linkExists, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_setLink: h5lexists_soft_f ('//trim(link_name)//')') + if (linkExists) then + call h5ldelete_f(loc_id,link_name, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_setLink: h5ldelete_soft_f ('//trim(link_name)//')') + endif + call h5lcreate_soft_f(target_name, loc_id, link_name, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_setLink: h5lcreate_soft_f ('//trim(target_name)//' '//trim(link_name)//')') end subroutine HDF5_setLink @@ -454,1140 +452,1218 @@ end subroutine HDF5_setLink !-------------------------------------------------------------------------------------------------- !> @brief read dataset of type pReal with 1 dimension !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal1(loc_id,dataset,datasetName,parallel) +subroutine HDF5_read_real1(loc_id,dataset,datasetName,parallel) - implicit none - real(pReal), intent(inout), dimension(:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel - - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & ! ToDo: Fortran 2018 size(shape(A)) = rank(A) - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr - -!--------------------------------------------------------------------------------------------------- -! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - -!--------------------------------------------------------------------------------------------------- -! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif - - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal1: h5dread_f') - - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - -end subroutine HDF5_read_pReal1 - -!-------------------------------------------------------------------------------------------------- -!> @brief read dataset of type pReal with 2 dimensions -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal2(loc_id,dataset,datasetName,parallel) - - implicit none - real(pReal), intent(inout), dimension(:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel - - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr - -!--------------------------------------------------------------------------------------------------- -! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - -!--------------------------------------------------------------------------------------------------- -! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif - - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal2: h5dread_f') - - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - -end subroutine HDF5_read_pReal2 - -!-------------------------------------------------------------------------------------------------- -!> @brief read dataset of type pReal with 2 dimensions -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal3(loc_id,dataset,datasetName,parallel) - - implicit none - real(pReal), intent(inout), dimension(:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel - - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr - -!--------------------------------------------------------------------------------------------------- -! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - -!--------------------------------------------------------------------------------------------------- -! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif + real(pReal), intent(inout), dimension(:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal3: h5dread_f') + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & ! ToDo: Fortran 2018 size(shape(A)) = rank(A) + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) +!--------------------------------------------------------------------------------------------------- +! determine shape of dataset + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) -end subroutine HDF5_read_pReal3 +!--------------------------------------------------------------------------------------------------- +! initialize HDF5 data structures + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif + + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real1: h5dread_f') + + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) + +end subroutine HDF5_read_real1 + +!-------------------------------------------------------------------------------------------------- +!> @brief read dataset of type pReal with 2 dimensions +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_read_real2(loc_id,dataset,datasetName,parallel) + + real(pReal), intent(inout), dimension(:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel + + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr + +!--------------------------------------------------------------------------------------------------- +! determine shape of dataset + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + +!--------------------------------------------------------------------------------------------------- +! initialize HDF5 data structures + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif + + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real2: h5dread_f') + + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) + +end subroutine HDF5_read_real2 + +!-------------------------------------------------------------------------------------------------- +!> @brief read dataset of type pReal with 2 dimensions +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_read_real3(loc_id,dataset,datasetName,parallel) + + real(pReal), intent(inout), dimension(:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel + + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr + +!--------------------------------------------------------------------------------------------------- +! determine shape of dataset + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + +!--------------------------------------------------------------------------------------------------- +! initialize HDF5 data structures + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif + + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real3: h5dread_f') + + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) + +end subroutine HDF5_read_real3 !-------------------------------------------------------------------------------------------------- !> @brief read dataset of type pReal with 4 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal4(loc_id,dataset,datasetName,parallel) +subroutine HDF5_read_real4(loc_id,dataset,datasetName,parallel) - implicit none - real(pReal), intent(inout), dimension(:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + real(pReal), intent(inout), dimension(:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) !--------------------------------------------------------------------------------------------------- ! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif + + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real4: h5dread_f') + + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal4: h5dread_f') - - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - -end subroutine HDF5_read_pReal4 +end subroutine HDF5_read_real4 !-------------------------------------------------------------------------------------------------- !> @brief read dataset of type pReal with 5 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal5(loc_id,dataset,datasetName,parallel) +subroutine HDF5_read_real5(loc_id,dataset,datasetName,parallel) - implicit none - real(pReal), intent(inout), dimension(:,:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel - - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr + real(pReal), intent(inout), dimension(:,:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel + + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) !--------------------------------------------------------------------------------------------------- ! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif + + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real5: h5dread_f') + + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal5: h5dread_f') - - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - -end subroutine HDF5_read_pReal5 +end subroutine HDF5_read_real5 !-------------------------------------------------------------------------------------------------- !> @brief read dataset of type pReal with 6 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal6(loc_id,dataset,datasetName,parallel) +subroutine HDF5_read_real6(loc_id,dataset,datasetName,parallel) - implicit none - real(pReal), intent(inout), dimension(:,:,:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel - - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr + real(pReal), intent(inout), dimension(:,:,:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel + + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) !--------------------------------------------------------------------------------------------------- ! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal6: h5dread_f') + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real6: h5dread_f') - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) -end subroutine HDF5_read_pReal6 +end subroutine HDF5_read_real6 !-------------------------------------------------------------------------------------------------- !> @brief read dataset of type pReal with 7 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal7(loc_id,dataset,datasetName,parallel) +subroutine HDF5_read_real7(loc_id,dataset,datasetName,parallel) - implicit none - real(pReal), intent(inout), dimension(:,:,:,:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + real(pReal), intent(inout), dimension(:,:,:,:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) !--------------------------------------------------------------------------------------------------- ! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif + + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real7: h5dread_f') + + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal7: h5dread_f') - - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - -end subroutine HDF5_read_pReal7 +end subroutine HDF5_read_real7 !-------------------------------------------------------------------------------------------------- !> @brief read dataset of type pInt with 1 dimension !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt1(loc_id,dataset,datasetName,parallel) +subroutine HDF5_read_int1(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel - - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr + integer, intent(inout), dimension(:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel + + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) !--------------------------------------------------------------------------------------------------- ! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif - call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt1: h5dread_f') + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int1: h5dread_f') - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) -end subroutine HDF5_read_pInt1 +end subroutine HDF5_read_int1 !-------------------------------------------------------------------------------------------------- !> @brief read dataset of type pInt with 2 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt2(loc_id,dataset,datasetName,parallel) +subroutine HDF5_read_int2(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + integer, intent(inout), dimension(:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) !--------------------------------------------------------------------------------------------------- ! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif - - call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt2: h5dread_f') + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int2: h5dread_f') + + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) -end subroutine HDF5_read_pInt2 +end subroutine HDF5_read_int2 !-------------------------------------------------------------------------------------------------- !> @brief read dataset of type pInt with 3 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt3(loc_id,dataset,datasetName,parallel) +subroutine HDF5_read_int3(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + integer, intent(inout), dimension(:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) !--------------------------------------------------------------------------------------------------- ! initialize HDF5 data structures - if (present(parallel)) then + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif + + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int3: h5dread_f') + + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt3: h5dread_f') - - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - -end subroutine HDF5_read_pInt3 +end subroutine HDF5_read_int3 !-------------------------------------------------------------------------------------------------- !> @brief read dataset of type pInt withh 4 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt4(loc_id,dataset,datasetName,parallel) +subroutine HDF5_read_int4(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel - - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr + integer, intent(inout), dimension(:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel + + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) !--------------------------------------------------------------------------------------------------- ! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif + + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int4: h5dread_f') + + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt4: h5dread_f') - - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - -end subroutine HDF5_read_pInt4 +end subroutine HDF5_read_int4 !-------------------------------------------------------------------------------------------------- !> @brief read dataset of type pInt with 5 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt5(loc_id,dataset,datasetName,parallel) +subroutine HDF5_read_int5(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:,:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel - - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr + integer, intent(inout), dimension(:,:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel + + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) !--------------------------------------------------------------------------------------------------- ! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif - call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt5: h5dread_f') - - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int5: h5dread_f') + + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) -end subroutine HDF5_read_pInt5 +end subroutine HDF5_read_int5 !-------------------------------------------------------------------------------------------------- !> @brief read dataset of type pInt with 6 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt6(loc_id,dataset,datasetName,parallel) +subroutine HDF5_read_int6(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:,:,:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + integer, intent(inout), dimension(:,:,:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) !--------------------------------------------------------------------------------------------------- ! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif + + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int6: h5dread_f') + + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt6: h5dread_f') - - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - -end subroutine HDF5_read_pInt6 +end subroutine HDF5_read_int6 !-------------------------------------------------------------------------------------------------- !> @brief read dataset of type pInt with 7 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt7(loc_id,dataset,datasetName,parallel) +subroutine HDF5_read_int7(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:,:,:,:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel - - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) - integer :: hdferr + integer, intent(inout), dimension(:,:,:,:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel + + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer :: hdferr !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) !--------------------------------------------------------------------------------------------------- ! initialize HDF5 data structures - if (present(parallel)) then - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,parallel) - else - call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & - myStart, globalShape, loc_id,localShape,datasetName,.false.) - endif + if (present(parallel)) then + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,parallel) + else + call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, & + myStart, totalShape, loc_id,myShape,datasetName,.false.) + endif + + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,& + file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int7: h5dread_f') + + call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,globalShape, hdferr,& - file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt7: h5dread_f') - - call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - -end subroutine HDF5_read_pInt7 +end subroutine HDF5_read_int7 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pReal with 1 dimension !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal1(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_real1(loc_id,dataset,datasetName,parallel) - implicit none - real(pReal), intent(inout), dimension(:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + real(pReal), intent(inout), dimension(:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape,loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape,loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,.false.) - endif + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,.false.) + endif - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5dread_f') + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real1: h5dwrite_f') + endif - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) -end subroutine HDF5_write_pReal1 +end subroutine HDF5_write_real1 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pReal with 2 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal2(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_real2(loc_id,dataset,datasetName,parallel) - implicit none - real(pReal), intent(inout), dimension(:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + real(pReal), intent(inout), dimension(:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,.false.) - endif + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,.false.) + endif - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5dread_f') + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real2: h5dwrite_f') + endif - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) -end subroutine HDF5_write_pReal2 +end subroutine HDF5_write_real2 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pReal with 3 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal3(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_real3(loc_id,dataset,datasetName,parallel) - implicit none - real(pReal), intent(inout), dimension(:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + real(pReal), intent(inout), dimension(:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,.false.) - endif + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,.false.) + endif - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5dread_f') + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real3: h5dwrite_f') + endif - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) -end subroutine HDF5_write_pReal3 +end subroutine HDF5_write_real3 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pReal with 4 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal4(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_real4(loc_id,dataset,datasetName,parallel) - implicit none - real(pReal), intent(inout), dimension(:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + real(pReal), intent(inout), dimension(:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,.false.) + endif + + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real4: h5dwrite_f') + endif + + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,.false.) - endif - - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5dread_f') - - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) - -end subroutine HDF5_write_pReal4 +end subroutine HDF5_write_real4 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pReal with 5 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal5(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_real5(loc_id,dataset,datasetName,parallel) - implicit none - real(pReal), intent(inout), dimension(:,:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + real(pReal), intent(inout), dimension(:,:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,.false.) + endif + + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real5: h5dwrite_f') + endif + + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,.false.) - endif - - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5dread_f') - - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) - -end subroutine HDF5_write_pReal5 +end subroutine HDF5_write_real5 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pReal with 6 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal6(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_real6(loc_id,dataset,datasetName,parallel) - implicit none - real(pReal), intent(inout), dimension(:,:,:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + real(pReal), intent(inout), dimension(:,:,:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,.false.) - endif + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,.false.) + endif - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5dread_f') + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real6: h5dwrite_f') + endif + + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) - -end subroutine HDF5_write_pReal6 +end subroutine HDF5_write_real6 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pReal with 7 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal7(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_real7(loc_id,dataset,datasetName,parallel) - implicit none - real(pReal), intent(inout), dimension(:,:,:,:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + real(pReal), intent(inout), dimension(:,:,:,:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_DOUBLE,.false.) - endif + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,.false.) + endif - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5dread_f') + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real7: h5dwrite_f') + endif + + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) - -end subroutine HDF5_write_pReal7 +end subroutine HDF5_write_real7 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pInt with 1 dimension !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt1(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_int1(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + integer, intent(inout), dimension(:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,.false.) - endif + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,.false.) + endif - call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pInt1: h5dread_f') + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int1: h5dwrite_f') + endif + + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) - -end subroutine HDF5_write_pInt1 +end subroutine HDF5_write_int1 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pInt with 2 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt2(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_int2(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + integer, intent(inout), dimension(:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,.false.) - endif + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,.false.) + endif - call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pInt2: h5dread_f') + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int2: h5dwrite_f') + endif - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) -end subroutine HDF5_write_pInt2 +end subroutine HDF5_write_int2 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pInt with 3 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt3(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_int3(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + integer, intent(inout), dimension(:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,.false.) - endif + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,.false.) + endif - call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pInt3: h5dread_f') + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int3: h5dwrite_f') + endif + + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) - -end subroutine HDF5_write_pInt3 +end subroutine HDF5_write_int3 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pInt with 4 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt4(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_int4(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + integer, intent(inout), dimension(:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,.false.) - endif + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,.false.) + endif - call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pInt4: h5dread_f') + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int4: h5dwrite_f') + endif - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) -end subroutine HDF5_write_pInt4 +end subroutine HDF5_write_int4 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pInt with 5 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt5(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_int5(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:,:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + integer, intent(inout), dimension(:,:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,.false.) - endif + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,.false.) + endif - call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pInt5: h5dread_f') + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int5: h5dwrite_f') + endif - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) -end subroutine HDF5_write_pInt5 +end subroutine HDF5_write_int5 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pInt with 6 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt6(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_int6(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:,:,:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + integer, intent(inout), dimension(:,:,:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - if (present(parallel)) then - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,.false.) - endif + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,.false.) + endif - call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pInt6: h5dread_f') + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int6: h5dwrite_f') + endif - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) -end subroutine HDF5_write_pInt6 +end subroutine HDF5_write_int6 !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type pInt with 7 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt7(loc_id,dataset,datasetName,parallel) +subroutine HDF5_write_int7(loc_id,dataset,datasetName,parallel) - implicit none - integer(pInt), intent(inout), dimension(:,:,:,:,:,:,:) :: dataset - integer(HID_T), intent(in) :: loc_id !< file or group handle - character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in), optional :: parallel + integer, intent(inout), dimension(:,:,:,:,:,:,:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id - integer(HSIZE_T), dimension(size(shape(dataset))) :: & - myStart, & - localShape, & !< shape of the dataset (this process) - globalShape !< shape of the dataset (all processes) + integer :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) !--------------------------------------------------------------------------------------------------- ! determine shape of dataset - localShape = int(shape(dataset),HSIZE_T) - if (any(localShape(1:size(localShape)-1) == 0)) return !< empty dataset (last dimension can be empty) + myShape = int(shape(dataset),HSIZE_T) + if (any(myShape(1:size(myShape)-1) == 0)) return !< empty dataset (last dimension can be empty) - if (present(parallel)) then + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel) + else + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,.false.) + endif + + if (product(totalShape) /= 0) then + call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int7: h5dwrite_f') + endif + + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) + +end subroutine HDF5_write_int7 + + +!-------------------------------------------------------------------------------------------------- +!> @brief writes a scalar orientation dataset +! ToDo: It might be possible to write the dataset as a whole +! ToDo: We could optionally write out other representations (axis angle, euler, ...) +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_write_rotation(loc_id,dataset,datasetName,parallel) + use rotations, only: & + rotation + + type(rotation), intent(in), dimension(:) :: dataset + integer(HID_T), intent(in) :: loc_id !< file or group handle + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + logical, intent(in), optional :: parallel + + integer :: hdferr + real(pReal), dimension(4,size(dataset)) :: dataset_asArray + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id,dtype_id,w_id,x_id,y_id,z_id + integer(HSIZE_T), dimension(size(shape(dataset))) :: & + myStart, & + myShape, & !< shape of the dataset (this process) + totalShape !< shape of the dataset (all processes) + integer(SIZE_T) :: type_size_real + integer :: i + + do i = 1, size(dataset) + dataset_asArray(1:4,i) = dataset(i)%asQuaternion() + enddo + +!--------------------------------------------------------------------------------------------------- +! determine shape of dataset + myShape = int(shape(dataset),HSIZE_T) + +!--------------------------------------------------------------------------------------------------- +! compound type: name of each quaternion component + call h5tget_size_f(H5T_NATIVE_DOUBLE, type_size_real, hdferr) + + call h5tcreate_f(H5T_COMPOUND_F, type_size_real*4_SIZE_T, dtype_id, hdferr) + call h5tinsert_f(dtype_id, "w", type_size_real*0_SIZE_T, H5T_NATIVE_DOUBLE, hdferr) + call h5tinsert_f(dtype_id, "x", type_size_real*1_SIZE_T, H5T_NATIVE_DOUBLE, hdferr) + call h5tinsert_f(dtype_id, "y", type_size_real*2_SIZE_T, H5T_NATIVE_DOUBLE, hdferr) + call h5tinsert_f(dtype_id, "z", type_size_real*3_SIZE_T, H5T_NATIVE_DOUBLE, hdferr) + + if (present(parallel)) then + call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & + myStart, totalShape, loc_id,myShape,datasetName,dtype_id,parallel) + else call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,parallel) - else - call initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, loc_id,localShape,datasetName,H5T_NATIVE_INTEGER,.false.) - endif + myStart, totalShape, loc_id,myShape,datasetName,dtype_id,.false.) + endif - call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), hdferr,& - file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pInt7: h5dread_f') + call h5pset_preserve_f(plist_id, .TRUE., hdferr) + + if (product(totalShape) /= 0) then + call h5tcreate_f(H5T_COMPOUND_F, type_size_real, x_id, hdferr) + call h5tinsert_f(x_id, "x", 0_SIZE_T, H5T_NATIVE_DOUBLE, hdferr) + call h5tcreate_f(H5T_COMPOUND_F, type_size_real, w_id, hdferr) + call h5tinsert_f(w_id, "w", 0_SIZE_T, H5T_NATIVE_DOUBLE, hdferr) + call h5tcreate_f(H5T_COMPOUND_F, type_size_real, y_id, hdferr) + call h5tinsert_f(y_id, "y", 0_SIZE_T, H5T_NATIVE_DOUBLE, hdferr) + call h5tcreate_f(H5T_COMPOUND_F, type_size_real, z_id, hdferr) + call h5tinsert_f(z_id, "z", 0_SIZE_T, H5T_NATIVE_DOUBLE, hdferr) + + call h5dwrite_f(dset_id, w_id,dataset_asArray(1,:),int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + call h5dwrite_f(dset_id, x_id,dataset_asArray(2,:),int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + call h5dwrite_f(dset_id, y_id,dataset_asArray(3,:),int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + call h5dwrite_f(dset_id, z_id,dataset_asArray(4,:),int(totalShape,HSIZE_T), hdferr,& + file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_rotation: h5dwrite_f') + endif - call finalize_write(plist_id, dset_id, filespace_id, memspace_id) + call finalize_write(plist_id, dset_id, filespace_id, memspace_id) -end subroutine HDF5_write_pInt7 +end subroutine HDF5_write_rotation !-------------------------------------------------------------------------------------------------- @@ -1600,7 +1676,6 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_ worldrank, & worldsize - implicit none integer(HID_T), intent(in) :: loc_id !< file or group handle character(len=*), intent(in) :: datasetName !< name of the dataset in the file logical, intent(in) :: parallel @@ -1670,12 +1745,13 @@ end subroutine initialize_read !-------------------------------------------------------------------------------------------------- subroutine finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id) - implicit none integer(HID_T), intent(in) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id integer :: hdferr call h5pclose_f(plist_id, hdferr) if (hdferr < 0) call IO_error(1_pInt,ext_msg='finalize_read: plist_id') + call h5pclose_f(aplist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='finalize_read: aplist_id') call h5dclose_f(dset_id, hdferr) if (hdferr < 0) call IO_error(1_pInt,ext_msg='finalize_read: h5dclose_f') call h5sclose_f(filespace_id, hdferr) @@ -1690,68 +1766,68 @@ end subroutine finalize_read !> @brief initialize HDF5 handles, determines global shape and start for parallel write !-------------------------------------------------------------------------------------------------- subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, & - myStart, globalShape, & - loc_id,localShape,datasetName,datatype,parallel) + myStart, totalShape, & + loc_id,myShape,datasetName,datatype,parallel) use numerics, only: & worldrank, & worldsize - implicit none integer(HID_T), intent(in) :: loc_id !< file or group handle character(len=*), intent(in) :: datasetName !< name of the dataset in the file - logical, intent(in) :: parallel - integer(HID_T), intent(in) :: datatype + logical, intent(in) :: parallel + integer(HID_T), intent(in) :: datatype integer(HSIZE_T), intent(in), dimension(:) :: & - localShape - integer(HSIZE_T), intent(out), dimension(size(localShape,1)):: & + myShape + integer(HSIZE_T), intent(out), dimension(size(myShape,1)):: & myStart, & - globalShape !< shape of the dataset (all processes) + totalShape !< shape of the dataset (all processes) integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id - integer(pInt), dimension(worldsize) :: & + integer, dimension(worldsize) :: & writeSize !< contribution of all processes integer :: ierr integer :: hdferr !------------------------------------------------------------------------------------------------- -! creating a property list for transfer properties +! creating a property list for transfer properties (is collective when reading in parallel) call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_write: h5pcreate_f') - -!-------------------------------------------------------------------------------------------------- - writeSize = 0_pInt - writeSize(worldrank+1) = int(localShape(ubound(localShape,1)),pInt) - #ifdef PETSc if (parallel) then call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_write: h5pset_dxpl_mpio_f') + endif +#endif + +!-------------------------------------------------------------------------------------------------- +! determine the global data layout among all processes + writeSize = 0_pInt + writeSize(worldrank+1) = int(myShape(ubound(myShape,1)),pInt) +#ifdef PETSc +if (parallel) then call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process if (ierr /= 0) call IO_error(894_pInt,ext_msg='initialize_write: MPI_allreduce') endif #endif - myStart = int(0,HSIZE_T) myStart(ubound(myStart)) = int(sum(writeSize(1:worldrank)),HSIZE_T) - globalShape = [localShape(1:ubound(localShape,1)-1),int(sum(writeSize),HSIZE_T)] + totalShape = [myShape(1:ubound(myShape,1)-1),int(sum(writeSize),HSIZE_T)] !-------------------------------------------------------------------------------------------------- ! create dataspace in memory (local shape) and in file (global shape) - call h5screate_simple_f(size(localShape), localShape, memspace_id, hdferr, localShape) + call h5screate_simple_f(size(myShape), myShape, memspace_id, hdferr, myShape) if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_write: h5dopen_f') - call h5screate_simple_f(size(globalShape), globalShape, filespace_id, hdferr, globalShape) + call h5screate_simple_f(size(totalShape), totalShape, filespace_id, hdferr, totalShape) if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_write: h5dget_space_f') !-------------------------------------------------------------------------------------------------- -! create dataset +! create dataset in the file and select a hyperslab from it (the portion of the current process) call h5dcreate_f(loc_id, trim(datasetName), datatype, filespace_id, dset_id, hdferr) if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_write: h5dcreate_f') - -!-------------------------------------------------------------------------------------------------- -! select a hyperslab (the portion of the current process) in the file - call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, localShape, hdferr) + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, myShape, hdferr) if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_write: h5sselect_hyperslab_f') + end subroutine initialize_write @@ -1760,7 +1836,6 @@ end subroutine initialize_write !-------------------------------------------------------------------------------------------------- subroutine finalize_write(plist_id, dset_id, filespace_id, memspace_id) - implicit none integer(HID_T), intent(in) :: dset_id, filespace_id, memspace_id, plist_id integer :: hdferr diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 5131eeaa9..0ae1323f0 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -9,10 +9,6 @@ #include "list.f90" #include "future.f90" #include "config.f90" -#ifdef DAMASKHDF5 -#include "HDF5_utilities.f90" -#include "results.f90" -#endif #include "math.f90" #include "quaternions.f90" #include "Lambert.f90" @@ -26,6 +22,10 @@ #ifdef Marc4DAMASK #include "mesh_marc.f90" #endif +#ifdef DAMASK_HDF5 +#include "HDF5_utilities.f90" +#include "results.f90" +#endif #include "material.f90" #include "lattice.f90" #include "source_thermal_dissipation.f90" diff --git a/src/config.f90 b/src/config.f90 index 23268d1de..f86057b25 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -231,15 +231,21 @@ end function read_materialConfig !-------------------------------------------------------------------------------------------------- subroutine parse_materialConfig(sectionNames,part,line, & fileContent) + use prec, only: & + pStringLen + use IO, only: & + IO_intOut + implicit none character(len=64), allocatable, dimension(:), intent(out) :: sectionNames type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part character(len=pStringLen), intent(inout) :: line character(len=pStringLen), dimension(:), intent(in) :: fileContent - integer, allocatable, dimension(:) :: partPosition ! position of [] tags + last line in section - integer :: i, j - logical :: echo + integer, allocatable, dimension(:) :: partPosition !< position of [] tags + last line in section + integer :: i, j + logical :: echo + character(len=pStringLen) :: section_ID echo = .false. @@ -263,7 +269,8 @@ subroutine parse_materialConfig(sectionNames,part,line, & partPosition = [partPosition, i] ! needed when actually storing content do i = 1, size(partPosition) -1 - sectionNames(i) = trim(adjustl(IO_getTag(fileContent(partPosition(i)),'[',']'))) + write(section_ID,'('//IO_intOut(size(partPosition))//')') i + sectionNames(i) = trim(section_ID)//'_'//trim(adjustl(IO_getTag(fileContent(partPosition(i)),'[',']'))) do j = partPosition(i) + 1, partPosition(i+1) -1 call part(i)%add(trim(adjustl(fileContent(j)))) enddo diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 23ae3f07b..1158ddc07 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -9,7 +9,7 @@ module constitutive implicit none private - integer(pInt), public, protected :: & + integer, public, protected :: & constitutive_plasticity_maxSizePostResults, & constitutive_plasticity_maxSizeDotState, & constitutive_source_maxSizePostResults, & @@ -37,7 +37,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates arrays pointing to array of the various constitutive modules !-------------------------------------------------------------------------------------------------- -subroutine constitutive_init() +subroutine constitutive_init use prec, only: & pReal use debug, only: & @@ -50,8 +50,7 @@ subroutine constitutive_init() IO_write_jobFile use config, only: & material_Nphase, & - phase_name, & - config_deallocate + phase_name use material, only: & material_phase, & phase_plasticity, & @@ -111,14 +110,14 @@ subroutine constitutive_init() use kinematics_thermal_expansion implicit none - integer(pInt), parameter :: FILEUNIT = 204_pInt - integer(pInt) :: & + integer, parameter :: FILEUNIT = 204 + integer :: & o, & !< counter in output loop ph, & !< counter in phase loop s, & !< counter in source loop ins !< instance of plasticity/source - integer(pInt), dimension(:,:), pointer :: thisSize + integer, dimension(:,:), pointer :: thisSize character(len=64), dimension(:,:), pointer :: thisOutput character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready logical :: knownPlasticity, knownSource, nonlocalConstitutionPresent @@ -149,15 +148,13 @@ subroutine constitutive_init() if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init - call config_deallocate('material.config/phase') - write(6,'(/,a)') ' <<<+- constitutive init -+>>>' mainProcess: if (worldrank == 0) then !-------------------------------------------------------------------------------------------------- ! write description file for constitutive output call IO_write_jobFile(FILEUNIT,'outputConstitutive') - PhaseLoop: do ph = 1_pInt,material_Nphase + PhaseLoop: do ph = 1,material_Nphase activePhase: if (any(material_phase == ph)) then ins = phase_plasticityInstance(ph) knownPlasticity = .true. ! assume valid @@ -197,14 +194,14 @@ subroutine constitutive_init() if (knownPlasticity) then write(FILEUNIT,'(a)') '(plasticity)'//char(9)//trim(outputName) if (phase_plasticity(ph) /= PLASTICITY_NONE_ID) then - OutputPlasticityLoop: do o = 1_pInt,size(thisOutput(:,ins)) - if(len(trim(thisOutput(o,ins))) > 0_pInt) & + OutputPlasticityLoop: do o = 1,size(thisOutput(:,ins)) + if(len(trim(thisOutput(o,ins))) > 0) & write(FILEUNIT,'(a,i4)') trim(thisOutput(o,ins))//char(9),thisSize(o,ins) enddo OutputPlasticityLoop endif endif - SourceLoop: do s = 1_pInt, phase_Nsources(ph) + SourceLoop: do s = 1, phase_Nsources(ph) knownSource = .true. ! assume valid sourceType: select case (phase_source(s,ph)) case (SOURCE_thermal_dissipation_ID) sourceType @@ -242,8 +239,8 @@ subroutine constitutive_init() end select sourceType if (knownSource) then write(FILEUNIT,'(a)') '(source)'//char(9)//trim(outputName) - OutputSourceLoop: do o = 1_pInt,size(thisOutput(:,ins)) - if(len(trim(thisOutput(o,ins))) > 0_pInt) & + OutputSourceLoop: do o = 1,size(thisOutput(:,ins)) + if(len(trim(thisOutput(o,ins))) > 0) & write(FILEUNIT,'(a,i4)') trim(thisOutput(o,ins))//char(9),thisSize(o,ins) enddo OutputSourceLoop endif @@ -253,17 +250,17 @@ subroutine constitutive_init() close(FILEUNIT) endif mainProcess - constitutive_plasticity_maxSizeDotState = 0_pInt - constitutive_plasticity_maxSizePostResults = 0_pInt - constitutive_source_maxSizeDotState = 0_pInt - constitutive_source_maxSizePostResults = 0_pInt + constitutive_plasticity_maxSizeDotState = 0 + constitutive_plasticity_maxSizePostResults = 0 + constitutive_source_maxSizeDotState = 0 + constitutive_source_maxSizePostResults = 0 - PhaseLoop2:do ph = 1_pInt,material_Nphase + PhaseLoop2:do ph = 1,material_Nphase !-------------------------------------------------------------------------------------------------- ! partition and inititalize state plasticState(ph)%partionedState0 = plasticState(ph)%state0 plasticState(ph)%state = plasticState(ph)%partionedState0 - forall(s = 1_pInt:phase_Nsources(ph)) + forall(s = 1:phase_Nsources(ph)) sourceState(ph)%p(s)%partionedState0 = sourceState(ph)%p(s)%state0 sourceState(ph)%p(s)%state = sourceState(ph)%p(s)%partionedState0 end forall @@ -302,7 +299,7 @@ function constitutive_homogenizedC(ipc,ip,el) implicit none real(pReal), dimension(6,6) :: constitutive_homogenizedC - integer(pInt), intent(in) :: & + integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element @@ -341,14 +338,14 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) plastic_disloUCLA_dependentState implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & Fe, & !< elastic deformation gradient Fp !< plastic deformation gradient - integer(pInt) :: & + integer :: & ho, & !< homogenization tme, & !< thermal member position instance, of @@ -412,7 +409,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & plastic_nonlocal_LpAndItsTangent implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element @@ -428,10 +425,10 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & dLp_dMp !< derivative of Lp with respect to Mandel stress real(pReal), dimension(3,3) :: & Mp !< Mandel stress work conjugate with Lp - integer(pInt) :: & + integer :: & ho, & !< homogenization tme !< thermal member position - integer(pInt) :: & + integer :: & i, j, instance, of ho = material_homogenizationAt(el) @@ -519,7 +516,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & kinematics_thermal_expansion_LiAndItsTangent implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element @@ -541,7 +538,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & my_dLi_dS real(pReal) :: & detFi - integer(pInt) :: & + integer :: & k, i, j, & instance, of @@ -562,7 +559,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS - KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el)) + KinematicsLoop: do k = 1, phase_Nkinematics(material_phase(ipc,ip,el)) kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el))) case (KINEMATICS_cleavage_opening_ID) kinematicsType call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ipc, ip, el) @@ -583,7 +580,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration temp_33 = matmul(FiInv,Li) - do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt + do i = 1,3; do j = 1,3 dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) @@ -612,22 +609,22 @@ pure function constitutive_initialFi(ipc, ip, el) kinematics_thermal_expansion_initialStrain implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), dimension(3,3) :: & constitutive_initialFi !< composite initial intermediate deformation gradient - integer(pInt) :: & + integer :: & k !< counter in kinematics loop - integer(pInt) :: & + integer :: & phase, & homog, offset constitutive_initialFi = math_I3 phase = material_phase(ipc,ip,el) - KinematicsLoop: do k = 1_pInt, phase_Nkinematics(phase) !< Warning: small initial strain assumption + KinematicsLoop: do k = 1, phase_Nkinematics(phase) !< Warning: small initial strain assumption kinematicsType: select case (phase_kinematics(k,phase)) case (KINEMATICS_thermal_expansion_ID) kinematicsType homog = material_homogenizationAt(el) @@ -650,7 +647,7 @@ subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) pReal implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element @@ -691,7 +688,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & STIFFNESS_DEGRADATION_damage_ID implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element @@ -705,19 +702,19 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient real(pReal), dimension(3,3) :: E real(pReal), dimension(3,3,3,3) :: C - integer(pInt) :: & + integer :: & ho, & !< homogenization d !< counter in degradation loop - integer(pInt) :: & + integer :: & i, j ho = material_homogenizationAt(el) C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el)) - DegradationLoop: do d = 1_pInt, phase_NstiffnessDegradations(material_phase(ipc,ip,el)) + DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phase(ipc,ip,el)) degradationType: select case(phase_stiffnessDegradation(d,material_phase(ipc,ip,el))) case (STIFFNESS_DEGRADATION_damage_ID) degradationType - C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2_pInt + C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2 end select degradationType enddo DegradationLoop @@ -725,7 +722,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration dS_dFe = 0.0_pReal - forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt) + forall (i=1:3, j=1:3) dS_dFe(i,j,1:3,1:3) = & matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn @@ -790,7 +787,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, source_thermal_externalheat_dotState implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element @@ -805,7 +802,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, S !< 2nd Piola Kirchhoff stress (vector notation) real(pReal), dimension(3,3) :: & Mp - integer(pInt) :: & + integer :: & ho, & !< homogenization tme, & !< thermal member position i, & !< counter in source loop @@ -848,7 +845,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, subdt,ip,el) end select plasticityType - SourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) + SourceLoop: do i = 1, phase_Nsources(material_phase(ipc,ip,el)) sourceType: select case (phase_source(i,material_phase(ipc,ip,el))) @@ -900,7 +897,7 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) source_damage_isoBrittle_deltaState implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element @@ -910,7 +907,7 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) Fi !< intermediate deformation gradient real(pReal), dimension(3,3) :: & Mp - integer(pInt) :: & + integer :: & i, & instance, of @@ -928,7 +925,7 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) end select plasticityType - sourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) + sourceLoop: do i = 1, phase_Nsources(material_phase(ipc,ip,el)) sourceType: select case (phase_source(i,material_phase(ipc,ip,el))) @@ -994,7 +991,7 @@ function constitutive_postResults(S, Fi, ipc, ip, el) source_damage_anisoDuctile_postResults implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element @@ -1007,9 +1004,9 @@ function constitutive_postResults(S, Fi, ipc, ip, el) S !< 2nd Piola Kirchhoff stress real(pReal), dimension(3,3) :: & Mp !< Mandel stress - integer(pInt) :: & + integer :: & startPos, endPos - integer(pInt) :: & + integer :: & ho, & !< homogenization tme, & !< thermal member position i, of, instance !< counter in source loop @@ -1021,7 +1018,7 @@ function constitutive_postResults(S, Fi, ipc, ip, el) ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) - startPos = 1_pInt + startPos = 1 endPos = plasticState(material_phase(ipc,ip,el))%sizePostResults of = phasememberAt(ipc,ip,el) @@ -1054,8 +1051,8 @@ function constitutive_postResults(S, Fi, ipc, ip, el) end select plasticityType - SourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) - startPos = endPos + 1_pInt + SourceLoop: do i = 1, phase_Nsources(material_phase(ipc,ip,el)) + startPos = endPos + 1 endPos = endPos + sourceState(material_phase(ipc,ip,el))%p(i)%sizePostResults of = phasememberAt(ipc,ip,el) sourceType: select case (phase_source(i,material_phase(ipc,ip,el))) @@ -1077,7 +1074,7 @@ end function constitutive_postResults !-------------------------------------------------------------------------------------------------- !> @brief writes constitutive results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine constitutive_results() +subroutine constitutive_results use material, only: & PLASTICITY_ISOTROPIC_ID, & PLASTICITY_PHENOPOWERLAW_ID, & @@ -1085,7 +1082,7 @@ subroutine constitutive_results() PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOUCLA_ID, & PLASTICITY_NONLOCAL_ID -#if defined(PETSc) || defined(DAMASKHDF5) +#if defined(PETSc) || defined(DAMASK_HDF5) use results use HDF5_utilities use config, only: & @@ -1108,36 +1105,40 @@ subroutine constitutive_results() use plastic_nonlocal, only: & plastic_nonlocal_results - implicit none - integer :: p - call HDF5_closeGroup(results_addGroup('current/phase')) - do p=1,size(config_name_phase) - call HDF5_closeGroup(results_addGroup('current/phase/'//trim(config_name_phase(p)))) + integer :: p + character(len=256) :: group + + call HDF5_closeGroup(results_addGroup('current/constitutive')) + + do p=1,size(config_name_phase) + group = trim('current/constitutive')//'/'//trim(config_name_phase(p)) + call HDF5_closeGroup(results_addGroup(group)) + group = trim(group)//'/'//'plastic' + + call HDF5_closeGroup(results_addGroup(group)) select case(material_phase_plasticity_type(p)) case(PLASTICITY_ISOTROPIC_ID) - call plastic_isotropic_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p))) + call plastic_isotropic_results(phase_plasticityInstance(p),group) case(PLASTICITY_PHENOPOWERLAW_ID) - call plastic_phenopowerlaw_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p))) + call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group) - case(PLASTICITY_KINEHARDENING_ID) - call plastic_kinehardening_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p))) + case(PLASTICITY_KINEHARDENING_ID) + call plastic_kinehardening_results(phase_plasticityInstance(p),group) - case(PLASTICITY_DISLOTWIN_ID) - call plastic_dislotwin_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p))) + case(PLASTICITY_DISLOTWIN_ID) + call plastic_dislotwin_results(phase_plasticityInstance(p),group) - case(PLASTICITY_DISLOUCLA_ID) - call plastic_disloUCLA_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p))) + case(PLASTICITY_DISLOUCLA_ID) + call plastic_disloUCLA_results(phase_plasticityInstance(p),group) - case(PLASTICITY_NONLOCAL_ID) - call plastic_nonlocal_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p))) - + case(PLASTICITY_NONLOCAL_ID) + call plastic_nonlocal_results(phase_plasticityInstance(p),group) end select - enddo - + enddo #endif diff --git a/src/crystallite.f90 b/src/crystallite.f90 index b234d133d..69c7839c7 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -10,7 +10,8 @@ module crystallite use prec, only: & - pReal + pReal, & + pStringLen use rotations, only: & rotation use FEsolving, only: & @@ -103,6 +104,13 @@ module crystallite end enum integer(kind(undefined_ID)),dimension(:,:), allocatable, private :: & crystallite_outputID !< ID of each post result output + + type, private :: tOutput !< new requested output (per phase) + character(len=65536), allocatable, dimension(:) :: & + label + end type tOutput + type(tOutput), allocatable, dimension(:), private :: output_constituent + procedure(), pointer :: integrateState public :: & @@ -111,7 +119,8 @@ module crystallite crystallite_stressTangent, & crystallite_orientations, & crystallite_push33ToRef, & - crystallite_postResults + crystallite_postResults, & + crystallite_results private :: & integrateStress, & integrateState, & @@ -156,6 +165,7 @@ subroutine crystallite_init use config, only: & config_deallocate, & config_crystallite, & + config_phase, & crystallite_name use constitutive, only: & constitutive_initialFi, & @@ -296,6 +306,18 @@ subroutine crystallite_init end select outputName enddo enddo + + allocate(output_constituent(size(config_phase))) + do c = 1, size(config_phase) +#if defined(__GFORTRAN__) + allocate(output_constituent(c)%label(1)) + output_constituent(c)%label(1)= 'GfortranBug86277' + output_constituent(c)%label = config_phase(c)%getStrings('(output)',defaultVal=output_constituent(c)%label ) + if (output_constituent(c)%label (1) == 'GfortranBug86277') output_constituent(c)%label = [character(len=pStringLen)::] +#else + output_constituent(c)%label = config_phase(c)%getStrings('(output)',defaultVal=[character(len=pStringLen)::]) +#endif + enddo do r = 1,size(config_crystallite) @@ -340,6 +362,7 @@ subroutine crystallite_init close(FILEUNIT) endif + call config_deallocate('material.config/phase') call config_deallocate('material.config/crystallite') !-------------------------------------------------------------------------------------------------- @@ -1053,6 +1076,158 @@ function crystallite_postResults(ipc, ip, el) end function crystallite_postResults +!-------------------------------------------------------------------------------------------------- +!> @brief writes constitutive results to HDF5 output file +!-------------------------------------------------------------------------------------------------- +subroutine crystallite_results +#if defined(PETSc) || defined(DAMASK_HDF5) + use lattice + use results + use HDF5_utilities + use rotations + use config, only: & + config_name_phase => phase_name ! anticipate logical name + + use material, only: & + material_phase_plasticity_type => phase_plasticity + + implicit none + integer :: p,o + real(pReal), allocatable, dimension(:,:,:) :: selected_tensors + type(rotation), allocatable, dimension(:) :: selected_rotations + character(len=256) :: group,lattice_label + + call HDF5_closeGroup(results_addGroup('current/constituent')) + + do p=1,size(config_name_phase) + group = trim('current/constituent')//'/'//trim(config_name_phase(p)) + call HDF5_closeGroup(results_addGroup(group)) + do o = 1, size(output_constituent(p)%label) + select case (output_constituent(p)%label(o)) + case('f') + selected_tensors = select_tensors(crystallite_partionedF,p) + call results_writeDataset(group,selected_tensors,'F',& + 'deformation gradient','1') + case('fe') + selected_tensors = select_tensors(crystallite_Fe,p) + call results_writeDataset(group,selected_tensors,'Fe',& + 'elastic deformation gradient','1') + case('fp') + selected_tensors = select_tensors(crystallite_Fp,p) + call results_writeDataset(group,selected_tensors,'Fp',& + 'plastic deformation gradient','1') + case('fi') + selected_tensors = select_tensors(crystallite_Fi,p) + call results_writeDataset(group,selected_tensors,'Fi',& + 'inelastic deformation gradient','1') + case('lp') + selected_tensors = select_tensors(crystallite_Lp,p) + call results_writeDataset(group,selected_tensors,'Lp',& + 'plastic velocity gradient','1/s') + case('li') + selected_tensors = select_tensors(crystallite_Li,p) + call results_writeDataset(group,selected_tensors,'Li',& + 'inelastic velocity gradient','1/s') + case('p') + selected_tensors = select_tensors(crystallite_P,p) + call results_writeDataset(group,selected_tensors,'P',& + '1st Piola-Kirchoff stress','Pa') + case('s') + selected_tensors = select_tensors(crystallite_S,p) + call results_writeDataset(group,selected_tensors,'S',& + '2nd Piola-Kirchoff stress','Pa') + case('orientation') + select case(lattice_structure(p)) + case(LATTICE_iso_ID) + lattice_label = 'iso' + case(LATTICE_fcc_ID) + lattice_label = 'fcc' + case(LATTICE_bcc_ID) + lattice_label = 'bcc' + case(LATTICE_bct_ID) + lattice_label = 'bct' + case(LATTICE_hex_ID) + lattice_label = 'hex' + case(LATTICE_ort_ID) + lattice_label = 'ort' + end select + selected_rotations = select_rotations(crystallite_orientation,p) + call results_writeDataset(group,selected_rotations,'orientation',& + 'crystal orientation as quaternion',lattice_label) + end select + enddo + enddo + + contains + +!-------------------------------------------------------------------------------------------------- +!> @brief select tensors for output +!-------------------------------------------------------------------------------------------------- + function select_tensors(dataset,instance) + + use material, only: & + homogenization_maxNgrains, & + material_phaseAt + + integer, intent(in) :: instance + real(pReal), dimension(:,:,:,:,:), intent(in) :: dataset + real(pReal), allocatable, dimension(:,:,:) :: select_tensors + integer :: e,i,c,j + + allocate(select_tensors(3,3,count(material_phaseAt==instance)*homogenization_maxNgrains)) + + j=1 + do e = 1, size(material_phaseAt,2) + do i = 1, homogenization_maxNgrains !ToDo: this needs to be changed for varying Ngrains + do c = 1, size(material_phaseAt,1) + if (material_phaseAt(c,e) == instance) then + select_tensors(1:3,1:3,j) = dataset(1:3,1:3,c,i,e) + j = j + 1 + endif + enddo + enddo + enddo + + + end function select_tensors + + +!-------------------------------------------------------------------------------------------------- +!> @brief select rotations for output +!-------------------------------------------------------------------------------------------------- + function select_rotations(dataset,instance) + + use material, only: & + homogenization_maxNgrains, & + material_phaseAt + + integer, intent(in) :: instance + type(rotation), dimension(:,:,:), intent(in) :: dataset + type(rotation), allocatable, dimension(:) :: select_rotations + integer :: e,i,c,j + + allocate(select_rotations(count(material_phaseAt==instance)*homogenization_maxNgrains)) + + j=1 + do e = 1, size(material_phaseAt,2) + do i = 1, homogenization_maxNgrains !ToDo: this needs to be changed for varying Ngrains + do c = 1, size(material_phaseAt,1) + if (material_phaseAt(c,e) == instance) then + select_rotations(j) = dataset(c,i,e) + j = j + 1 + endif + enddo + enddo + enddo + + + end function select_rotations +#endif + + +end subroutine crystallite_results + + !-------------------------------------------------------------------------------------------------- !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction diff --git a/src/lattice.f90 b/src/lattice.f90 index 1b844c31f..d11932c29 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -507,10 +507,12 @@ module lattice public :: & lattice_init, & lattice_qDisorientation, & + LATTICE_iso_ID, & LATTICE_fcc_ID, & LATTICE_bcc_ID, & LATTICE_bct_ID, & LATTICE_hex_ID, & + LATTICE_ort_ID, & lattice_SchmidMatrix_slip, & lattice_SchmidMatrix_twin, & lattice_SchmidMatrix_trans, & @@ -581,18 +583,18 @@ subroutine lattice_init do p = 1, size(config_phase) tag = config_phase(p)%getString('lattice_structure') - select case(trim(tag)) - case('iso','isotropic') + select case(trim(tag(1:3))) + case('iso') lattice_structure(p) = LATTICE_iso_ID case('fcc') lattice_structure(p) = LATTICE_fcc_ID case('bcc') lattice_structure(p) = LATTICE_bcc_ID - case('hex','hexagonal') + case('hex') lattice_structure(p) = LATTICE_hex_ID case('bct') lattice_structure(p) = LATTICE_bct_ID - case('ort','orthorhombic') + case('ort') lattice_structure(p) = LATTICE_ort_ID end select diff --git a/src/material.f90 b/src/material.f90 index d35cfebd4..0b749c8ef 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -147,16 +147,14 @@ module material damage_initialPhi !< initial damage per each homogenization ! NEW MAPPINGS - integer(pInt), dimension(:), allocatable, public, protected :: & - material_homogenizationAt, & !< homogenization ID of each element (copy of mesh_homogenizationAt) - material_homogenizationMemberAt, & !< position of the element within its homogenization instance - material_aggregateAt, & !< aggregate ID of each element FUTURE USE FOR OUTPUT - material_aggregatMemberAt !< position of the element within its aggregate instance FUTURE USE FOR OUTPUT - integer(pInt), dimension(:,:), allocatable, public, protected :: & - material_phaseAt, & !< phase ID of each element - material_phaseMemberAt, & !< position of the element within its phase instance - material_crystalliteAt, & !< crystallite ID of each element CURRENTLY NOT PER CONSTITUTENT - material_crystalliteMemberAt !< position of the element within its crystallite instance CURRENTLY NOT PER CONSTITUTENT + integer, dimension(:), allocatable, public, protected :: & ! (elem) + material_homogenizationAt !< homogenization ID of each element (copy of mesh_homogenizationAt) + integer, dimension(:,:), allocatable, public, protected :: & ! (ip,elem) + material_homogenizationMemberAt !< position of the element within its homogenization instance + integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) + material_phaseAt !< phase ID of each element + integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,ip,elem) + material_phaseMemberAt !< position of the element within its phase instance ! END NEW MAPPINGS ! DEPRECATED: use material_phaseAt @@ -275,7 +273,10 @@ contains !> @details figures out if solverJobName.materialConfig is present, if not looks for !> material.config !-------------------------------------------------------------------------------------------------- -subroutine material_init() +subroutine material_init +#if defined(PETSc) || defined(DAMASK_HDF5) + use results +#endif use IO, only: & IO_error use debug, only: & @@ -382,21 +383,57 @@ subroutine material_init() endif debugOut call material_populateGrains + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! new mappings + allocate(material_homogenizationAt,source=theMesh%homogenizationAt) + allocate(material_homogenizationMemberAt(theMesh%elem%nIPs,theMesh%Nelems),source=0) + allocate(CounterHomogenization(size(config_homogenization)),source=0) + do e = 1, theMesh%Nelems + do i = 1, theMesh%elem%nIPs + CounterHomogenization(material_homogenizationAt(e)) = & + CounterHomogenization(material_homogenizationAt(e)) + 1 + material_homogenizationMemberAt(i,e) = CounterHomogenization(material_homogenizationAt(e)) + enddo + enddo + + + allocate(material_phaseAt(homogenization_maxNgrains,theMesh%Nelems), source=material_phase(:,1,:)) + allocate(material_phaseMemberAt(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0) + + allocate(CounterPhase(size(config_phase)),source=0) + do e = 1, theMesh%Nelems + do i = 1, theMesh%elem%nIPs + do c = 1, homogenization_maxNgrains + CounterPhase(material_phaseAt(c,e)) = & + CounterPhase(material_phaseAt(c,e)) + 1 + material_phaseMemberAt(c,i,e) = CounterPhase(material_phaseAt(c,e)) + enddo + enddo + enddo + +#if defined(PETSc) || defined(DAMASK_HDF5) + call results_openJobFile + call results_mapping_constituent(material_phaseAt,material_phaseMemberAt,phase_name) + call results_mapping_materialpoint(material_homogenizationAt,material_homogenizationMemberAt,homogenization_name) + call results_closeJobFile +#endif + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN DEPRECATED allocate(phaseAt ( homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt) allocate(phasememberAt ( homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt) allocate(mappingHomogenization (2, theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt) allocate(mappingHomogenizationConst( theMesh%elem%nIPs,theMesh%Nelems),source=1_pInt) -! END DEPRECATED - - allocate(material_homogenizationAt,source=theMesh%homogenizationAt) - allocate(material_AggregateAt, source=theMesh%homogenizationAt) - allocate(CounterPhase (size(config_phase)), source=0_pInt) - allocate(CounterHomogenization(size(config_homogenization)),source=0_pInt) + CounterHomogenization=0 + CounterPhase =0 + -! BEGIN DEPRECATED do e = 1_pInt,theMesh%Nelems myHomog = theMesh%homogenizationAt(e) do i = 1_pInt, theMesh%elem%nIPs diff --git a/src/numerics.f90 b/src/numerics.f90 index 955696219..f7c603c60 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -21,7 +21,7 @@ module numerics pert_method = 1_pInt, & !< method used in perturbation technique for tangent randomSeed = 0_pInt, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed worldrank = 0_pInt, & !< MPI worldrank (/=0 for MPI simulations only) - worldsize = 0_pInt, & !< MPI worldsize (/=0 for MPI simulations only) + worldsize = 1_pInt, & !< MPI worldsize (/=1 for MPI simulations only) numerics_integrator = 1_pInt !< method used for state integration Default 1: fix-point iteration integer(4), protected, public :: & DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index 88aa27432..19df4bdce 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -560,23 +560,41 @@ end function plastic_disloUCLA_postResults !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine plastic_disloUCLA_results(instance,group) -#if defined(PETSc) || defined(DAMASKHDF5) - use results +#if defined(PETSc) || defined(DAMASK_HDF5) + use results, only: & + results_writeDataset implicit none - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*), intent(in) :: group + integer :: o - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) outputsLoop: do o = 1,size(prm%outputID) select case(prm%outputID(o)) + case (rho_mob_ID) + call results_writeDataset(group,stt%rho_mob,'rho_mob',& + 'mobile dislocation density','1/m²') + case (rho_dip_ID) + call results_writeDataset(group,stt%rho_dip,'rho_dip',& + 'dislocation dipole density''1/m²') + case (dot_gamma_sl_ID) + call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& + 'plastic shear','1') + case (Lambda_sl_ID) + call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& + 'mean free path for slip','m') + case (thresholdstress_ID) + call results_writeDataset(group,dst%threshold_stress,'threshold_stress',& + 'threshold stress for slip','Pa') end select enddo outputsLoop end associate + #else - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*), intent(in) :: group #endif end subroutine plastic_disloUCLA_results diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90 index cb13265b4..8e52b3f41 100644 --- a/src/plastic_dislotwin.f90 +++ b/src/plastic_dislotwin.f90 @@ -148,7 +148,7 @@ module plastic_dislotwin type(tDislotwinState), allocatable, dimension(:), private :: & dotState, & state - type(tDislotwinMicrostructure), allocatable, dimension(:), private :: microstructure + type(tDislotwinMicrostructure), allocatable, dimension(:), private :: dependentState public :: & plastic_dislotwin_init, & @@ -241,14 +241,14 @@ subroutine plastic_dislotwin_init allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) - allocate(microstructure(Ninstance)) + allocate(dependentState(Ninstance)) do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_DISLOTWIN_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), & - dst => microstructure(phase_plasticityInstance(p)), & + dst => dependentState(phase_plasticityInstance(p)), & config => config_phase(p)) prm%aTol_rho = config%getFloat('atol_rho', defaultVal=0.0_pReal) @@ -801,7 +801,7 @@ subroutine plastic_dislotwin_dotState(Mp,T,instance,of) dot_gamma_tr associate(prm => param(instance), stt => state(instance), & - dot => dotstate(instance), dst => microstructure(instance)) + dot => dotstate(instance), dst => dependentState(instance)) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,of)) & @@ -897,7 +897,7 @@ subroutine plastic_dislotwin_dependentState(T,instance,of) associate(prm => param(instance),& stt => state(instance),& - dst => microstructure(instance)) + dst => dependentState(instance)) sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,of)) sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,of)) @@ -1002,7 +1002,7 @@ function plastic_dislotwin_postResults(Mp,T,instance,of) result(postResults) integer :: & o,c,j - associate(prm => param(instance), stt => state(instance), dst => microstructure(instance)) + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) c = 0 @@ -1063,20 +1063,53 @@ end function plastic_dislotwin_postResults !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine plastic_dislotwin_results(instance,group) -#if defined(PETSc) || defined(DAMASKHDF5) - use results +#if defined(PETSc) || defined(DAMASK_HDF5) + use results, only: & + results_writeDataset implicit none integer, intent(in) :: instance character(len=*) :: group integer :: o - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) outputsLoop: do o = 1,size(prm%outputID) select case(prm%outputID(o)) + + case (rho_mob_ID) + call results_writeDataset(group,stt%rho_mob,'rho_mob',& + 'mobile dislocation density','1/m²') + case (rho_dip_ID) + call results_writeDataset(group,stt%rho_dip,'rho_dip',& + 'dislocation dipole density''1/m²') + case (dot_gamma_sl_ID) + call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& + 'plastic shear','1') + case (Lambda_sl_ID) + call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& + 'mean free path for slip','m') + case (threshold_stress_slip_ID) + call results_writeDataset(group,dst%tau_pass,'tau_pass',& + 'passing stress for slip','Pa') + + case (f_tw_ID) + call results_writeDataset(group,stt%f_tw,'f_tw',& + 'twinned volume fraction','m³/m³') + case (Lambda_tw_ID) + call results_writeDataset(group,dst%Lambda_tw,'Lambda_tw',& + 'mean free path for twinning','m') + case (tau_hat_tw_ID) + call results_writeDataset(group,dst%tau_hat_tw,'tau_hat_tw',& + 'threshold stress for twinning','Pa') + + case (f_tr_ID) + call results_writeDataset(group,stt%f_tr,'f_tr',& + 'martensite volume fraction','m³/m³') + end select enddo outputsLoop end associate + #else integer, intent(in) :: instance character(len=*) :: group @@ -1130,7 +1163,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & tau_eff !< effective resolved stress integer :: i - associate(prm => param(instance), stt => state(instance), dst => microstructure(instance)) + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) do i = 1, prm%sum_N_sl tau(i) = math_mul33xx33(Mp,prm%P_sl(1:3,1:3,i)) @@ -1203,7 +1236,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& integer :: i,s1,s2 - associate(prm => param(instance), stt => state(instance), dst => microstructure(instance)) + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) do i = 1, prm%sum_N_tw tau(i) = math_mul33xx33(Mp,prm%P_tw(1:3,1:3,i)) @@ -1275,7 +1308,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& integer :: i,s1,s2 - associate(prm => param(instance), stt => state(instance), dst => microstructure(instance)) + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) do i = 1, prm%sum_N_tr tau(i) = math_mul33xx33(Mp,prm%P_tr(1:3,1:3,i)) diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index d7b071baa..facfa6d80 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -410,8 +410,7 @@ subroutine plastic_isotropic_dotState(Mp,instance,of) xi_inf_star = prm%xi_inf else xi_inf_star = prm%xi_inf & - + asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2) & - )**(1.0_pReal / prm%c_3) & + + asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) & / prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n) endif dot%xi(of) = dot_gamma & @@ -470,7 +469,7 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults) c = c + 1 case (dot_gamma_ID) postResults(c+1) = prm%dot_gamma_0 & - * (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n + * (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n c = c + 1 end select @@ -485,23 +484,28 @@ end function plastic_isotropic_postResults !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_results(instance,group) -#if defined(PETSc) || defined(DAMASKHDF5) - use results +#if defined(PETSc) || defined(DAMASK_HDF5) + use results, only: & + results_writeDataset implicit none - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*), intent(in) :: group + integer :: o associate(prm => param(instance), stt => state(instance)) outputsLoop: do o = 1,size(prm%outputID) select case(prm%outputID(o)) + case (xi_ID) + call results_writeDataset(group,stt%xi,'xi','resistance against plastic flow','Pa') end select enddo outputsLoop end associate + #else - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*), intent(in) :: group #endif end subroutine plastic_isotropic_results diff --git a/src/plastic_kinematichardening.f90 b/src/plastic_kinematichardening.f90 index ed21b09f7..04927c85b 100644 --- a/src/plastic_kinematichardening.f90 +++ b/src/plastic_kinematichardening.f90 @@ -551,8 +551,9 @@ end function plastic_kinehardening_postResults !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine plastic_kinehardening_results(instance,group) -#if defined(PETSc) || defined(DAMASKHDF5) - use results +#if defined(PETSc) || defined(DAMASK_HDF5) + use results, only: & + results_writeDataset implicit none integer, intent(in) :: instance @@ -562,6 +563,27 @@ subroutine plastic_kinehardening_results(instance,group) associate(prm => param(instance), stt => state(instance)) outputsLoop: do o = 1,size(prm%outputID) select case(prm%outputID(o)) + case (crss_ID) + call results_writeDataset(group,stt%crss,'xi_sl', & + 'resistance against plastic slip','Pa') + + case(crss_back_ID) + call results_writeDataset(group,stt%crss_back,'tau_back', & + 'back stress against plastic slip','Pa') + + case (sense_ID) + call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1') + + case (chi0_ID) + call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa') + + case (gamma0_ID) + call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1') + + case (accshear_ID) + call results_writeDataset(group,stt%accshear,'gamma_sl', & + 'plastic shear','1') + end select enddo outputsLoop end associate diff --git a/src/plastic_none.f90 b/src/plastic_none.f90 index f8a64b55b..b73bd20ab 100644 --- a/src/plastic_none.f90 +++ b/src/plastic_none.f90 @@ -5,13 +5,13 @@ !> @brief Dummy plasticity for purely elastic material !-------------------------------------------------------------------------------------------------- module plastic_none - - implicit none - private - - public :: & - plastic_none_init - + + implicit none + private + + public :: & + plastic_none_init + contains !-------------------------------------------------------------------------------------------------- @@ -19,41 +19,39 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_none_init - use debug, only: & - debug_level, & - debug_constitutive, & - debug_levelBasic - use material, only: & - phase_plasticity, & - material_allocatePlasticState, & - PLASTICITY_NONE_label, & - PLASTICITY_NONE_ID, & - material_phase, & - plasticState - - implicit none - integer :: & - Ninstance, & - p, & - NipcMyPhase - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' - - Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID) - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - do p = 1, size(phase_plasticity) - if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle - -!-------------------------------------------------------------------------------------------------- -! allocate state arrays - NipcMyPhase = count(material_phase == p) - call material_allocatePlasticState(p,NipcMyPhase,0,0,0, & - 0,0,0) - plasticState(p)%sizePostResults = 0 - - enddo + use debug, only: & + debug_level, & + debug_constitutive, & + debug_levelBasic + use material, only: & + phase_plasticity, & + material_allocatePlasticState, & + PLASTICITY_NONE_label, & + PLASTICITY_NONE_ID, & + material_phase, & + plasticState + + implicit none + integer :: & + Ninstance, & + p, & + NipcMyPhase + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' + + Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID) + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + do p = 1, size(phase_plasticity) + if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle + + NipcMyPhase = count(material_phase == p) + call material_allocatePlasticState(p,NipcMyPhase,0,0,0, & + 0,0,0) + plasticState(p)%sizePostResults = 0 + + enddo end subroutine plastic_none_init diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index f5f48ed11..a9ef98b06 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -2316,7 +2316,7 @@ outputsLoop: do o = 1,size(param(instance)%outputID) case (rho_dot_ann_ath_ID) postResults(cs+1:cs+ns) = results(instance)%rhoDotAthermalAnnihilation(1:ns,1,of) & - + results(instance)%rhoDotAthermalAnnihilation(1:ns,2,of) + + results(instance)%rhoDotAthermalAnnihilation(1:ns,2,of) cs = cs + ns case (rho_dot_ann_the_edge_ID) @@ -2402,8 +2402,9 @@ end function getRho !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_results(instance,group) -#if defined(PETSc) || defined(DAMASKHDF5) - use results +#if defined(PETSc) || defined(DAMASK_HDF5) + use results, only: & + results_writeDataset implicit none integer, intent(in) :: instance @@ -2413,6 +2414,39 @@ subroutine plastic_nonlocal_results(instance,group) associate(prm => param(instance), stt => state(instance)) outputsLoop: do o = 1,size(prm%outputID) select case(prm%outputID(o)) + case (rho_sgl_mob_edg_pos_ID) + call results_writeDataset(group,stt%rho_sgl_mob_edg_pos, 'rho_sgl_mob_edg_pos', & + 'positive mobile edge density','1/m²') + case (rho_sgl_imm_edg_pos_ID) + call results_writeDataset(group,stt%rho_sgl_imm_edg_pos, 'rho_sgl_imm_edg_pos',& + 'positive immobile edge density','1/m²') + case (rho_sgl_mob_edg_neg_ID) + call results_writeDataset(group,stt%rho_sgl_mob_edg_neg, 'rho_sgl_mob_edg_neg',& + 'negative mobile edge density','1/m²') + case (rho_sgl_imm_edg_neg_ID) + call results_writeDataset(group,stt%rho_sgl_imm_edg_neg, 'rho_sgl_imm_edg_neg',& + 'negative immobile edge density','1/m²') + case (rho_dip_edg_ID) + call results_writeDataset(group,stt%rho_dip_edg, 'rho_dip_edg',& + 'edge dipole density','1/m²') + case (rho_sgl_mob_scr_pos_ID) + call results_writeDataset(group,stt%rho_sgl_mob_scr_pos, 'rho_sgl_mob_scr_pos',& + 'positive mobile screw density','1/m²') + case (rho_sgl_imm_scr_pos_ID) + call results_writeDataset(group,stt%rho_sgl_imm_scr_pos, 'rho_sgl_imm_scr_pos',& + 'positive immobile screw density','1/m²') + case (rho_sgl_mob_scr_neg_ID) + call results_writeDataset(group,stt%rho_sgl_mob_scr_neg, 'rho_sgl_mob_scr_neg',& + 'negative mobile screw density','1/m²') + case (rho_sgl_imm_scr_neg_ID) + call results_writeDataset(group,stt%rho_sgl_imm_scr_neg, 'rho_sgl_imm_scr_neg',& + 'negative immobile screw density','1/m²') + case (rho_dip_scr_ID) + call results_writeDataset(group,stt%rho_dip_scr, 'rho_dip_scr',& + 'screw dipole density','1/m²') + case (rho_forest_ID) + call results_writeDataset(group,stt%rho_forest, 'rho_forest',& + 'forest density','1/m²') end select enddo outputsLoop end associate diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 4124856d1..272c4d631 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -563,28 +563,43 @@ end function plastic_phenopowerlaw_postResults !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine plastic_phenopowerlaw_results(instance,group) -#if defined(PETSc) || defined(DAMASKHDF5) - use results +#if defined(PETSc) || defined(DAMASK_HDF5) + use results, only: & + results_writeDataset - implicit none - integer, intent(in) :: instance - character(len=*) :: group - integer :: o + implicit none + integer, intent(in) :: instance + character(len=*), intent(in) :: group + + integer :: o - associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (resistance_slip_ID) - call results_writeVectorDataset(group,stt%xi_slip,'xi_slip','Pa') - case (accumulatedshear_slip_ID) - call results_writeVectorDataset(group,stt%gamma_slip,'gamma_slip','-') - end select - enddo outputsLoop - end associate + associate(prm => param(instance), stt => state(instance)) + outputsLoop: do o = 1,size(prm%outputID) + select case(prm%outputID(o)) + + case (resistance_slip_ID) + call results_writeDataset(group,stt%xi_slip, 'xi_sl', & + 'resistance against plastic slip','Pa') + case (accumulatedshear_slip_ID) + call results_writeDataset(group,stt%gamma_slip,'gamma_sl', & + 'plastic shear','1') + + case (resistance_twin_ID) + call results_writeDataset(group,stt%xi_twin, 'xi_tw', & + 'resistance against twinning','Pa') + case (accumulatedshear_twin_ID) + call results_writeDataset(group,stt%gamma_twin,'gamma_tw', & + 'twinning shear','1') + + end select + enddo outputsLoop + end associate + #else - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*), intent(in) :: group #endif + end subroutine plastic_phenopowerlaw_results diff --git a/src/results.f90 b/src/results.f90 index f70e124f8..4ed5cc751 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -5,9 +5,6 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !-------------------------------------------------------------------------------------------------- module results - use prec - use IO - use HDF5 use HDF5_utilities #ifdef PETSc use PETSC @@ -18,6 +15,18 @@ module results integer(HID_T), public, protected :: tempCoordinates, tempResults integer(HID_T), private :: resultsFile, currentIncID, plist_id + interface results_writeDataset + + module procedure results_writeTensorDataset_real + module procedure results_writeVectorDataset_real + module procedure results_writeScalarDataset_real + + module procedure results_writeTensorDataset_int + module procedure results_writeVectorDataset_int + + module procedure results_writeScalarDataset_rotation + + end interface results_writeDataset public :: & results_init, & @@ -26,23 +35,35 @@ module results results_addIncrement, & results_addGroup, & results_openGroup, & - results_writeVectorDataset, & + results_writeDataset, & results_setLink, & - results_removeLink + results_addAttribute, & + results_removeLink, & + results_mapping_constituent, & + results_mapping_materialpoint contains subroutine results_init use DAMASK_interface, only: & getSolverJobName - implicit none + character(len=pStringLen) :: commandLine write(6,'(/,a)') ' <<<+- results init -+>>>' write(6,'(/,a)') ' Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017' write(6,'(a)') ' https://doi.org/10.1007/s40192-018-0118-7' - call HDF5_closeFile(HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.)) + resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.) + call HDF5_addAttribute(resultsFile,'DADF5-version',0.1) + call HDF5_addAttribute(resultsFile,'DADF5-major',0) + call HDF5_addAttribute(resultsFile,'DADF5-minor',1) + call HDF5_addAttribute(resultsFile,'DAMASK',DAMASKVERSION) + call get_command(commandLine) + call HDF5_addAttribute(resultsFile,'call',trim(commandLine)) + call HDF5_closeGroup(results_addGroup('mapping')) + call HDF5_closeGroup(results_addGroup('mapping/cellResults')) + call HDF5_closeFile(resultsFile) end subroutine results_init @@ -50,18 +71,12 @@ end subroutine results_init !-------------------------------------------------------------------------------------------------- !> @brief opens the results file to append data !-------------------------------------------------------------------------------------------------- -subroutine results_openJobFile() - use DAMASK_interface, only: & - getSolverJobName +subroutine results_openJobFile + use DAMASK_interface, only: & + getSolverJobName - implicit none - character(len=pStringLen) :: commandLine - resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','a',.true.) - call HDF5_addAttribute(resultsFile,'DADF5',0.1_pReal) - call HDF5_addAttribute(resultsFile,'DAMASK',DAMASKVERSION) - call get_command(commandLine) - call HDF5_addAttribute(resultsFile,'call',trim(commandLine)) + resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','a',.true.) end subroutine results_openJobFile @@ -69,10 +84,9 @@ end subroutine results_openJobFile !-------------------------------------------------------------------------------------------------- !> @brief closes the results file !-------------------------------------------------------------------------------------------------- -subroutine results_closeJobFile() - implicit none +subroutine results_closeJobFile - call HDF5_closeFile(resultsFile) + call HDF5_closeFile(resultsFile) end subroutine results_closeJobFile @@ -82,15 +96,14 @@ end subroutine results_closeJobFile !-------------------------------------------------------------------------------------------------- subroutine results_addIncrement(inc,time) - implicit none - integer(pInt), intent(in) :: inc - real(pReal), intent(in) :: time - character(len=pStringLen) :: incChar + integer(pInt), intent(in) :: inc + real(pReal), intent(in) :: time + character(len=pStringLen) :: incChar - write(incChar,*) inc - call HDF5_closeGroup(results_addGroup(trim('inc'//trim(adjustl(incChar))))) - call results_setLink(trim('inc'//trim(adjustl(incChar))),'current') - call HDF5_addAttribute(resultsFile,'time/s',time,trim('inc'//trim(adjustl(incChar)))) + write(incChar,'(i5.5)') inc ! allow up to 99999 increments + call HDF5_closeGroup(results_addGroup(trim('inc'//trim(adjustl(incChar))))) + call results_setLink(trim('inc'//trim(adjustl(incChar))),'current') + call HDF5_addAttribute(resultsFile,'time/s',time,trim('inc'//trim(adjustl(incChar)))) end subroutine results_addIncrement @@ -99,10 +112,9 @@ end subroutine results_addIncrement !-------------------------------------------------------------------------------------------------- integer(HID_T) function results_openGroup(groupName) - implicit none - character(len=*), intent(in) :: groupName - - results_openGroup = HDF5_openGroup(resultsFile,groupName) + character(len=*), intent(in) :: groupName + + results_openGroup = HDF5_openGroup(resultsFile,groupName) end function results_openGroup @@ -112,10 +124,9 @@ end function results_openGroup !-------------------------------------------------------------------------------------------------- integer(HID_T) function results_addGroup(groupName) - implicit none - character(len=*), intent(in) :: groupName - - results_addGroup = HDF5_addGroup(resultsFile,groupName) + character(len=*), intent(in) :: groupName + + results_addGroup = HDF5_addGroup(resultsFile,groupName) end function results_addGroup @@ -124,852 +135,733 @@ end function results_addGroup !> @brief set link to object in results file !-------------------------------------------------------------------------------------------------- subroutine results_setLink(path,link) - use hdf5_utilities, only: & - HDF5_setLink - implicit none - character(len=*), intent(in) :: path, link + character(len=*), intent(in) :: path, link - call HDF5_setLink(resultsFile,path,link) + call HDF5_setLink(resultsFile,path,link) end subroutine results_setLink +!-------------------------------------------------------------------------------------------------- +!> @brief adds an attribute to an object +!-------------------------------------------------------------------------------------------------- +subroutine results_addAttribute(attrLabel,attrValue,path) + + character(len=*), intent(in) :: attrLabel, attrValue, path + + call HDF5_addAttribute_str(resultsFile,attrLabel, attrValue, path) + +end subroutine results_addAttribute + + !-------------------------------------------------------------------------------------------------- !> @brief remove link to an object !-------------------------------------------------------------------------------------------------- subroutine results_removeLink(link) - use hdf5 - implicit none - character(len=*), intent(in) :: link - integer :: hdferr + character(len=*), intent(in) :: link + integer :: hdferr - call h5ldelete_f(resultsFile,link, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'results_removeLink: h5ldelete_soft_f ('//trim(link)//')') + call h5ldelete_f(resultsFile,link, hdferr) + if (hdferr < 0) call IO_error(1,ext_msg = 'results_removeLink: h5ldelete_soft_f ('//trim(link)//')') end subroutine results_removeLink +!-------------------------------------------------------------------------------------------------- +!> @brief stores a scalar dataset in a group +!-------------------------------------------------------------------------------------------------- +subroutine results_writeScalarDataset_real(group,dataset,label,description,SIunit) + + character(len=*), intent(in) :: label,group,description + character(len=*), intent(in), optional :: SIunit + real(pReal), intent(inout), dimension(:) :: dataset + + integer(HID_T) :: groupHandle + + groupHandle = results_openGroup(group) + +#ifdef PETSc + call HDF5_write(groupHandle,dataset,label,.true.) +#endif + + if (HDF5_objectExists(groupHandle,label)) & + call HDF5_addAttribute(groupHandle,'Description',description,label) + if (HDF5_objectExists(groupHandle,label) .and. present(SIunit)) & + call HDF5_addAttribute(groupHandle,'Unit',SIunit,label) + call HDF5_closeGroup(groupHandle) + +end subroutine results_writeScalarDataset_real + !-------------------------------------------------------------------------------------------------- !> @brief stores a vector dataset in a group !-------------------------------------------------------------------------------------------------- -subroutine results_writeVectorDataset(group,dataset,label,SIunit) +subroutine results_writeVectorDataset_real(group,dataset,label,description,SIunit) - implicit none - character(len=*), intent(in) :: SIunit,label,group - real(pReal), intent(inout), dimension(:,:) :: dataset - integer(HID_T) :: groupHandle + character(len=*), intent(in) :: label,group,description + character(len=*), intent(in), optional :: SIunit + real(pReal), intent(inout), dimension(:,:) :: dataset + + integer(HID_T) :: groupHandle - groupHandle = results_openGroup(group) - call HDF5_write(groupHandle,dataset,label) - if (HDF5_objectExists(groupHandle,label)) & - call HDF5_addAttribute(groupHandle,'Unit',SIunit,label) - call HDF5_closeGroup(groupHandle) + groupHandle = results_openGroup(group) + +#ifdef PETSc + call HDF5_write(groupHandle,dataset,label,.true.) +#endif + + if (HDF5_objectExists(groupHandle,label)) & + call HDF5_addAttribute(groupHandle,'Description',description,label) + if (HDF5_objectExists(groupHandle,label) .and. present(SIunit)) & + call HDF5_addAttribute(groupHandle,'Unit',SIunit,label) + call HDF5_closeGroup(groupHandle) -end subroutine results_writeVectorDataset +end subroutine results_writeVectorDataset_real + + +!-------------------------------------------------------------------------------------------------- +!> @brief stores a tensor dataset in a group +!-------------------------------------------------------------------------------------------------- +subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit) + + character(len=*), intent(in) :: label,group,description + character(len=*), intent(in), optional :: SIunit + real(pReal), intent(inout), dimension(:,:,:) :: dataset + + integer(HID_T) :: groupHandle + + groupHandle = results_openGroup(group) + +#ifdef PETSc + call HDF5_write(groupHandle,dataset,label,.true.) +#endif + + if (HDF5_objectExists(groupHandle,label)) & + call HDF5_addAttribute(groupHandle,'Description',description,label) + if (HDF5_objectExists(groupHandle,label) .and. present(SIunit)) & + call HDF5_addAttribute(groupHandle,'Unit',SIunit,label) + call HDF5_closeGroup(groupHandle) + +end subroutine results_writeTensorDataset_real + + +!-------------------------------------------------------------------------------------------------- +!> @brief stores a vector dataset in a group +!-------------------------------------------------------------------------------------------------- +subroutine results_writeVectorDataset_int(group,dataset,label,description,SIunit) + + character(len=*), intent(in) :: label,group,description + character(len=*), intent(in), optional :: SIunit + integer, intent(inout), dimension(:,:) :: dataset + + integer(HID_T) :: groupHandle + + groupHandle = results_openGroup(group) + +#ifdef PETSc + call HDF5_write(groupHandle,dataset,label,.true.) +#endif + + if (HDF5_objectExists(groupHandle,label)) & + call HDF5_addAttribute(groupHandle,'Description',description,label) + if (HDF5_objectExists(groupHandle,label) .and. present(SIunit)) & + call HDF5_addAttribute(groupHandle,'Unit',SIunit,label) + call HDF5_closeGroup(groupHandle) + +end subroutine results_writeVectorDataset_int + + +!-------------------------------------------------------------------------------------------------- +!> @brief stores a vector dataset in a group +!-------------------------------------------------------------------------------------------------- +subroutine results_writeTensorDataset_int(group,dataset,label,description,SIunit) + + character(len=*), intent(in) :: label,group,description + character(len=*), intent(in), optional :: SIunit + integer, intent(inout), dimension(:,:,:) :: dataset + + integer(HID_T) :: groupHandle + + groupHandle = results_openGroup(group) + +#ifdef PETSc + call HDF5_write(groupHandle,dataset,label,.true.) +#endif + + if (HDF5_objectExists(groupHandle,label)) & + call HDF5_addAttribute(groupHandle,'Description',description,label) + if (HDF5_objectExists(groupHandle,label) .and. present(SIunit)) & + call HDF5_addAttribute(groupHandle,'Unit',SIunit,label) + call HDF5_closeGroup(groupHandle) + +end subroutine results_writeTensorDataset_int + + +!-------------------------------------------------------------------------------------------------- +!> @brief stores a vector dataset in a group +!-------------------------------------------------------------------------------------------------- +subroutine results_writeScalarDataset_rotation(group,dataset,label,description,lattice_structure) + use rotations, only: & + rotation + + character(len=*), intent(in) :: label,group,description + character(len=*), intent(in), optional :: lattice_structure + type(rotation), intent(inout), dimension(:) :: dataset + + integer(HID_T) :: groupHandle + + groupHandle = results_openGroup(group) + +#ifdef PETSc + call HDF5_write(groupHandle,dataset,label,.true.) +#endif + + if (HDF5_objectExists(groupHandle,label)) & + call HDF5_addAttribute(groupHandle,'Description',description,label) + if (HDF5_objectExists(groupHandle,label) .and. present(lattice_structure)) & + call HDF5_addAttribute(groupHandle,'Lattice',lattice_structure,label) + call HDF5_closeGroup(groupHandle) + +end subroutine results_writeScalarDataset_rotation !-------------------------------------------------------------------------------------------------- !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- -subroutine HDF5_mappingPhase(mapping,mapping2,Nconstituents,material_phase,phase_name,dataspace_size,mpiOffset,mpiOffset_phase) - use hdf5 +subroutine results_mapping_constituent(phaseAt,memberAt,label) + use numerics, only: & + worldrank, & + worldsize + + integer, dimension(:,:), intent(in) :: phaseAt !< phase section at (constituent,element) + integer, dimension(:,:,:), intent(in) :: memberAt !< phase member at (constituent,IP,element) + character(len=64), dimension(:), intent(in) :: label !< label of each phase section + + integer, dimension(size(memberAt,1),size(memberAt,2),size(memberAt,3)) :: & + phaseAt_perIP, & + memberAt_total + integer, dimension(size(label),0:worldsize-1) :: memberOffset !< offset in member counting per process + integer, dimension(0:worldsize-1) :: writeSize !< amount of data written per process + integer(HSIZE_T), dimension(2) :: & + myShape, & !< shape of the dataset (this process) + myOffset, & + totalShape !< shape of the dataset (all processes) + + integer(HID_T) :: & + loc_id, & !< identifier of group in file + dtype_id, & !< identifier of compound data type + name_id, & !< identifier of name (string) in compound data type + position_id, & !< identifier of position/index (integer) in compound data type + dset_id, & + memspace_id, & + filespace_id, & + plist_id, & + dt_id - implicit none - integer(pInt), intent(in) :: Nconstituents, dataspace_size, mpiOffset - integer(pInt), intent(in), dimension(:) :: mapping, mapping2 - character(len=*), intent(in), dimension(:) :: phase_name - integer(pInt), intent(in), dimension(:) :: mpiOffset_phase - integer(pInt), intent(in), dimension(:,:,:) :: material_phase - - character(len=len(phase_name(1))), dimension(:), allocatable :: namesNA - character(len=len(phase_name(1))) :: a - character(len=*), parameter :: n = "NULL" - - integer(pInt) :: hdferr, NmatPoints, i, j, k - integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id, name_id, position_id, plist_id, memspace - integer(HID_T) :: dt5_id ! Memory datatype identifier - integer(SIZE_T) :: typesize, type_sizec, type_sizei, type_size - - integer(HSIZE_T), dimension(2) :: counter - integer(HSSIZE_T), dimension(2) :: fileOffset - integer(pInt), dimension(:,:), allocatable :: arrOffset - - a = n - allocate(namesNA(0:size(phase_name)),source=[a,phase_name]) - NmatPoints = size(mapping,1)/Nconstituents - mapping_ID = results_openGroup("current/mapGeometry") - - allocate(arrOffset(Nconstituents,NmatPoints)) - do i=1_pInt, NmatPoints - do k=1_pInt, Nconstituents - do j=1_pInt, size(phase_name) - if(material_phase(k,1,i) == j) & - arrOffset(k,i) = mpiOffset_phase(j) - enddo - enddo - enddo + + integer(SIZE_T) :: type_size_string, type_size_int + integer :: ierr, i + +!--------------------------------------------------------------------------------------------------- +! compound type: name of phase section + position/index within results array + call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, ierr) + call h5tset_size_f(dt_id, int(len(label(1)),SIZE_T), ierr) + call h5tget_size_f(dt_id, type_size_string, ierr) + + call h5tget_size_f(H5T_NATIVE_INTEGER, type_size_int, ierr) + + call h5tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, ierr) + call h5tinsert_f(dtype_id, "Name", 0_SIZE_T, dt_id,ierr) + call h5tinsert_f(dtype_id, "Position", type_size_string, H5T_NATIVE_INTEGER, ierr) + +!-------------------------------------------------------------------------------------------------- +! create memory types for each component of the compound type + call h5tcreate_f(H5T_COMPOUND_F, type_size_string, name_id, ierr) + call h5tinsert_f(name_id, "Name", 0_SIZE_T, dt_id, ierr) + + call h5tcreate_f(H5T_COMPOUND_F, type_size_int, position_id, ierr) + call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_NATIVE_INTEGER, ierr) + + call h5tclose_f(dt_id, ierr) !-------------------------------------------------------------------------------------------------- -! create dataspace - call h5screate_simple_f(2, int([Nconstituents,dataspace_size],HSIZE_T), space_id, hdferr, & - int([Nconstituents,dataspace_size],HSIZE_T)) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeMapping') +! prepare MPI communication (transparent for non-MPI runs) + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr) + memberOffset = 0 + do i=1, size(label) + memberOffset(i,worldrank) = count(phaseAt == i)*size(memberAt,2) ! number of points/instance of this process + enddo + writeSize = 0 + writeSize(worldrank) = size(memberAt(1,:,:)) ! total number of points by this process !-------------------------------------------------------------------------------------------------- -! compound type - ! First calculate total size by calculating sizes of each member - ! - CALL h5tcopy_f(H5T_NATIVE_CHARACTER, dt5_id, hdferr) - typesize = len(phase_name(1)) - CALL h5tset_size_f(dt5_id, typesize, hdferr) - CALL h5tget_size_f(dt5_id, type_sizec, hdferr) - CALL h5tget_size_f(H5T_STD_I32LE,type_sizei, hdferr) - type_size = type_sizec + type_sizei - call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeMapping: h5tcreate_f dtype_id') - - call h5tinsert_f(dtype_id, "Name", 0_SIZE_T, dt5_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tinsert_f 0') - call h5tinsert_f(dtype_id, "Position", type_sizec, H5T_STD_I32LE, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tinsert_f 2') - -!-------------------------------------------------------------------------------------------------- -! create Dataset - call h5dcreate_f(mapping_id, 'constitutive', dtype_id, space_id, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase') - -!-------------------------------------------------------------------------------------------------- -! Create memory types (one compound datatype for each member) - call h5tcreate_f(H5T_COMPOUND_F, int(type_sizec,SIZE_T), name_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tcreate_f instance_id') - call h5tinsert_f(name_id, "Name", 0_SIZE_T, dt5_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tinsert_f instance_id') - - call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tcreate_f position_id') - call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tinsert_f position_id') - -!-------------------------------------------------------------------------------------------------- -! Define and select hyperslabs - counter(1) = Nconstituents ! how big i am - counter(2) = NmatPoints - fileOffset(1) = 0 ! where i start to write my data - fileOffset(2) = mpiOffset - - call h5screate_simple_f(2, counter, memspace, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5screate_simple_f') - call h5dget_space_f(dset_id, space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5dget_space_f') - call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5sselect_hyperslab_f') - -!-------------------------------------------------------------------------------------------------- - ! Create property list for collective dataset write +! MPI settings and communication #ifdef PETSc - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5pcreate_f') - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5pset_dxpl_mpio_f') + call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, ierr) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5pset_dxpl_mpio_f') + + call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process + if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_constituent: MPI_allreduce/writeSize') + + call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process + if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_constituent: MPI_allreduce/memberOffset') #endif + myShape = int([size(phaseAt,1),writeSize(worldrank)], HSIZE_T) + myOffset = int([0,sum(writeSize(0:worldrank-1))], HSIZE_T) + totalShape = int([size(phaseAt,1),sum(writeSize)], HSIZE_T) + !-------------------------------------------------------------------------------------------------- -! write data by fields in the datatype. Fields order is not important. - call h5dwrite_f(dset_id, name_id, reshape(namesNA(mapping),[Nconstituents,NmatPoints]), & - int([Nconstituents, dataspace_size],HSIZE_T), hdferr, & - file_space_id = space_id, mem_space_id = memspace, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5dwrite_f position_id') +! create dataspace in memory (local shape = hyperslab) and in file (global shape) + call h5screate_simple_f(2,myShape,memspace_id,ierr,myShape) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5screate_simple_f/memspace_id') + + call h5screate_simple_f(2,totalShape,filespace_id,ierr,totalShape) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5screate_simple_f/filespace_id') + + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, ierr) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5sselect_hyperslab_f') - call h5dwrite_f(dset_id, position_id, reshape(mapping2-1_pInt,[Nconstituents,NmatPoints])+arrOffset, & - int([Nconstituents, dataspace_size],HSIZE_T), hdferr, & - file_space_id = space_id, mem_space_id = memspace, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5dwrite_f instance_id') +!--------------------------------------------------------------------------------------------------- +! expand phaseAt to consider IPs (is not stored per IP) + do i = 1, size(phaseAt_perIP,2) + phaseAt_perIP(:,i,:) = phaseAt + enddo + +!--------------------------------------------------------------------------------------------------- +! renumber member from my process to all processes + do i = 1, size(label) + where(phaseAt_perIP == i) memberAt_total = memberAt + sum(memberOffset(i,0:worldrank-1)) + enddo !-------------------------------------------------------------------------------------------------- -! close types, dataspaces - call h5tclose_f(dtype_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tclose_f dtype_id') - call h5tclose_f(position_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tclose_f position_id') - call h5tclose_f(name_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tclose_f name_id ') - call h5tclose_f(dt5_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tclose_f dt5_id') - call h5dclose_f(dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5dclose_f') - call h5sclose_f(space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5sclose_f space_id') - call h5sclose_f(memspace, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5sclose_f memspace') - call h5pclose_f(plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5pclose_f') - call HDF5_closeGroup(mapping_ID) - -end subroutine HDF5_mappingPhase +! write the components of the compound type individually + call h5pset_preserve_f(plist_id, .TRUE., ierr) + + loc_id = results_openGroup('/mapping/cellResults') + call h5dcreate_f(loc_id, 'constituent', dtype_id, filespace_id, dset_id, ierr) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5dcreate_f') + + call h5dwrite_f(dset_id, name_id, reshape(label(pack(phaseAt_perIP,.true.)),myShape), & + myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5dwrite_f/name_id') + call h5dwrite_f(dset_id, position_id, reshape(pack(memberAt_total,.true.),myShape), & + myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5dwrite_f/position_id') !-------------------------------------------------------------------------------------------------- -!> @brief adds the backward mapping from spatial position and constituent ID to results -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_backwardMappingPhase(material_phase,phasememberat,phase_name,dataspace_size,mpiOffset,mpiOffset_phase) - use hdf5 +! close all + call HDF5_closeGroup(loc_id) + call h5pclose_f(plist_id, ierr) + call h5sclose_f(filespace_id, ierr) + call h5sclose_f(memspace_id, ierr) + call h5dclose_f(dset_id, ierr) + call h5tclose_f(dtype_id, ierr) + call h5tclose_f(name_id, ierr) + call h5tclose_f(position_id, ierr) - implicit none - integer(pInt), intent(in), dimension(:,:,:) :: material_phase, phasememberat - character(len=*), intent(in), dimension(:) :: phase_name - integer(pInt), intent(in), dimension(:) :: dataspace_size, mpiOffset_phase - integer(pInt), intent(in) :: mpiOffset +end subroutine results_mapping_constituent - integer(pInt) :: hdferr, NmatPoints, Nconstituents, i, j - integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id, position_id, plist_id, memspace - integer(SIZE_T) :: type_size - - integer(pInt), dimension(:,:), allocatable :: arr - - integer(HSIZE_T), dimension(1) :: counter - integer(HSSIZE_T), dimension(1) :: fileOffset - - character(len=64) :: phaseID - - Nconstituents = size(phasememberat,1) - NmatPoints = count(material_phase /=0_pInt)/Nconstituents - - allocate(arr(2,NmatPoints*Nconstituents)) - - do i=1_pInt, NmatPoints - do j=Nconstituents-1_pInt, 0_pInt, -1_pInt - arr(1,Nconstituents*i-j) = i-1_pInt - enddo - enddo - arr(2,:) = pack(material_phase,material_phase/=0_pInt) - - do i=1_pInt, size(phase_name) - write(phaseID, '(i0)') i - mapping_ID = results_openGroup('/current/constitutive/'//trim(phaseID)//'_'//phase_name(i)) - NmatPoints = count(material_phase == i) - -!-------------------------------------------------------------------------------------------------- - ! create dataspace - call h5screate_simple_f(1, int([dataspace_size(i)],HSIZE_T), space_id, hdferr, & - int([dataspace_size(i)],HSIZE_T)) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeBackwardMapping') - -!-------------------------------------------------------------------------------------------------- - ! compound type - call h5tget_size_f(H5T_STD_I32LE, type_size, hdferr) - call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeBackwardMapping: h5tcreate_f dtype_id') - - call h5tinsert_f(dtype_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5tinsert_f 0') - -!-------------------------------------------------------------------------------------------------- - ! create Dataset - call h5dcreate_f(mapping_id, 'mapGeometry', dtype_id, space_id, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase') - -!-------------------------------------------------------------------------------------------------- - ! Create memory types (one compound datatype for each member) - call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5tcreate_f position_id') - call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5tinsert_f position_id') - -!-------------------------------------------------------------------------------------------------- - ! Define and select hyperslabs - counter = NmatPoints ! how big i am - fileOffset = mpiOffset_phase(i) ! where i start to write my data - - call h5screate_simple_f(1, counter, memspace, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5screate_simple_f') - call h5dget_space_f(dset_id, space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5dget_space_f') - call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5sselect_hyperslab_f') - -!-------------------------------------------------------------------------------------------------- - ! Create property list for collective dataset write -#ifdef PETSc - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5pcreate_f') - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5pset_dxpl_mpio_f') -#endif - -!-------------------------------------------------------------------------------------------------- - ! write data by fields in the datatype. Fields order is not important. - call h5dwrite_f(dset_id, position_id, pack(arr(1,:),arr(2,:)==i)+mpiOffset, int([dataspace_size(i)],HSIZE_T),& - hdferr, file_space_id = space_id, mem_space_id = memspace, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5dwrite_f instance_id') - -!-------------------------------------------------------------------------------------------------- - !close types, dataspaces - call h5tclose_f(dtype_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5tclose_f dtype_id') - call h5tclose_f(position_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5tclose_f position_id') - call h5dclose_f(dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5dclose_f') - call h5sclose_f(space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5sclose_f space_id') - call h5sclose_f(memspace, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5sclose_f memspace') - call h5pclose_f(plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5pclose_f') - call HDF5_closeGroup(mapping_ID) - - enddo - -end subroutine HDF5_backwardMappingPhase !-------------------------------------------------------------------------------------------------- !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- -subroutine HDF5_mappingHomog(material_homog,homogmemberat,homogenization_name,dataspace_size,mpiOffset,mpiOffset_homog) - use hdf5 +subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label) + use numerics, only: & + worldrank, & + worldsize + + integer, dimension(:), intent(in) :: homogenizationAt !< homogenization section at (element) + integer, dimension(:,:), intent(in) :: memberAt !< homogenization member at (IP,element) + character(len=64), dimension(:), intent(in) :: label !< label of each homogenization section + + integer, dimension(size(memberAt,1),size(memberAt,2)) :: & + homogenizationAt_perIP, & + memberAt_total + integer, dimension(size(label),0:worldsize-1) :: memberOffset !< offset in member counting per process + integer, dimension(0:worldsize-1) :: writeSize !< amount of data written per process + integer(HSIZE_T), dimension(1) :: & + myShape, & !< shape of the dataset (this process) + myOffset, & + totalShape !< shape of the dataset (all processes) + + integer(HID_T) :: & + loc_id, & !< identifier of group in file + dtype_id, & !< identifier of compound data type + name_id, & !< identifier of name (string) in compound data type + position_id, & !< identifier of position/index (integer) in compound data type + dset_id, & + memspace_id, & + filespace_id, & + plist_id, & + dt_id - implicit none - integer(pInt), intent(in), dimension(:,:) :: material_homog, homogmemberat - character(len=*), intent(in), dimension(:) :: homogenization_name - integer(pInt), intent(in), dimension(:) :: mpiOffset_homog - integer(pInt), intent(in) :: dataspace_size, mpiOffset - - integer(pInt) :: hdferr, NmatPoints, i, j - integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id, name_id, position_id, plist_id, memspace - - integer(HID_T) :: dt5_id ! Memory datatype identifier - integer(SIZE_T) :: typesize, type_sizec, type_sizei, type_size - - integer(HSIZE_T), dimension(1) :: counter - integer(HSSIZE_T), dimension(1) :: fileOffset - integer(pInt), dimension(:), allocatable :: arrOffset - - NmatPoints = count(material_homog /=0_pInt) - mapping_ID = results_openGroup("current/mapGeometry") - - allocate(arrOffset(NmatPoints)) - do i=1_pInt, NmatPoints - do j=1_pInt, size(homogenization_name) - if(material_homog(1,i) == j) & - arrOffset(i) = mpiOffset_homog(j) - enddo - enddo + + integer(SIZE_T) :: type_size_string, type_size_int + integer :: ierr, i + +!--------------------------------------------------------------------------------------------------- +! compound type: name of phase section + position/index within results array + call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, ierr) + call h5tset_size_f(dt_id, int(len(label(1)),SIZE_T), ierr) + call h5tget_size_f(dt_id, type_size_string, ierr) + + call h5tget_size_f(H5T_NATIVE_INTEGER, type_size_int, ierr) + + call h5tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, ierr) + call h5tinsert_f(dtype_id, "Name", 0_SIZE_T, dt_id,ierr) + call h5tinsert_f(dtype_id, "Position", type_size_string, H5T_NATIVE_INTEGER, ierr) + +!-------------------------------------------------------------------------------------------------- +! create memory types for each component of the compound type + call h5tcreate_f(H5T_COMPOUND_F, type_size_string, name_id, ierr) + call h5tinsert_f(name_id, "Name", 0_SIZE_T, dt_id, ierr) + + call h5tcreate_f(H5T_COMPOUND_F, type_size_int, position_id, ierr) + call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_NATIVE_INTEGER, ierr) + + call h5tclose_f(dt_id, ierr) !-------------------------------------------------------------------------------------------------- -! create dataspace - call h5screate_simple_f(1, int([dataspace_size],HSIZE_T), space_id, hdferr, & - int([dataspace_size],HSIZE_T)) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeMapping') +! prepare MPI communication (transparent for non-MPI runs) + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr) + memberOffset = 0 + do i=1, size(label) + memberOffset(i,worldrank) = count(homogenizationAt == i)*size(memberAt,1) ! number of points/instance of this process + enddo + writeSize = 0 + writeSize(worldrank) = size(memberAt) ! total number of points by this process !-------------------------------------------------------------------------------------------------- -! compound type - ! First calculate total size by calculating sizes of each member - ! - CALL h5tcopy_f(H5T_NATIVE_CHARACTER, dt5_id, hdferr) - typesize = len(homogenization_name(1)) - CALL h5tset_size_f(dt5_id, typesize, hdferr) - CALL h5tget_size_f(dt5_id, type_sizec, hdferr) - CALL h5tget_size_f(H5T_STD_I32LE,type_sizei, hdferr) - type_size = type_sizec + type_sizei - call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeMapping: h5tcreate_f dtype_id') - - call h5tinsert_f(dtype_id, "Name", 0_SIZE_T, dt5_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tinsert_f 0') - call h5tinsert_f(dtype_id, "Position", type_sizec, H5T_STD_I32LE, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tinsert_f 2') - -!-------------------------------------------------------------------------------------------------- -! create Dataset - call h5dcreate_f(mapping_id, 'homogenization', dtype_id, space_id, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog') - -!-------------------------------------------------------------------------------------------------- -! Create memory types (one compound datatype for each member) - call h5tcreate_f(H5T_COMPOUND_F, int(type_sizec,SIZE_T), name_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tcreate_f instance_id') - call h5tinsert_f(name_id, "Name", 0_SIZE_T, dt5_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tinsert_f instance_id') - - call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tcreate_f position_id') - call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tinsert_f position_id') - -!-------------------------------------------------------------------------------------------------- -! Define and select hyperslabs - counter = NmatPoints ! how big i am - fileOffset = mpiOffset ! where i start to write my data - - call h5screate_simple_f(1, counter, memspace, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5screate_simple_f') - call h5dget_space_f(dset_id, space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5dget_space_f') - call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5sselect_hyperslab_f') - -!-------------------------------------------------------------------------------------------------- -! Create property list for collective dataset write +! MPI settings and communication #ifdef PETSc - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5pcreate_f') - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5pset_dxpl_mpio_f') + call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, ierr) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5pset_dxpl_mpio_f') + + call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process + if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_materialpoint: MPI_allreduce/writeSize') + + call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process + if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_materialpoint: MPI_allreduce/memberOffset') #endif + myShape = int([writeSize(worldrank)], HSIZE_T) + myOffset = int([sum(writeSize(0:worldrank-1))], HSIZE_T) + totalShape = int([sum(writeSize)], HSIZE_T) + !-------------------------------------------------------------------------------------------------- -! write data by fields in the datatype. Fields order is not important. - call h5dwrite_f(dset_id, name_id, homogenization_name(pack(material_homog,material_homog/=0_pInt)), & - int([dataspace_size],HSIZE_T), hdferr, file_space_id = space_id, & - mem_space_id = memspace, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5dwrite_f position_id') +! create dataspace in memory (local shape = hyperslab) and in file (global shape) + call h5screate_simple_f(1,myShape,memspace_id,ierr,myShape) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5screate_simple_f/memspace_id') + + call h5screate_simple_f(1,totalShape,filespace_id,ierr,totalShape) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5screate_simple_f/filespace_id') + + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, ierr) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5sselect_hyperslab_f') - call h5dwrite_f(dset_id, position_id, pack(homogmemberat-1_pInt,homogmemberat/=0_pInt) + arrOffset, & - int([dataspace_size],HSIZE_T), hdferr, file_space_id = space_id, & - mem_space_id = memspace, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5dwrite_f instance_id') +!--------------------------------------------------------------------------------------------------- +! expand phaseAt to consider IPs (is not stored per IP) + do i = 1, size(homogenizationAt_perIP,1) + homogenizationAt_perIP(i,:) = homogenizationAt + enddo + +!--------------------------------------------------------------------------------------------------- +! renumber member from my process to all processes + do i = 1, size(label) + where(homogenizationAt_perIP == i) memberAt_total = memberAt + sum(memberOffset(i,0:worldrank-1)) + enddo !-------------------------------------------------------------------------------------------------- -!close types, dataspaces -call h5tclose_f(dtype_id, hdferr) -if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tclose_f dtype_id') -call h5tclose_f(position_id, hdferr) -if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tclose_f position_id') -call h5tclose_f(name_id, hdferr) -if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tclose_f name_id ') -call h5tclose_f(dt5_id, hdferr) -if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tclose_f dt5_id') -call h5dclose_f(dset_id, hdferr) -if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5dclose_f') -call h5sclose_f(space_id, hdferr) -if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5sclose_f space_id') -call h5sclose_f(memspace, hdferr) -if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5sclose_f memspace') -call h5pclose_f(plist_id, hdferr) -if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5pclose_f') -call HDF5_closeGroup(mapping_ID) - -end subroutine HDF5_mappingHomog +! write the components of the compound type individually + call h5pset_preserve_f(plist_id, .TRUE., ierr) + + loc_id = results_openGroup('/mapping/cellResults') + call h5dcreate_f(loc_id, 'materialpoint', dtype_id, filespace_id, dset_id, ierr) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dcreate_f') + + call h5dwrite_f(dset_id, name_id, reshape(label(pack(homogenizationAt_perIP,.true.)),myShape), & + myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dwrite_f/name_id') + call h5dwrite_f(dset_id, position_id, reshape(pack(memberAt_total,.true.),myShape), & + myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dwrite_f/position_id') !-------------------------------------------------------------------------------------------------- -!> @brief adds the backward mapping from spatial position and constituent ID to results -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_backwardMappingHomog(material_homog,homogmemberat,homogenization_name,dataspace_size,mpiOffset,mpiOffset_homog) - use hdf5 - - implicit none - integer(pInt), intent(in), dimension(:,:) :: material_homog, homogmemberat - character(len=*), intent(in), dimension(:) :: homogenization_name - integer(pInt), intent(in), dimension(:) :: dataspace_size, mpiOffset_homog - integer(pInt), intent(in) :: mpiOffset - - integer(pInt) :: hdferr, NmatPoints, i - integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id, position_id, plist_id, memspace - integer(SIZE_T) :: type_size - - integer(pInt), dimension(:,:), allocatable :: arr - - integer(HSIZE_T), dimension(1) :: counter - integer(HSSIZE_T), dimension(1) :: fileOffset - - character(len=64) :: homogID - - NmatPoints = count(material_homog /=0_pInt) - allocate(arr(2,NmatPoints)) - - arr(1,:) = (/(i, i=0_pint,NmatPoints-1_pInt)/) - arr(2,:) = pack(material_homog,material_homog/=0_pInt) - - do i=1_pInt, size(homogenization_name) - write(homogID, '(i0)') i - mapping_ID = results_openGroup('/current/homogenization/'//trim(homogID)//'_'//homogenization_name(i)) - -!-------------------------------------------------------------------------------------------------- - ! create dataspace - call h5screate_simple_f(1, int([dataspace_size(i)],HSIZE_T), space_id, hdferr, & - int([dataspace_size(i)],HSIZE_T)) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeBackwardMapping') - -!-------------------------------------------------------------------------------------------------- - ! compound type - call h5tget_size_f(H5T_STD_I32LE, type_size, hdferr) - call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeBackwardMapping: h5tcreate_f dtype_id') - - call h5tinsert_f(dtype_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5tinsert_f 0') - -!-------------------------------------------------------------------------------------------------- - ! create Dataset - call h5dcreate_f(mapping_id, 'mapGeometry', dtype_id, space_id, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog') - -!-------------------------------------------------------------------------------------------------- - ! Create memory types (one compound datatype for each member) - call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5tcreate_f position_id') - call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5tinsert_f position_id') - -!-------------------------------------------------------------------------------------------------- - ! Define and select hyperslabs - counter = NmatPoints ! how big i am - fileOffset = mpiOffset_homog(i) ! where i start to write my data - - call h5screate_simple_f(1, counter, memspace, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5screate_simple_f') - call h5dget_space_f(dset_id, space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5dget_space_f') - call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5sselect_hyperslab_f') - -!-------------------------------------------------------------------------------------------------- - ! Create property list for collective dataset write -#ifdef PETSc - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5pcreate_f') - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5pset_dxpl_mpio_f') -#endif - -!-------------------------------------------------------------------------------------------------- - ! write data by fields in the datatype. Fields order is not important. - call h5dwrite_f(dset_id, position_id, pack(arr(1,:),arr(2,:)==i)+mpiOffset,int([dataspace_size(i)],HSIZE_T),& - hdferr, file_space_id = space_id, mem_space_id = memspace, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5dwrite_f instance_id') - -!-------------------------------------------------------------------------------------------------- - !close types, dataspaces - call h5tclose_f(dtype_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5tclose_f dtype_id') - call h5tclose_f(position_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5tclose_f position_id') - call h5dclose_f(dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5dclose_f') - call h5sclose_f(space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5sclose_f space_id') - call h5sclose_f(memspace, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5sclose_f memspace') - call h5pclose_f(plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5pclose_f') - call HDF5_closeGroup(mapping_ID) - - enddo - -end subroutine HDF5_backwardMappingHomog - -!-------------------------------------------------------------------------------------------------- -!> @brief adds the unique mapping from spatial position and constituent ID to results -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_mappingCrystallite(crystalliteAt,crystmemberAt,crystallite_name,dataspace_size,mpiOffset,mpiOffset_cryst) - use hdf5 - - implicit none - integer(pInt), intent(in), dimension(:,:) :: crystalliteAt - integer(pInt), intent(in), dimension(:,:,:) :: crystmemberAt - character(len=*), intent(in), dimension(:) :: crystallite_name - integer(pInt), intent(in), dimension(:) :: mpiOffset_cryst - integer(pInt), intent(in) :: dataspace_size, mpiOffset - - integer :: hdferr - integer(pInt) :: NmatPoints, Nconstituents, i, j - integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id, name_id, plist_id, memspace - - integer(HID_T), dimension(:), allocatable :: position_id - - integer(HID_T) :: dt5_id ! Memory datatype identifier - integer(SIZE_T) :: typesize, type_sizec, type_sizei, type_size - - integer(HSIZE_T), dimension(1) :: counter - integer(HSSIZE_T), dimension(1) :: fileOffset - integer(pInt), dimension(:), allocatable :: arrOffset - - character(len=64) :: m - - Nconstituents = size(crystmemberAt,1) - NmatPoints = count(crystalliteAt /=0_pInt) - mapping_ID = results_openGroup("current/mapGeometry") - - allocate(position_id(Nconstituents)) - - allocate(arrOffset(NmatPoints)) - do i=1_pInt, NmatPoints - do j=1_pInt, size(crystallite_name) - if(crystalliteAt(1,i) == j) & - arrOffset(i) = Nconstituents*mpiOffset_cryst(j) - enddo - enddo - -!-------------------------------------------------------------------------------------------------- -! create dataspace - call h5screate_simple_f(1, int([dataspace_size],HSIZE_T), space_id, hdferr, & - int([dataspace_size],HSIZE_T)) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeMapping') - -!-------------------------------------------------------------------------------------------------- -! compound type - ! First calculate total size by calculating sizes of each member - ! - CALL h5tcopy_f(H5T_NATIVE_CHARACTER, dt5_id, hdferr) - typesize = len(crystallite_name(1)) - CALL h5tset_size_f(dt5_id, typesize, hdferr) - CALL h5tget_size_f(dt5_id, type_sizec, hdferr) - CALL h5tget_size_f(H5T_STD_I32LE, type_sizei, hdferr) - type_size = type_sizec + type_sizei*Nconstituents - call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeMapping: h5tcreate_f dtype_id') - - call h5tinsert_f(dtype_id, "Name", 0_SIZE_T, dt5_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f 0') - do i=1_pInt, Nconstituents - write(m, '(i0)') i - call h5tinsert_f(dtype_id, "Position "//trim(m), type_sizec+(i-1)*type_sizei, H5T_STD_I32LE, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f 2 '//trim(m)) - enddo - -!-------------------------------------------------------------------------------------------------- -! create Dataset - call h5dcreate_f(mapping_id, 'crystallite', dtype_id, space_id, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite') - -!-------------------------------------------------------------------------------------------------- -! Create memory types (one compound datatype for each member) - call h5tcreate_f(H5T_COMPOUND_F, int(type_sizec,SIZE_T), name_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f instance_id') - call h5tinsert_f(name_id, "Name", 0_SIZE_T, dt5_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f instance_id') - - do i=1_pInt, Nconstituents - write(m, '(i0)') i - call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id(i), hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f position_id') - call h5tinsert_f(position_id(i), "Position "//trim(m), 0_SIZE_T, H5T_STD_I32LE, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f position_id') - enddo - -!-------------------------------------------------------------------------------------------------- -! Define and select hyperslabs - counter = NmatPoints ! how big i am - fileOffset = mpiOffset ! where i start to write my data - - call h5screate_simple_f(1, counter, memspace, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5screate_simple_f') - call h5dget_space_f(dset_id, space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dget_space_f') - call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5sselect_hyperslab_f') - -!-------------------------------------------------------------------------------------------------- - ! Create property list for collective dataset write -#ifdef PETSc - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5pcreate_f') - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5pset_dxpl_mpio_f') -#endif - -!-------------------------------------------------------------------------------------------------- -! write data by fields in the datatype. Fields order is not important. - call h5dwrite_f(dset_id, name_id, crystallite_name(pack(crystalliteAt,crystalliteAt/=0_pInt)), & - int([dataspace_size],HSIZE_T), hdferr, file_space_id = space_id, & - mem_space_id = memspace, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dwrite_f position_id') - - do i=1_pInt, Nconstituents - call h5dwrite_f(dset_id, position_id(i), pack(crystmemberAt(i,:,:)-1_pInt,crystmemberAt(i,:,:)/=0_pInt)+arrOffset,& - int([dataspace_size],HSIZE_T), hdferr, file_space_id = space_id, & - mem_space_id = memspace, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dwrite_f instance_id') - enddo - -!-------------------------------------------------------------------------------------------------- -!close types, dataspaces - call h5tclose_f(dtype_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f dtype_id') - do i=1_pInt, Nconstituents - call h5tclose_f(position_id(i), hdferr) - enddo - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f position_id') - call h5tclose_f(name_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f name_id') - call h5tclose_f(dt5_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f dt5_id') - call h5dclose_f(dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dclose_f') - call h5sclose_f(space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5sclose_f space_id') - call h5sclose_f(memspace, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5sclose_f memspace') - call h5pclose_f(plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5pclose_f') - call HDF5_closeGroup(mapping_ID) - -end subroutine HDF5_mappingCrystallite - - -!-------------------------------------------------------------------------------------------------- -!> @brief adds the backward mapping from spatial position and constituent ID to results -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_backwardMappingCrystallite(crystalliteAt,crystmemberAt,crystallite_name,dataspace_size,mpiOffset,mpiOffset_cryst) - use hdf5 - - implicit none - integer(pInt), intent(in), dimension(:,:) :: crystalliteAt - integer(pInt), intent(in), dimension(:,:,:) :: crystmemberAt - character(len=*), intent(in), dimension(:) :: crystallite_name - integer(pInt), intent(in), dimension(:) :: dataspace_size, mpiOffset_cryst - integer(pInt), intent(in) :: mpiOffset - - integer :: hdferr - integer(pInt) :: NmatPoints, Nconstituents, i, j - integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id, position_id, plist_id, memspace - integer(SIZE_T) :: type_size - - integer(pInt), dimension(:,:), allocatable :: h_arr, arr - - integer(HSIZE_T), dimension(1) :: counter - integer(HSSIZE_T), dimension(1) :: fileOffset - - character(len=64) :: crystallID - - Nconstituents = size(crystmemberAt,1) - NmatPoints = count(crystalliteAt /=0_pInt) - - allocate(h_arr(2,NmatPoints)) - allocate(arr(2,Nconstituents*NmatPoints)) - - h_arr(1,:) = (/(i, i=0_pInt,NmatPoints-1_pInt)/) - h_arr(2,:) = pack(crystalliteAt,crystalliteAt/=0_pInt) - - do i=1_pInt, NmatPoints - do j=Nconstituents-1_pInt, 0_pInt, -1_pInt - arr(1,Nconstituents*i-j) = h_arr(1,i) - arr(2,Nconstituents*i-j) = h_arr(2,i) - enddo - enddo - - do i=1_pInt, size(crystallite_name) - if (crystallite_name(i) == 'none') cycle - write(crystallID, '(i0)') i - mapping_ID = results_openGroup('/current/crystallite/'//trim(crystallID)//'_'//crystallite_name(i)) - NmatPoints = count(crystalliteAt == i) - -!-------------------------------------------------------------------------------------------------- - ! create dataspace - call h5screate_simple_f(1, int([Nconstituents*dataspace_size(i)],HSIZE_T), space_id, hdferr, & - int([Nconstituents*dataspace_size(i)],HSIZE_T)) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeBackwardMapping') - -!-------------------------------------------------------------------------------------------------- - ! compound type - call h5tget_size_f(H5T_STD_I32LE, type_size, hdferr) - call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeBackwardMapping: h5tcreate_f dtype_id') - - call h5tinsert_f(dtype_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5tinsert_f 0') - -!-------------------------------------------------------------------------------------------------- - ! create Dataset - call h5dcreate_f(mapping_id, 'mapGeometry', dtype_id, space_id, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite') - -!-------------------------------------------------------------------------------------------------- - ! Create memory types - call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5tcreate_f position_id') - call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5tinsert_f position_id') - -!-------------------------------------------------------------------------------------------------- - ! Define and select hyperslabs - counter = Nconstituents*NmatPoints ! how big i am - fileOffset = Nconstituents*mpiOffset_cryst(i) ! where i start to write my data - - call h5screate_simple_f(1, counter, memspace, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5screate_simple_f') - call h5dget_space_f(dset_id, space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5dget_space_f') - call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5sselect_hyperslab_f') - -!-------------------------------------------------------------------------------------------------- - ! Create property list for collective dataset write -#ifdef PETSc - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5pcreate_f') - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5pset_dxpl_mpio_f') -#endif - -!-------------------------------------------------------------------------------------------------- - ! write data by fields in the datatype. Fields order is not important. - call h5dwrite_f(dset_id, position_id, pack(arr(1,:),arr(2,:)==i) + mpiOffset,& - int([Nconstituents*dataspace_size(i)],HSIZE_T), hdferr, file_space_id = space_id, & - mem_space_id = memspace, xfer_prp = plist_id) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5dwrite_f instance_id') - -!-------------------------------------------------------------------------------------------------- - !close types, dataspaces - call h5tclose_f(dtype_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5tclose_f dtype_id') - call h5tclose_f(position_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5tclose_f position_id') - call h5dclose_f(dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5dclose_f') - call h5sclose_f(space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5sclose_f space_id') - call h5sclose_f(memspace, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5sclose_f memspace') - call h5pclose_f(plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5pclose_f') - call HDF5_closeGroup(mapping_ID) - - enddo - -end subroutine HDF5_backwardMappingCrystallite - -!-------------------------------------------------------------------------------------------------- -!> @brief adds the unique cell to node mapping -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_mappingCells(mapping) - use hdf5 - - implicit none - integer(pInt), intent(in), dimension(:) :: mapping - - integer :: hdferr, Nnodes - integer(HID_T) :: mapping_id, dset_id, space_id - - Nnodes=size(mapping) - mapping_ID = results_openGroup("mapping") - -!-------------------------------------------------------------------------------------------------- -! create dataspace - call h5screate_simple_f(1, int([Nnodes],HSIZE_T), space_id, hdferr, & - int([Nnodes],HSIZE_T)) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells: h5screate_simple_f') - -!-------------------------------------------------------------------------------------------------- -! create Dataset - call h5dcreate_f(mapping_id, "Cell",H5T_NATIVE_INTEGER, space_id, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells') - -!-------------------------------------------------------------------------------------------------- -! write data by fields in the datatype. Fields order is not important. - call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, mapping, int([Nnodes],HSIZE_T), hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells: h5dwrite_f instance_id') - -!-------------------------------------------------------------------------------------------------- -!close types, dataspaces - call h5dclose_f(dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dclose_f') - call h5sclose_f(space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5sclose_f') - call HDF5_closeGroup(mapping_ID) - -end subroutine HDF5_mappingCells +! close all + call HDF5_closeGroup(loc_id) + call h5pclose_f(plist_id, ierr) + call h5sclose_f(filespace_id, ierr) + call h5sclose_f(memspace_id, ierr) + call h5dclose_f(dset_id, ierr) + call h5tclose_f(dtype_id, ierr) + call h5tclose_f(name_id, ierr) + call h5tclose_f(position_id, ierr) + +end subroutine results_mapping_materialpoint + + +!!-------------------------------------------------------------------------------------------------- +!!> @brief adds the backward mapping from spatial position and constituent ID to results +!!-------------------------------------------------------------------------------------------------- +!subroutine HDF5_backwardMappingPhase(material_phase,phasememberat,phase_name,dataspace_size,mpiOffset,mpiOffset_phase) +! use hdf5 + +! integer(pInt), intent(in), dimension(:,:,:) :: material_phase, phasememberat +! character(len=*), intent(in), dimension(:) :: phase_name +! integer(pInt), intent(in), dimension(:) :: dataspace_size, mpiOffset_phase +! integer(pInt), intent(in) :: mpiOffset + +! integer(pInt) :: hdferr, NmatPoints, Nconstituents, i, j +! integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id, position_id, plist_id, memspace +! integer(SIZE_T) :: type_size + +! integer(pInt), dimension(:,:), allocatable :: arr + +! integer(HSIZE_T), dimension(1) :: counter +! integer(HSSIZE_T), dimension(1) :: fileOffset + +! character(len=64) :: phaseID + +! Nconstituents = size(phasememberat,1) +! NmatPoints = count(material_phase /=0)/Nconstituents + +! allocate(arr(2,NmatPoints*Nconstituents)) + +! do i=1, NmatPoints +! do j=Nconstituents-1, 0, -1 +! arr(1,Nconstituents*i-j) = i-1 +! enddo +! enddo +! arr(2,:) = pack(material_phase,material_phase/=0) + +! do i=1, size(phase_name) +! write(phaseID, '(i0)') i +! mapping_ID = results_openGroup('/current/constitutive/'//trim(phaseID)//'_'//phase_name(i)) +! NmatPoints = count(material_phase == i) + +!!-------------------------------------------------------------------------------------------------- +! ! create dataspace +! call h5screate_simple_f(1, int([dataspace_size(i)],HSIZE_T), space_id, hdferr, & +! int([dataspace_size(i)],HSIZE_T)) +! if (hdferr < 0) call IO_error(1,ext_msg='HDF5_writeBackwardMapping') + +!!-------------------------------------------------------------------------------------------------- +! ! compound type +! call h5tget_size_f(H5T_STD_I32LE, type_size, hdferr) +! call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='HDF5_writeBackwardMapping: h5tcreate_f dtype_id') + +! call h5tinsert_f(dtype_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5tinsert_f 0') + +!!-------------------------------------------------------------------------------------------------- +! ! create Dataset +! call h5dcreate_f(mapping_id, 'mapGeometry', dtype_id, space_id, dset_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase') + +!!-------------------------------------------------------------------------------------------------- +! ! Create memory types (one compound datatype for each member) +! call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5tcreate_f position_id') +! call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5tinsert_f position_id') + +!!-------------------------------------------------------------------------------------------------- +! ! Define and select hyperslabs +! counter = NmatPoints ! how big i am +! fileOffset = mpiOffset_phase(i) ! where i start to write my data + +! call h5screate_simple_f(1, counter, memspace, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5screate_simple_f') +! call h5dget_space_f(dset_id, space_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5dget_space_f') +! call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5sselect_hyperslab_f') + +!!-------------------------------------------------------------------------------------------------- +! ! Create property list for collective dataset write +!#ifdef PETSc +! call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5pcreate_f') +! call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5pset_dxpl_mpio_f') +!#endif + +!!-------------------------------------------------------------------------------------------------- +! ! write data by fields in the datatype. Fields order is not important. +! call h5dwrite_f(dset_id, position_id, pack(arr(1,:),arr(2,:)==i)+mpiOffset, int([dataspace_size(i)],HSIZE_T),& +! hdferr, file_space_id = space_id, mem_space_id = memspace, xfer_prp = plist_id) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5dwrite_f instance_id') + +!!-------------------------------------------------------------------------------------------------- +! !close types, dataspaces +! call h5tclose_f(dtype_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5tclose_f dtype_id') +! call h5tclose_f(position_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5tclose_f position_id') +! call h5dclose_f(dset_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5dclose_f') +! call h5sclose_f(space_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5sclose_f space_id') +! call h5sclose_f(memspace, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5sclose_f memspace') +! call h5pclose_f(plist_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingPhase: h5pclose_f') +! call HDF5_closeGroup(mapping_ID) + +! enddo + +!end subroutine HDF5_backwardMappingPhase + + +!!-------------------------------------------------------------------------------------------------- +!!> @brief adds the backward mapping from spatial position and constituent ID to results +!!-------------------------------------------------------------------------------------------------- +!subroutine HDF5_backwardMappingHomog(material_homog,homogmemberat,homogenization_name,dataspace_size,mpiOffset,mpiOffset_homog) +! use hdf5 + +! integer(pInt), intent(in), dimension(:,:) :: material_homog, homogmemberat +! character(len=*), intent(in), dimension(:) :: homogenization_name +! integer(pInt), intent(in), dimension(:) :: dataspace_size, mpiOffset_homog +! integer(pInt), intent(in) :: mpiOffset + +! integer(pInt) :: hdferr, NmatPoints, i +! integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id, position_id, plist_id, memspace +! integer(SIZE_T) :: type_size + +! integer(pInt), dimension(:,:), allocatable :: arr + +! integer(HSIZE_T), dimension(1) :: counter +! integer(HSSIZE_T), dimension(1) :: fileOffset + +! character(len=64) :: homogID + +! NmatPoints = count(material_homog /=0) +! allocate(arr(2,NmatPoints)) + +! arr(1,:) = (/(i, i=0,NmatPoints-1)/) +! arr(2,:) = pack(material_homog,material_homog/=0) + +! do i=1, size(homogenization_name) +! write(homogID, '(i0)') i +! mapping_ID = results_openGroup('/current/homogenization/'//trim(homogID)//'_'//homogenization_name(i)) + +!!-------------------------------------------------------------------------------------------------- +! ! create dataspace +! call h5screate_simple_f(1, int([dataspace_size(i)],HSIZE_T), space_id, hdferr, & +! int([dataspace_size(i)],HSIZE_T)) +! if (hdferr < 0) call IO_error(1,ext_msg='HDF5_writeBackwardMapping') + +!!-------------------------------------------------------------------------------------------------- +! ! compound type +! call h5tget_size_f(H5T_STD_I32LE, type_size, hdferr) +! call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='HDF5_writeBackwardMapping: h5tcreate_f dtype_id') + +! call h5tinsert_f(dtype_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5tinsert_f 0') + +!!-------------------------------------------------------------------------------------------------- +! ! create Dataset +! call h5dcreate_f(mapping_id, 'mapGeometry', dtype_id, space_id, dset_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog') + +!!-------------------------------------------------------------------------------------------------- +! ! Create memory types (one compound datatype for each member) +! call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5tcreate_f position_id') +! call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5tinsert_f position_id') + +!!-------------------------------------------------------------------------------------------------- +! ! Define and select hyperslabs +! counter = NmatPoints ! how big i am +! fileOffset = mpiOffset_homog(i) ! where i start to write my data + +! call h5screate_simple_f(1, counter, memspace, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5screate_simple_f') +! call h5dget_space_f(dset_id, space_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5dget_space_f') +! call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5sselect_hyperslab_f') + +!!-------------------------------------------------------------------------------------------------- +! ! Create property list for collective dataset write +!#ifdef PETSc +! call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5pcreate_f') +! call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5pset_dxpl_mpio_f') +!#endif + +!!-------------------------------------------------------------------------------------------------- +! ! write data by fields in the datatype. Fields order is not important. +! call h5dwrite_f(dset_id, position_id, pack(arr(1,:),arr(2,:)==i)+mpiOffset,int([dataspace_size(i)],HSIZE_T),& +! hdferr, file_space_id = space_id, mem_space_id = memspace, xfer_prp = plist_id) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5dwrite_f instance_id') + +!!-------------------------------------------------------------------------------------------------- +! !close types, dataspaces +! call h5tclose_f(dtype_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5tclose_f dtype_id') +! call h5tclose_f(position_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5tclose_f position_id') +! call h5dclose_f(dset_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5dclose_f') +! call h5sclose_f(space_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5sclose_f space_id') +! call h5sclose_f(memspace, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5sclose_f memspace') +! call h5pclose_f(plist_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_backwardMappingHomog: h5pclose_f') +! call HDF5_closeGroup(mapping_ID) + +! enddo + +!end subroutine HDF5_backwardMappingHomog + + +!!-------------------------------------------------------------------------------------------------- +!!> @brief adds the unique cell to node mapping +!!-------------------------------------------------------------------------------------------------- +!subroutine HDF5_mappingCells(mapping) +! use hdf5 + +! integer(pInt), intent(in), dimension(:) :: mapping + +! integer :: hdferr, Nnodes +! integer(HID_T) :: mapping_id, dset_id, space_id + +! Nnodes=size(mapping) +! mapping_ID = results_openGroup("mapping") + +!!-------------------------------------------------------------------------------------------------- +!! create dataspace +! call h5screate_simple_f(1, int([Nnodes],HSIZE_T), space_id, hdferr, & +! int([Nnodes],HSIZE_T)) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_mappingCells: h5screate_simple_f') + +!!-------------------------------------------------------------------------------------------------- +!! create Dataset +! call h5dcreate_f(mapping_id, "Cell",H5T_NATIVE_INTEGER, space_id, dset_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_mappingCells') + +!!-------------------------------------------------------------------------------------------------- +!! write data by fields in the datatype. Fields order is not important. +! call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, mapping, int([Nnodes],HSIZE_T), hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_mappingCells: h5dwrite_f instance_id') + +!!-------------------------------------------------------------------------------------------------- +!!close types, dataspaces +! call h5dclose_f(dset_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_mappingConstitutive: h5dclose_f') +! call h5sclose_f(space_id, hdferr) +! if (hdferr < 0) call IO_error(1,ext_msg='IO_mappingConstitutive: h5sclose_f') +! call HDF5_closeGroup(mapping_ID) + +!end subroutine HDF5_mappingCells end module results diff --git a/src/rotations.f90 b/src/rotations.f90 index b899adacb..7c29b03bc 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -777,12 +777,12 @@ pure function qu2ax(qu) result(ax) real(pReal) :: omega, s - omega = 2.0 * acos(math_clip(qu%w,-1.0_pReal,1.0_pReal)) ! if the angle equals zero, then we return the rotation axis as [001] - if (dEq0(omega)) then - ax = [ 0.0, 0.0, 1.0, 0.0 ] + if (dEq0(sqrt(qu%x**2+qu%y**2+qu%z**2))) then + ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] elseif (dNeq0(qu%w)) then s = sign(1.0_pReal,qu%w)/sqrt(qu%x**2+qu%y**2+qu%z**2) + omega = 2.0_pReal * acos(math_clip(qu%w,-1.0_pReal,1.0_pReal)) ax = [ qu%x*s, qu%y*s, qu%z*s, omega ] else ax = [ qu%x, qu%y, qu%z, PI ]